From dcd62e8e7943d1327daeb16d90751ec1046ff90d Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@uliege.be>
Date: Mon, 3 Jan 2022 10:22:08 +0100
Subject: [PATCH] empty initializers are not allowed in ANSI C (#1695)
---
CHANGELOG.txt | 8 +-
contrib/MathEx/mathex.cpp | 2 +-
doc/texinfo/gmsh.texi | 7 +-
doc/texinfo/opt_fields.texi | 41 ++++-
doc/texinfo/opt_mesh.texi | 2 +-
examples/boolean/extend_field.geo | 38 +++++
src/common/DefaultOptions.h | 4 +-
src/mesh/BackgroundMeshTools.cpp | 6 +-
src/mesh/Field.cpp | 272 +++++++++++++++++++++++-------
src/mesh/meshGRegionHxt.cpp | 2 +-
src/mesh/meshGRegionMMG.cpp | 6 +-
src/parser/Gmsh.tab.cpp | 4 +-
src/parser/Gmsh.y | 4 +-
tutorials/c++/t10.cpp | 6 +-
tutorials/c/t2.c | 4 +-
tutorials/python/t10.py | 5 +-
tutorials/t10.geo | 7 +-
17 files changed, 320 insertions(+), 98 deletions(-)
create mode 100644 examples/boolean/extend_field.geo
diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index 50b21bd193..136b6eac95 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,7 +1,9 @@
4.9.3 (Work-in-progress): improved handling of degenerate curves in periodic
-surfaces and boundary layer extrusion; moved all kernel sources to src/
-subdirectory; renamed demos/ as examples/ and tutorial/ as tutorials/; extended
-Mesh.SaveGroupsOfElements capabilities for INP export; small bug fixes.
+surfaces and boundary layer extrusion; extended Mesh.SaveGroupsOfElements
+capabilities for INP export; extended Mesh.MeshSizeExtendFromBoundary + new
+"Extend" mesh size field to enable alternative mesh size extensions from
+boundary; moved all kernel sources to src/ subdirectory; renamed demos/ as
+examples/ and tutorial/ as tutorials/; small bug fixes.
4.9.2 (December 23, 2021): faster projection on OCC entities; extended
Mesh.SaveGroupsOfNodes capabilities for INP export; improved transfinite meshing
diff --git a/contrib/MathEx/mathex.cpp b/contrib/MathEx/mathex.cpp
index 81a609d825..76522c3553 100644
--- a/contrib/MathEx/mathex.cpp
+++ b/contrib/MathEx/mathex.cpp
@@ -182,7 +182,7 @@
// in future, put more precisery checking, for overflow ?
{
if (y == 0)
- throw mathex::error("Error [binary_divide()]: divisin by zero");
+ throw mathex::error("Error [binary_divide()]: division by zero");
else
return (x/y);
} // binary_divide()
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index 981778d01c..52b17375bf 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -2882,10 +2882,9 @@ the smallest element size is selected at any given point. Using the
Gmsh API, the resulting value can be further modified using a C++, C,
Python or Julia mesh size callback function provided via
@code{gmsh/model/mesh/setSizeCallback} (@pxref{Namespace
-gmsh/model/mesh}). In addition, boundary mesh sizes (on curves or
-surfaces) are interpolated inside the enclosed entity (surface or
-volume, respectively) if the option
-@code{Mesh.MeshSizeExtendFromBoundary} is set (it is by default).
+gmsh/model/mesh}). In addition, boundary mesh sizes are interpolated
+inside surfaces and/or volumes depending on the value of
+@code{Mesh.MeshSizeExtendFromBoundary}.
All element sizes are further constrained in the interval [
@code{Mesh.MeshSizeMin}, @code{Mesh.MeshSizeMax} ] (which can also be
diff --git a/doc/texinfo/opt_fields.texi b/doc/texinfo/opt_fields.texi
index d6580e73cb..4100c7e692 100644
--- a/doc/texinfo/opt_fields.texi
+++ b/doc/texinfo/opt_fields.texi
@@ -356,6 +356,37 @@ type: list@*
default value: @code{@{@}}
@end table
+@item Extend
+Compute an extension of the mesh sizes from the given boundary curves (resp. surfaces) inside the surfaces (resp. volumes) being meshed. If the mesh size on the boundary, computed as the local average length of the edges of the boundary elements, is denoted by SizeBnd, the extension is computed as:@*
+@*
+ F = f * SizeBnd + (1 - f) * SizeMax, if d < DistMax@*
+@*
+ F = SizeMax, if d >= DistMax@*
+where d denotes the distance to the boundary and f = ((DistMax - d) / DistMax)^Power.@*
+Options:@*
+@table @code
+@item CurvesList
+Tags of curves in the geometric model@*
+type: list@*
+default value: @code{@{@}}
+@item DistMax
+Maximum distance of the size extension@*
+type: float@*
+default value: @code{1}
+@item Power
+Power exponent used to interpolate the mesh size@*
+type: float@*
+default value: @code{1}
+@item SizeMax
+Mesh size outside DistMax@*
+type: float@*
+default value: @code{1}
+@item SurfacesList
+Tags of surfaces in the geometric model@*
+type: list@*
+default value: @code{@{@}}
+@end table
+
@item ExternalProcess
**This Field is experimental**@*
Call an external process that received coordinates triple (x,y,z) as binary double precision numbers on stdin and is supposed to write the field value on stdout as a binary double precision number.@*
@@ -800,15 +831,15 @@ Return F = SizeMin if Field[InField] <= DistMin, F = SizeMax if Field[InField] >
Options:@*
@table @code
@item DistMax
-Distance from entity after which the mesh size will be SizeMax@*
+Value after which the mesh size will be SizeMax@*
type: float@*
default value: @code{10}
@item DistMin
-Distance from entity up to which the mesh size will be SizeMin@*
+Value up to which the mesh size will be SizeMin@*
type: float@*
default value: @code{1}
@item InField
-Tag of the field to evaluate@*
+Tag of the field computing the input value, usually a distance@*
type: integer@*
default value: @code{0}
@item Sigmoid
@@ -816,11 +847,11 @@ True to interpolate between SizeMin and LcMax using a sigmoid, false to interpol
type: boolean@*
default value: @code{0}
@item SizeMax
-Mesh size outside DistMax@*
+Mesh size when value > DistMax@*
type: float@*
default value: @code{1}
@item SizeMin
-Mesh size inside DistMin@*
+Mesh size when value < DistMin@*
type: float@*
default value: @code{0.1}
@item StopAtDistMax
diff --git a/doc/texinfo/opt_mesh.texi b/doc/texinfo/opt_mesh.texi
index 5ed24b9588..25ca5c6525 100644
--- a/doc/texinfo/opt_mesh.texi
+++ b/doc/texinfo/opt_mesh.texi
@@ -327,7 +327,7 @@ Default value: @code{0}@*
Saved in: @code{General.OptionsFileName}
@item Mesh.MeshSizeExtendFromBoundary
-Extend computation of mesh element sizes from the boundaries into the interior (for 3D Delaunay, use 1: longest or 2: shortest surface edge length)@*
+Extend computation of mesh element sizes from the boundaries into the interior (0: never; 1: for surfaces and volumes; 2: for surfaces and volumes, but use smallest surface element edge length instead of longest length in 3D Delaunay; -2: only for surfaces; -3: only for volumes)@*
Default value: @code{1}@*
Saved in: @code{General.OptionsFileName}
diff --git a/examples/boolean/extend_field.geo b/examples/boolean/extend_field.geo
new file mode 100644
index 0000000000..80a3fe627d
--- /dev/null
+++ b/examples/boolean/extend_field.geo
@@ -0,0 +1,38 @@
+SetFactory("OpenCASCADE");
+Box(1) = {0,0,0, 1,1,0.5};
+Sphere(2) = {1,1,0.5,0.4};
+BooleanDifference(3) = { Volume{1}; Delete; }{ Volume{2}; Delete; };
+
+Box(10) = {0.3,0.3,0.5, 0.1,0.1,0.05};
+Box(11) = {0.5,0.5,0.5, 0.1,0.1,0.05};
+BooleanFragments{ Volume{10, 11}; Delete; }{ Volume{3}; Delete; }
+
+MeshSize{ : } = 0.04;
+//MeshSize{ : } = 0.002;
+MeshSize{ PointsOf{ Volume{10, 11}; } } = 0.002;
+Mesh.Algorithm3D = 10;
+
+DefineConstant[
+ use_field = {1, Choices{0, 1},
+ Name "Parameters/0Use Extend field?"}
+ size_max = {0.04, Min 0.002, Max 0.2, Step 0.001,
+ Name "Parameters/Mesh size after Extend", Visible use_field}
+ dist_max = {0.1, Min 0.01, Max 1, Step 0.01,
+ Name "Parameters/Extend distance", Visible use_field}
+ power = {2, Min 0.1, Max 10, Step 0.1,
+ Name "Parameters/Extend power", Visible use_field}
+];
+
+If(use_field == 0)
+ Mesh.MeshSizeExtendFromBoundary = 1;
+Else
+ // set to -2 or -3 to only extend in 2D or 3D
+ Mesh.MeshSizeExtendFromBoundary = 0;
+ Field[1] = Extend;
+ Field[1].SurfacesList = {Surface{:}};
+ Field[1].CurvesList = {Curve{:}};
+ Field[1].DistMax = dist_max;
+ Field[1].SizeMax = size_max;
+ Field[1].Power = power;
+ Background Field = 1;
+EndIf
diff --git a/src/common/DefaultOptions.h b/src/common/DefaultOptions.h
index 46c9a96413..30b443e7e9 100644
--- a/src/common/DefaultOptions.h
+++ b/src/common/DefaultOptions.h
@@ -1246,7 +1246,9 @@ StringXNumber MeshOptions_Number[] = {
"Mesh only entities that have no existing mesh" },
{ F|O, "MeshSizeExtendFromBoundary" , opt_mesh_lc_extend_from_boundary, 1. ,
"Extend computation of mesh element sizes from the boundaries into the interior "
- "(for 3D Delaunay, use 1: longest or 2: shortest surface edge length)"},
+ "(0: never; 1: for surfaces and volumes; 2: for surfaces and volumes, but use "
+ "smallest surface element edge length instead of longest length in 3D Delaunay; "
+ "-2: only for surfaces; -3: only for volumes)"},
{ F|O, "MeshSizeFactor" , opt_mesh_lc_factor , 1.0 ,
"Factor applied to all mesh element sizes" },
{ F|O, "MeshSizeMin" , opt_mesh_lc_min, 0.0 ,
diff --git a/src/mesh/BackgroundMeshTools.cpp b/src/mesh/BackgroundMeshTools.cpp
index 1d19d77f19..3f51f05117 100644
--- a/src/mesh/BackgroundMeshTools.cpp
+++ b/src/mesh/BackgroundMeshTools.cpp
@@ -338,12 +338,14 @@ SMetric3 BGM_MeshMetric(GEntity *ge, double U, double V, double X, double Y,
bool Extend1dMeshIn2dSurfaces(GFace *gf)
{
- return gf->getMeshSizeFromBoundary();
+ int val = gf->getMeshSizeFromBoundary();
+ return (val > 0 || val == -2);
}
bool Extend2dMeshIn3dVolumes()
{
- return CTX::instance()->mesh.lcExtendFromBoundary;
+ int val = CTX::instance()->mesh.lcExtendFromBoundary;
+ return (val > 0 || val == -3);
}
SMetric3 max_edge_curvature_metric(const GVertex *gv)
diff --git a/src/mesh/Field.cpp b/src/mesh/Field.cpp
index 04cd871a57..a5ca502ce3 100644
--- a/src/mesh/Field.cpp
+++ b/src/mesh/Field.cpp
@@ -633,15 +633,16 @@ public:
_stopAtDistMax = false;
options["InField"] =
- new FieldOptionInt(_inField, "Tag of the field to evaluate");
+ new FieldOptionInt(_inField, "Tag of the field computing the input "
+ "value, usually a distance");
options["DistMin"] = new FieldOptionDouble(
- _dMin, "Distance from entity up to which the mesh size will be SizeMin");
+ _dMin, "Value up to which the mesh size will be SizeMin");
options["DistMax"] = new FieldOptionDouble(
- _dMax, "Distance from entity after which the mesh size will be SizeMax");
+ _dMax, "Value after which the mesh size will be SizeMax");
options["SizeMin"] =
- new FieldOptionDouble(_lcMin, "Mesh size inside DistMin");
+ new FieldOptionDouble(_lcMin, "Mesh size when value < DistMin");
options["SizeMax"] =
- new FieldOptionDouble(_lcMax, "Mesh size outside DistMax");
+ new FieldOptionDouble(_lcMax, "Mesh size when value > DistMax");
options["Sigmoid"] = new FieldOptionBool(
_sigmoid,
"True to interpolate between SizeMin and LcMax using a sigmoid, "
@@ -2718,83 +2719,231 @@ public:
}
void update()
{
- _infos.clear();
- _pc.pts.clear();
- if(_kdtree) delete _kdtree;
+ if(updateNeeded) {
+ _infos.clear();
+ _pc.pts.clear();
+ if(_kdtree) delete _kdtree;
- for(auto it = _pointTags.begin(); it != _pointTags.end(); ++it) {
- GVertex *gv = GModel::current()->getVertexByTag(*it);
- if(gv) {
- _pc.pts.push_back(SPoint3(gv->x(), gv->y(), gv->z()));
- _infos.push_back(AttractorInfo(*it, 0, 0, 0));
- }
- else {
- Msg::Warning("Unknown point %d", *it);
+ for(auto it = _pointTags.begin(); it != _pointTags.end(); ++it) {
+ GVertex *gv = GModel::current()->getVertexByTag(*it);
+ if(gv) {
+ _pc.pts.push_back(SPoint3(gv->x(), gv->y(), gv->z()));
+ _infos.push_back(AttractorInfo(*it, 0, 0, 0));
+ }
+ else {
+ Msg::Warning("Unknown point %d", *it);
+ }
}
- }
- for(auto it = _curveTags.begin(); it != _curveTags.end(); ++it) {
- GEdge *e = GModel::current()->getEdgeByTag(*it);
- if(e) {
- if(e->mesh_vertices.size()) {
- for(std::size_t i = 0; i < e->mesh_vertices.size(); i++) {
- _pc.pts.push_back(SPoint3(e->mesh_vertices[i]->x(),
- e->mesh_vertices[i]->y(),
- e->mesh_vertices[i]->z()));
- double t = 0.;
- e->mesh_vertices[i]->getParameter(0, t);
+ for(auto it = _curveTags.begin(); it != _curveTags.end(); ++it) {
+ GEdge *e = GModel::current()->getEdgeByTag(*it);
+ if(e) {
+ if(e->mesh_vertices.size()) {
+ for(std::size_t i = 0; i < e->mesh_vertices.size(); i++) {
+ _pc.pts.push_back(SPoint3(e->mesh_vertices[i]->x(),
+ e->mesh_vertices[i]->y(),
+ e->mesh_vertices[i]->z()));
+ double t = 0.;
+ e->mesh_vertices[i]->getParameter(0, t);
+ _infos.push_back(AttractorInfo(*it, 1, t, 0));
+ }
+ }
+ int NNN = _sampling - e->mesh_vertices.size();
+ for(int i = 1; i < NNN - 1; i++) {
+ double u = (double)i / (NNN - 1);
+ Range<double> b = e->parBounds(0);
+ double t = b.low() + u * (b.high() - b.low());
+ GPoint gp = e->point(t);
+ _pc.pts.push_back(SPoint3(gp.x(), gp.y(), gp.z()));
_infos.push_back(AttractorInfo(*it, 1, t, 0));
}
}
- int NNN = _sampling - e->mesh_vertices.size();
- for(int i = 1; i < NNN - 1; i++) {
- double u = (double)i / (NNN - 1);
- Range<double> b = e->parBounds(0);
- double t = b.low() + u * (b.high() - b.low());
- GPoint gp = e->point(t);
- _pc.pts.push_back(SPoint3(gp.x(), gp.y(), gp.z()));
- _infos.push_back(AttractorInfo(*it, 1, t, 0));
+ else {
+ Msg::Warning("Unknown curve %d", *it);
}
}
- else {
- Msg::Warning("Unknown curve %d", *it);
+
+ for(auto it = _surfaceTags.begin(); it != _surfaceTags.end(); ++it) {
+ GFace *f = GModel::current()->getFaceByTag(*it);
+ if(f) {
+ double maxDist = f->bounds().diag() / _sampling;
+ std::vector<SPoint2> uvpoints;
+ f->fillPointCloud(maxDist, &_pc.pts, &uvpoints);
+ for(std::size_t i = 0; i < uvpoints.size(); i++)
+ _infos.push_back
+ (AttractorInfo(*it, 2, uvpoints[i].x(), uvpoints[i].y()));
+ }
+ else {
+ Msg::Warning("Unknown surface %d", *it);
+ }
}
+
+ // construct a kd-tree index:
+ _kdtree = new SPoint3KDTree(3, _pc2kdtree,
+ nanoflann::KDTreeSingleIndexAdaptorParams(10));
+ _kdtree->buildIndex();
+ updateNeeded = false;
}
+ }
+ using Field::operator();
+ virtual double operator()(double X, double Y, double Z, GEntity *ge = nullptr)
+ {
+ if(!_kdtree) return MAX_LC;
+ double pt[3] = {X, Y, Z};
+ nanoflann::KNNResultSet<double> res(1);
+ double outDistSqr;
+ res.init(&_outIndex, &outDistSqr);
+ _kdtree->findNeighbors(res, &pt[0], nanoflann::SearchParams(10));
+ return sqrt(outDistSqr);
+ }
+};
- for(auto it = _surfaceTags.begin(); it != _surfaceTags.end(); ++it) {
- GFace *f = GModel::current()->getFaceByTag(*it);
- if(f) {
- double maxDist = f->bounds().diag() / _sampling;
- std::vector<SPoint2> uvpoints;
- f->fillPointCloud(maxDist, &_pc.pts, &uvpoints);
- for(std::size_t i = 0; i < uvpoints.size(); i++)
- _infos.push_back
- (AttractorInfo(*it, 2, uvpoints[i].x(), uvpoints[i].y()));
+class ExtendField : public Field {
+ std::list<int> _tagCurves, _tagSurfaces;
+ std::vector<double> _sizeCurves, _sizeSurfaces;
+ SPoint3Cloud _pcCurves, _pcSurfaces;
+ SPoint3CloudAdaptor<SPoint3Cloud> _pc2kdtreeCurves, _pc2kdtreeSurfaces;
+ SPoint3KDTree *_kdtreeCurves, *_kdtreeSurfaces;
+ double _distMax, _sizeMax, _power;
+public:
+ ExtendField()
+ : _pc2kdtreeCurves(_pcCurves), _pc2kdtreeSurfaces(_pcSurfaces),
+ _kdtreeCurves(nullptr), _kdtreeSurfaces(nullptr)
+ {
+ _distMax = 1.;
+ _sizeMax = 1.;
+ _power = 1.;
+ options["CurvesList"] = new FieldOptionList(
+ _tagCurves, "Tags of curves in the geometric model", &updateNeeded);
+ options["SurfacesList"] = new FieldOptionList(
+ _tagSurfaces, "Tags of surfaces in the geometric model", &updateNeeded);
+ options["DistMax"] = new FieldOptionDouble(
+ _distMax, "Maximum distance of the size extension");
+ options["SizeMax"] = new FieldOptionDouble(
+ _sizeMax, "Mesh size outside DistMax");
+ options["Power"] = new FieldOptionDouble(
+ _power, "Power exponent used to interpolate the mesh size");
+ }
+ ~ExtendField()
+ {
+ if(_kdtreeCurves) delete _kdtreeCurves;
+ if(_kdtreeSurfaces) delete _kdtreeSurfaces;
+ }
+ const char *getName() { return "Extend"; }
+ std::string getDescription()
+ {
+ return "Compute an extension of the mesh sizes from the given boundary "
+ "curves (resp. surfaces) inside the surfaces (resp. volumes) "
+ "being meshed. If the mesh size on the boundary, computed "
+ "as the local average length of the edges of the boundary elements, "
+ "is denoted by SizeBnd, the extension is computed as:\n\n"
+ " F = f * SizeBnd + (1 - f) * SizeMax, if d < DistMax\n\n"
+ " F = SizeMax, if d >= DistMax\n"
+ "where d denotes the distance to the boundary and "
+ "f = ((DistMax - d) / DistMax)^Power.";
+ }
+ void recomputeCurves()
+ {
+ _sizeCurves.clear();
+ _pcCurves.pts.clear();
+ if(_kdtreeCurves) delete _kdtreeCurves;
+ for(auto t : _tagCurves) {
+ GEdge *ge = GModel::current()->getEdgeByTag(t);
+ if(ge) {
+ for(auto e : ge->lines) {
+ _pcCurves.pts.push_back(e->barycenter());
+ _sizeCurves.push_back(e->getEdge(0).length());
+ }
}
else {
- Msg::Warning("Unknown surface %d", *it);
+ Msg::Warning("Unknown curve %d", t);
}
}
-
- // construct a kd-tree index:
- _kdtree = new SPoint3KDTree(3, _pc2kdtree,
- nanoflann::KDTreeSingleIndexAdaptorParams(10));
- _kdtree->buildIndex();
- updateNeeded = false;
+ if(_sizeCurves.empty()) {
+ _kdtreeCurves = nullptr;
+ }
+ else {
+ _kdtreeCurves = new SPoint3KDTree(3, _pc2kdtreeCurves,
+ nanoflann::KDTreeSingleIndexAdaptorParams(10));
+ _kdtreeCurves->buildIndex();
+ }
+ }
+ void recomputeSurfaces()
+ {
+ _sizeSurfaces.clear();
+ _pcSurfaces.pts.clear();
+ if(_kdtreeSurfaces) delete _kdtreeSurfaces;
+ for(auto t : _tagSurfaces) {
+ GFace *gf = GModel::current()->getFaceByTag(t);
+ if(gf) {
+ for(auto e : gf->triangles) {
+ _pcSurfaces.pts.push_back(e->barycenter());
+ double s = e->getEdge(0).length() + e->getEdge(1).length() +
+ e->getEdge(2).length();
+ _sizeSurfaces.push_back(s / 3.);
+ }
+ }
+ else {
+ Msg::Warning("Unknown surface %d", t);
+ }
+ }
+ if(_sizeSurfaces.empty()) {
+ _kdtreeSurfaces = nullptr;
+ }
+ else {
+ _kdtreeSurfaces = new SPoint3KDTree(3, _pc2kdtreeSurfaces,
+ nanoflann::KDTreeSingleIndexAdaptorParams(10));
+ _kdtreeSurfaces->buildIndex();
+ }
}
using Field::operator();
virtual double operator()(double X, double Y, double Z, GEntity *ge = nullptr)
{
+ if(!ge) return MAX_LC;
+ if(ge->dim() != 2 && ge->dim() != 3) return MAX_LC;
+ if(ge->dim() == 2 && _tagCurves.empty()) return MAX_LC;
+ if(ge->dim() == 3 && _tagSurfaces.empty()) return MAX_LC;
#pragma omp critical
- {
- if(!_kdtree || updateNeeded) update();
+ if(updateNeeded ||
+ (ge->dim() == 2 && _tagCurves.size() && _sizeCurves.empty())) {
+ // we are meshing our first surface; recompute distance to the elements on
+ // curves, and invalidate the distance to surfaces
+ recomputeCurves();
+ _sizeSurfaces.clear();
+ updateNeeded = false;
}
- double query_pt[3] = {X, Y, Z};
- nanoflann::KNNResultSet<double> resultSet(1);
- double outDistSqr;
- resultSet.init(&_outIndex, &outDistSqr);
- _kdtree->findNeighbors(resultSet, &query_pt[0], nanoflann::SearchParams(10));
- return sqrt(outDistSqr);
+#pragma omp critical
+ if(updateNeeded ||
+ (ge->dim() == 3 && _tagSurfaces.size() && _sizeSurfaces.empty())) {
+ // we are meshing our first volume; recompute distance to the elements on
+ // surfaces, and invalidate the distance to curves (to be ready for
+ // subsequent surface meshing pass)
+ recomputeSurfaces();
+ _sizeCurves.clear();
+ updateNeeded = false;
+ }
+ double pt[3] = {X, Y, Z};
+ nanoflann::KNNResultSet<double> res(1);
+ std::size_t index = 0;
+ double dist2 = 0., sbnd = 0.;
+ res.init(&index, &dist2);
+ if(ge->dim() == 2 && _kdtreeCurves) {
+ _kdtreeCurves->findNeighbors(res, &pt[0], nanoflann::SearchParams(10));
+ sbnd = _sizeCurves[index];
+ }
+ else if(ge->dim() == 3 && _kdtreeSurfaces) {
+ _kdtreeSurfaces->findNeighbors(res, &pt[0], nanoflann::SearchParams(10));
+ sbnd = _sizeSurfaces[index];
+ }
+ if(sbnd <= 0 || _distMax <= 0 || _sizeMax <= 0) return MAX_LC;
+ double d = sqrt(dist2);
+ if(d >= _distMax) return _sizeMax;
+ double f = (_distMax - d) / _distMax;
+ if(_power != 1.) {
+ f = std::pow(f, _power);
+ }
+ double s = f * sbnd + (1. - f) * _sizeMax;
+ return s;
}
};
@@ -3225,6 +3374,7 @@ FieldManager::FieldManager()
mapTypeName["Gradient"] = new FieldFactoryT<GradientField>();
mapTypeName["Octree"] = new FieldFactoryT<OctreeField>();
mapTypeName["Distance"] = new FieldFactoryT<DistanceField>();
+ mapTypeName["Extend"] = new FieldFactoryT<ExtendField>();
mapTypeName["Restrict"] = new FieldFactoryT<RestrictField>();
mapTypeName["Constant"] = new FieldFactoryT<ConstantField>();
mapTypeName["Min"] = new FieldFactoryT<MinField>();
diff --git a/src/mesh/meshGRegionHxt.cpp b/src/mesh/meshGRegionHxt.cpp
index b4d6bda973..ab87c8a36c 100644
--- a/src/mesh/meshGRegionHxt.cpp
+++ b/src/mesh/meshGRegionHxt.cpp
@@ -41,7 +41,7 @@ static HXTStatus nodalSizesCallBack(double *pts, uint32_t *volume,
std::vector<GRegion *> *allGR = (std::vector<GRegion *> *)userData;
double lcGlob = CTX::instance()->lc;
- int useInterpolatedSize = CTX::instance()->mesh.lcExtendFromBoundary;
+ int useInterpolatedSize = Extend2dMeshIn3dVolumes();
HXT_INFO("Computing %smesh sizes...", useInterpolatedSize ? "interpolated " : "");
diff --git a/src/mesh/meshGRegionMMG.cpp b/src/mesh/meshGRegionMMG.cpp
index d5995068e2..852079f33c 100644
--- a/src/mesh/meshGRegionMMG.cpp
+++ b/src/mesh/meshGRegionMMG.cpp
@@ -172,7 +172,7 @@ static void gmsh2MMG(GRegion *gr, MMG5_pMesh mmg, MMG5_pSol sol,
auto itv = LCS.find(v);
if(itv != LCS.end()) {
mmg2gmsh[(*it)->getNum()] = *it;
- // if (CTX::instance()->mesh.lcExtendFromBoundary){
+ // if (Extend2dMeshIn3dVolumes()){
double LL = itv->second.first / itv->second.second;
SMetric3 l4(1. / (LL * LL));
SMetric3 MM = intersection_conserve_mostaniso(l4, m);
@@ -218,7 +218,7 @@ static void updateSizes(GRegion *gr, MMG5_pMesh mmg, MMG5_pSol sol,
std::vector<GFace *> f = gr->faces();
std::map<MVertex *, std::pair<double, int> > LCS;
- // if (CTX::instance()->mesh.lcExtendFromBoundary){
+ // if (Extend2dMeshIn3dVolumes()){
for(auto it = f.begin(); it != f.end(); ++it) {
for(unsigned int i = 0; i < (*it)->triangles.size(); i++) {
MTriangle *t = (*it)->triangles[i];
@@ -250,7 +250,7 @@ static void updateSizes(GRegion *gr, MMG5_pMesh mmg, MMG5_pSol sol,
auto it = mmg2gmsh.find(k);
- if(it != mmg2gmsh.end() && CTX::instance()->mesh.lcExtendFromBoundary) {
+ if(it != mmg2gmsh.end() && Extend2dMeshIn3dVolumes()) {
auto itv = LCS.find(it->second);
if(itv != LCS.end()) {
double LL = itv->second.first / itv->second.second;
diff --git a/src/parser/Gmsh.tab.cpp b/src/parser/Gmsh.tab.cpp
index 7fd34f5822..8778165852 100644
--- a/src/parser/Gmsh.tab.cpp
+++ b/src/parser/Gmsh.tab.cpp
@@ -12243,8 +12243,8 @@ yyreduce:
case 376:
#line 4771 "Gmsh.y"
{
- // lcExtendFromBoundary onstraints are stored in GEO internals in addition
- // to GModel, as they can be copied around during GEO operations
+ // mesh size from boundary onstraints are stored in GEO internals in
+ // addition to GModel, as they can be copied around during GEO operations
if(GModel::current()->getOCCInternals() &&
GModel::current()->getOCCInternals()->getChanged())
GModel::current()->getOCCInternals()->synchronize(GModel::current());
diff --git a/src/parser/Gmsh.y b/src/parser/Gmsh.y
index 900cd2644b..04e3a42972 100644
--- a/src/parser/Gmsh.y
+++ b/src/parser/Gmsh.y
@@ -4769,8 +4769,8 @@ Constraints :
}
| tMeshSizeFromBoundary tSurface '{' RecursiveListOfDouble '}' tAFFECT FExpr tEND
{
- // lcExtendFromBoundary onstraints are stored in GEO internals in addition
- // to GModel, as they can be copied around during GEO operations
+ // mesh size from boundary onstraints are stored in GEO internals in
+ // addition to GModel, as they can be copied around during GEO operations
if(GModel::current()->getOCCInternals() &&
GModel::current()->getOCCInternals()->getChanged())
GModel::current()->getOCCInternals()->synchronize(GModel::current());
diff --git a/tutorials/c++/t10.cpp b/tutorials/c++/t10.cpp
index d475662abe..4c0d1061c2 100644
--- a/tutorials/c++/t10.cpp
+++ b/tutorials/c++/t10.cpp
@@ -129,9 +129,9 @@ int main(int argc, char **argv)
// The value can then be further modified by the mesh size callback, if any,
// before being constrained in the interval [`Mesh.MeshSizeMin',
// `Mesh.MeshSizeMax'] and multiplied by `Mesh.MeshSizeFactor'. In addition,
- // boundary mesh sizes (on curves or surfaces) are interpolated inside the
- // enclosed entity (surface or volume, respectively) if the option
- // `Mesh.MeshSizeExtendFromBoundary' is set (which is the case by default).
+ // boundary mesh sizes are interpolated inside surfaces and/or volumes
+ // depending on the value of `Mesh.MeshSizeExtendFromBoundary' (which is set
+ // by default).
//
// When the element size is fully specified by a background mesh (as it is in
// this example), it is thus often desirable to set
diff --git a/tutorials/c/t2.c b/tutorials/c/t2.c
index 28cda7ef5a..1df37d3122 100644
--- a/tutorials/c/t2.c
+++ b/tutorials/c/t2.c
@@ -172,8 +172,8 @@ int main(int argc, char **argv)
int *ov2;
size_t ov2_n;
- int numElements[] = {};
- double heights[] = {};
+ int *numElements;
+ double *heights;
size_t numElements_n = 0;
size_t heights_n = 0;
int recombine = 0;
diff --git a/tutorials/python/t10.py b/tutorials/python/t10.py
index b11d8f6a87..c92b277565 100644
--- a/tutorials/python/t10.py
+++ b/tutorials/python/t10.py
@@ -121,9 +121,8 @@ gmsh.model.mesh.setSizeCallback(meshSizeCallback)
# The value can then be further modified by the mesh size callback, if any,
# before being constrained in the interval [`Mesh.MeshSizeMin',
# `Mesh.MeshSizeMax'] and multiplied by `Mesh.MeshSizeFactor'. In addition,
-# boundary mesh sizes (on curves or surfaces) are interpolated inside the
-# enclosed entity (surface or volume, respectively) if the option
-# `Mesh.MeshSizeExtendFromBoundary' is set (which is the case by default).
+# boundary mesh sizes are interpolated inside surfaces and/or volumes depending
+# on the value of `Mesh.MeshSizeExtendFromBoundary' (which is set by default).
#
# When the element size is fully specified by a background mesh (as it is in
# this example), it is thus often desirable to set
diff --git a/tutorials/t10.geo b/tutorials/t10.geo
index cf92f2dad3..8e39ea3ba7 100644
--- a/tutorials/t10.geo
+++ b/tutorials/t10.geo
@@ -96,10 +96,9 @@ Background Field = 7;
// 5) any per-entity mesh size constraint.
//
// This value is then constrained in the interval [`Mesh.MeshSizeMin',
-// `Mesh.MeshSizeMax'] and multiplied by `Mesh.MeshSizeFactor'. In addition,
-// boundary mesh sizes (on curves or surfaces) are interpolated inside the
-// enclosed entity (surface or volume, respectively) if the option
-// `Mesh.MeshSizeExtendFromBoundary' is set (which is the case by default).
+// `Mesh.MeshSizeMax'] and multiplied by `Mesh.MeshSizeFactor'. In addition,
+// boundary mesh sizes are interpolated inside surfaces and/or volumes depending
+// on the value of `Mesh.MeshSizeExtendFromBoundary' (which is set by default).
//
// When the element size is fully specified by a background mesh size field (as
// it is in this example), it is thus often desirable to set
--
GitLab