1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5
6 #include "BackgroundMeshTools.h"
7 #include "GFace.h"
8 #include "GVertex.h"
9 #include "GEdge.h"
10 #include "GEntity.h"
11 #include "Context.h"
12 #include "Field.h"
13 #include "GModel.h"
14
max_surf_curvature(const GEdge * ge,double u)15 static double max_surf_curvature(const GEdge *ge, double u)
16 {
17 double val = 0;
18 std::vector<GFace *> faces = ge->faces();
19 auto it = faces.begin();
20 while(it != faces.end()) {
21 SPoint2 par = ge->reparamOnFace((*it), u, 1);
22 double cc = (*it)->curvatureMax(par);
23 val = std::max(cc, val);
24 ++it;
25 }
26 return val;
27 }
28
max_edge_curvature(const GVertex * gv)29 static double max_edge_curvature(const GVertex *gv)
30 {
31 double val = 0;
32 std::vector<GEdge *> const &l_edges = gv->edges();
33 for(auto ite = l_edges.begin(); ite != l_edges.end(); ++ite) {
34 GEdge *_myGEdge = *ite;
35 Range<double> range = _myGEdge->parBounds(0);
36 double t = (gv == _myGEdge->getBeginVertex()) ? range.low() : range.high();
37 double EC = _myGEdge->curvature(t);
38 double FC = max_surf_curvature(_myGEdge, t);
39 val = std::max(std::max(val, EC), FC);
40 }
41 return val;
42 }
43
44 // the mesh vertex is classified on a model vertex. we compute the maximum of
45 // the curvature of model faces surrounding this point if it is classified on a
46 // model edge, we do the same for all model faces surrounding it if it is on a
47 // model face, we compute the curvature at this location
LC_MVertex_CURV(GEntity * ge,double U,double V)48 static double LC_MVertex_CURV(GEntity *ge, double U, double V)
49 {
50 double Crv = 0;
51 switch(ge->dim()) {
52 case 0:
53 Crv = max_edge_curvature((const GVertex *)ge);
54 // Crv = std::max(max_surf_curvature_vertex((const GVertex *)ge), Crv);
55 // Crv = max_surf_curvature((const GVertex *)ge);
56 break;
57 case 1: {
58 GEdge *ged = (GEdge *)ge;
59 Crv = ged->curvature(U);
60 Crv = std::max(Crv, max_surf_curvature(ged, U));
61 } break;
62 case 2: {
63 GFace *gf = (GFace *)ge;
64 Crv = gf->curvatureMax(SPoint2(U, V));
65 } break;
66 }
67 double ne = CTX::instance()->mesh.lcFromCurvature;
68 if(ne < 1) {
69 Msg::Warning("Invalid number of elements per 2*pi curvature %g", ne);
70 ne = 1;
71 }
72 double lc = Crv > 0 ? 2 * M_PI / Crv / ne : MAX_LC;
73 return lc;
74 }
75
max_edge_curvature_metric(const GEdge * ge,double u)76 SMetric3 max_edge_curvature_metric(const GEdge *ge, double u)
77 {
78 SVector3 t = ge->firstDer(u);
79 t.normalize();
80 double ne = CTX::instance()->mesh.lcFromCurvature;
81 if(ne < 1) {
82 Msg::Warning("Invalid number of elements per 2*pi curvature %g", ne);
83 ne = 1;
84 }
85 double l_t = (2 * M_PI) / (fabs(ge->curvature(u)) * ne);
86 double l_n = 1.e12;
87 return buildMetricTangentToCurve(t, l_t, l_n);
88 }
89
metric_based_on_surface_curvature(const GEdge * ge,double u,bool iso_surf)90 static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u,
91 bool iso_surf)
92 {
93 SMetric3 mesh_size(1.e-12);
94 std::vector<GFace *> faces = ge->faces();
95 auto it = faces.begin();
96 // we choose the metric eigenvectors to be the ones
97 // related to the edge ...
98 SMetric3 curvMetric = max_edge_curvature_metric(ge, u);
99 while(it != faces.end()) {
100 SPoint2 par = ge->reparamOnFace((*it), u, 1);
101 SMetric3 m =
102 metric_based_on_surface_curvature(*it, par.x(), par.y(), iso_surf);
103 curvMetric = intersection_conserveM1(curvMetric, m);
104 ++it;
105 }
106 return curvMetric;
107 }
108
metric_based_on_surface_curvature(const GVertex * gv,bool iso_surf)109 static SMetric3 metric_based_on_surface_curvature(const GVertex *gv,
110 bool iso_surf)
111 {
112 SMetric3 mesh_size(1.e-15);
113 return mesh_size;
114 }
115
LC_MVertex_CURV_ANISO(GEntity * ge,double U,double V)116 SMetric3 LC_MVertex_CURV_ANISO(GEntity *ge, double U, double V)
117 {
118 bool iso_surf = CTX::instance()->mesh.lcFromCurvatureIso;
119
120 switch(ge->dim()) {
121 case 0:
122 return metric_based_on_surface_curvature((const GVertex *)ge, iso_surf);
123 case 1:
124 return metric_based_on_surface_curvature((const GEdge *)ge, U, iso_surf);
125 case 2:
126 return metric_based_on_surface_curvature((const GFace *)ge, U, V, iso_surf);
127 }
128 Msg::Error("Curvature control impossible to compute for a volume");
129 return SMetric3();
130 }
131
132 // compute the mesh size at a given vertex due to prescribed sizes at
133 // mesh vertices
LC_MVertex_PNTS(GEntity * ge,double U,double V)134 static double LC_MVertex_PNTS(GEntity *ge, double U, double V)
135 {
136 switch(ge->dim()) {
137 case 0: {
138 GVertex *gv = (GVertex *)ge;
139 double lc = gv->prescribedMeshSizeAtVertex();
140 // FIXME we might want to remove this to make all lc treatment consistent
141 if(lc >= MAX_LC) return CTX::instance()->lc / 10.;
142 return lc;
143 }
144 case 1: {
145 GEdge *ged = (GEdge *)ge;
146 GVertex *v1 = ged->getBeginVertex();
147 GVertex *v2 = ged->getEndVertex();
148 if(v1 && v2) {
149 double lc1 = v1->prescribedMeshSizeAtVertex();
150 double lc2 = v2->prescribedMeshSizeAtVertex();
151 if(lc1 >= MAX_LC && lc2 >= MAX_LC) {
152 // FIXME we might want to remove this to make all lc treatment
153 // consistent
154 return CTX::instance()->lc / 10.;
155 }
156 else {
157 Range<double> range = ged->parBounds(0);
158 double a = (U - range.low()) / (range.high() - range.low());
159 return (1 - a) * lc1 + (a)*lc2;
160 }
161 }
162 else
163 return MAX_LC;
164 }
165 default: return MAX_LC;
166 }
167 }
168
buildMetricTangentToCurve(SVector3 & t,double l_t,double l_n)169 SMetric3 buildMetricTangentToCurve(SVector3 &t, double l_t, double l_n)
170 {
171 if(l_t == 0.0) return SMetric3(1.e-22);
172 SVector3 a;
173 if(fabs(t(0)) <= fabs(t(1)) && fabs(t(0)) <= fabs(t(2))) {
174 a = SVector3(1, 0, 0);
175 }
176 else if(fabs(t(1)) <= fabs(t(0)) && fabs(t(1)) <= fabs(t(2))) {
177 a = SVector3(0, 1, 0);
178 }
179 else {
180 a = SVector3(0, 0, 1);
181 }
182 SVector3 b = crossprod(t, a);
183 SVector3 c = crossprod(b, t);
184 b.normalize();
185 c.normalize();
186 t.normalize();
187 SMetric3 Metric(1. / (l_t * l_t), 1. / (l_n * l_n), 1. / (l_n * l_n), t, b,
188 c);
189 // printf("bmttc %g %g %g %g
190 // %g\n",l_t,l_n,Metric(0,0),Metric(0,1),Metric(1,1));
191 return Metric;
192 }
193
buildMetricTangentToSurface(SVector3 & t1,SVector3 & t2,double l_t1,double l_t2,double l_n)194 SMetric3 buildMetricTangentToSurface(SVector3 &t1, SVector3 &t2, double l_t1,
195 double l_t2, double l_n)
196 {
197 t1.normalize();
198 t2.normalize();
199 SVector3 n = crossprod(t1, t2);
200 n.normalize();
201
202 l_t1 = std::max(l_t1, CTX::instance()->mesh.lcMin);
203 l_t2 = std::max(l_t2, CTX::instance()->mesh.lcMin);
204 l_t1 = std::min(l_t1, CTX::instance()->mesh.lcMax);
205 l_t2 = std::min(l_t2, CTX::instance()->mesh.lcMax);
206 SMetric3 Metric(1. / (l_t1 * l_t1), 1. / (l_t2 * l_t2), 1. / (l_n * l_n), t1,
207 t2, n);
208 return Metric;
209 }
210
BGM_MeshSizeWithoutScaling(GEntity * ge,double U,double V,double X,double Y,double Z)211 double BGM_MeshSizeWithoutScaling(GEntity *ge, double U, double V, double X,
212 double Y, double Z)
213 {
214 // lc from points
215 double l1 = MAX_LC;
216 if(ge && CTX::instance()->mesh.lcFromPoints && ge->dim() < 2)
217 l1 = LC_MVertex_PNTS(ge, U, V);
218
219 // lc from curvature
220 double l2 = MAX_LC;
221 if(ge && CTX::instance()->mesh.lcFromCurvature > 0 && ge->dim() < 3)
222 l2 = LC_MVertex_CURV(ge, U, V);
223
224 // lc from fields
225 double l3 = MAX_LC;
226 if(ge) {
227 FieldManager *fields = ge->model()->getFields();
228 if(fields->getBackgroundField() > 0) {
229 Field *f = fields->get(fields->getBackgroundField());
230 if(f) l3 = (*f)(X, Y, Z, ge);
231 }
232 }
233
234 // global lc from entity
235 double l4 = ge ? ge->getMeshSize() : MAX_LC;
236
237 double l5 = (ge && ge->dim() == 1) ?
238 ((GEdge *)ge)->prescribedMeshSizeAtParam(U) :
239 MAX_LC;
240
241 // take the minimum
242 double lc = std::min(std::min(std::min(std::min(l1, l2), l3), l4), l5);
243
244 // lc from callback
245 if(GModel::current()->lcCallback) {
246 int dim = (ge ? ge->dim() : -1);
247 int tag = (ge ? ge->tag() : -1);
248 lc = GModel::current()->lcCallback(dim, tag, X, Y, Z, lc);
249 }
250
251 return lc;
252 }
253
254 // This is the only function that is used by the meshers
BGM_MeshSize(GEntity * ge,double U,double V,double X,double Y,double Z)255 double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y,
256 double Z)
257 {
258 if(!ge) Msg::Warning("No entity in background mesh size evaluation");
259
260 // default size to size of model
261 double lc = CTX::instance()->lc;
262
263 // min of all sizes
264 lc = std::min(lc, BGM_MeshSizeWithoutScaling(ge, U, V, X, Y, Z));
265
266 // constrain by lcMin and lcMax
267 lc = std::max(lc, CTX::instance()->mesh.lcMin);
268 lc = std::min(lc, CTX::instance()->mesh.lcMax);
269
270 if(lc <= 0.) {
271 Msg::Error("Wrong mesh element size lc = %g (lcmin = %g, lcmax = %g)", lc,
272 CTX::instance()->mesh.lcMin, CTX::instance()->mesh.lcMax);
273 lc = CTX::instance()->lc;
274 }
275
276 // size factor from entity
277 if(ge && ge->getMeshSizeFactor() != 1.0) lc *= ge->getMeshSizeFactor();
278
279 // global size factor
280 return lc * CTX::instance()->mesh.lcFactor;
281 }
282
283 // anisotropic version of the background field
BGM_MeshMetric(GEntity * ge,double U,double V,double X,double Y,double Z)284 SMetric3 BGM_MeshMetric(GEntity *ge, double U, double V, double X, double Y,
285 double Z)
286 {
287 // default size to size of model
288 double lc = CTX::instance()->lc;
289
290 // lc from points
291 if(CTX::instance()->mesh.lcFromPoints && ge->dim() < 2)
292 lc = std::min(lc, LC_MVertex_PNTS(ge, U, V));
293
294 // global lc from entity
295 lc = std::min(lc, ge->getMeshSize());
296
297 // constrain by global min/max
298 lc = std::max(lc, CTX::instance()->mesh.lcMin);
299 lc = std::min(lc, CTX::instance()->mesh.lcMax);
300
301 if(lc <= 0.) {
302 Msg::Error("Wrong mesh element size lc = %g (lcmin = %g, lcmax = %g)", lc,
303 CTX::instance()->mesh.lcMin, CTX::instance()->mesh.lcMax);
304 lc = CTX::instance()->lc;
305 }
306
307 // isotropic base metric
308 SMetric3 m0(1. / (lc * lc));
309
310 // intersect with metrics from fields if applicable
311 FieldManager *fields = ge->model()->getFields();
312 SMetric3 m1 = m0;
313 if(fields->getBackgroundField() > 0) {
314 Field *f = fields->get(fields->getBackgroundField());
315 if(f) {
316 SMetric3 l4;
317 if(!f->isotropic()) { (*f)(X, Y, Z, l4, ge); }
318 else {
319 const double L = (*f)(X, Y, Z, ge);
320 l4 = SMetric3(1 / (L * L));
321 }
322 m1 = intersection(l4, m0);
323 }
324 }
325
326 // intersect with metrics from curvature if applicable
327 SMetric3 m = (CTX::instance()->mesh.lcFromCurvature > 0 && ge->dim() < 3) ?
328 intersection(m1, LC_MVertex_CURV_ANISO(ge, U, V)) :
329 m1;
330
331 // apply global size factor
332 if(CTX::instance()->mesh.lcFactor != 0 &&
333 CTX::instance()->mesh.lcFactor != 1.)
334 m *= 1. / (CTX::instance()->mesh.lcFactor * CTX::instance()->mesh.lcFactor);
335
336 return m;
337 }
338
Extend1dMeshIn2dSurfaces(GFace * gf)339 bool Extend1dMeshIn2dSurfaces(GFace *gf)
340 {
341 return gf->getMeshSizeFromBoundary();
342 }
343
Extend2dMeshIn3dVolumes()344 bool Extend2dMeshIn3dVolumes()
345 {
346 return CTX::instance()->mesh.lcExtendFromBoundary;
347 }
348
max_edge_curvature_metric(const GVertex * gv)349 SMetric3 max_edge_curvature_metric(const GVertex *gv)
350 {
351 SMetric3 val(1.e-12);
352 std::vector<GEdge *> const &l_edges = gv->edges();
353 for(auto ite = l_edges.begin(); ite != l_edges.end(); ++ite) {
354 GEdge *_myGEdge = *ite;
355 Range<double> range = _myGEdge->parBounds(0);
356 SMetric3 cc;
357 double ne = CTX::instance()->mesh.lcFromCurvature;
358 if(ne < 1) {
359 Msg::Warning("Invalid number of elements per 2*pi curvature %g", ne);
360 ne = 1;
361 }
362 if(gv == _myGEdge->getBeginVertex()) {
363 SVector3 t = _myGEdge->firstDer(range.low());
364 t.normalize();
365 double l_t = 2 * M_PI / (fabs(_myGEdge->curvature(range.low())) * ne);
366 double l_n = 1.e12;
367 cc = buildMetricTangentToCurve(t, l_t, l_n);
368 }
369 else {
370 SVector3 t = _myGEdge->firstDer(range.high());
371 t.normalize();
372 double l_t = 2 * M_PI / (fabs(_myGEdge->curvature(range.high())) * ne);
373 double l_n = 1.e12;
374 cc = buildMetricTangentToCurve(t, l_t, l_n);
375 }
376 val = intersection(val, cc);
377 }
378 return val;
379 }
380
metric_based_on_surface_curvature(const GFace * gf,double u,double v,bool surface_isotropic,double d_normal,double d_tangent_max)381 SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v,
382 bool surface_isotropic,
383 double d_normal,
384 double d_tangent_max)
385 {
386 if(gf->geomType() == GEntity::Plane) return SMetric3(1.e-12);
387 double cmax, cmin;
388 SVector3 dirMax, dirMin;
389 cmax = gf->curvatures(SPoint2(u, v), dirMax, dirMin, cmax, cmin);
390 if(cmin == 0) cmin = 1.e-12;
391 if(cmax == 0) cmax = 1.e-12;
392 double ne = CTX::instance()->mesh.lcFromCurvature;
393 if(ne < 1) {
394 Msg::Warning("Invalid number of elements per 2*pi curvature %g", ne);
395 ne = 1;
396 }
397 double lambda1 = 2 * M_PI / (fabs(cmin) * ne);
398 double lambda2 = 2 * M_PI / (fabs(cmax) * ne);
399 SVector3 Z = crossprod(dirMax, dirMin);
400 if(surface_isotropic) lambda2 = lambda1 = std::min(lambda2, lambda1);
401 dirMin.normalize();
402 dirMax.normalize();
403 Z.normalize();
404 lambda1 = std::max(lambda1, CTX::instance()->mesh.lcMin);
405 lambda2 = std::max(lambda2, CTX::instance()->mesh.lcMin);
406 lambda1 = std::min(lambda1, CTX::instance()->mesh.lcMax);
407 lambda2 = std::min(lambda2, CTX::instance()->mesh.lcMax);
408 double lambda3 = std::min(d_normal, CTX::instance()->mesh.lcMax);
409 lambda3 = std::max(lambda3, CTX::instance()->mesh.lcMin);
410 lambda1 = std::min(lambda1, d_tangent_max);
411 lambda2 = std::min(lambda2, d_tangent_max);
412
413 SMetric3 curvMetric(1. / (lambda1 * lambda1), 1. / (lambda2 * lambda2),
414 1. / (lambda3 * lambda3), dirMin, dirMax, Z);
415 return curvMetric;
416 }
417