1 //ff-c++-LIBRARY-dep: mmg scotch
2 //ff-c++-cpp-dep:
3
4 #include "ff++.hpp"
5 #include "memory.h"
6 #include "mmg/libmmg.h"
7 #include "GenericMesh.hpp"
8
9 using namespace Fem2D;
10
11 template<class ffmesh> int ffmesh_to_MMG5_pMesh(const ffmesh &, MMG5_pMesh&);
12
13 template<>
ffmesh_to_MMG5_pMesh(const Mesh3 & Th,MMG5_pMesh & mesh)14 int ffmesh_to_MMG5_pMesh<Mesh3>(const Mesh3 &Th, MMG5_pMesh& mesh) {
15 int nVertices = Th.nv;
16 int nTetrahedra = Th.nt;
17 int nPrisms = 0;
18 int nTriangles = Th.nbe;
19 int nQuadrilaterals = 0;
20 int nEdges = 0;
21
22 if ( MMG3D_Set_meshSize(mesh,nVertices,nTetrahedra,nPrisms,nTriangles,
23 nQuadrilaterals,nEdges) != 1 ) {
24 exit(EXIT_FAILURE);
25 }
26
27 for (int k = 0; k < Th.nv; k++) {
28 if ( MMG3D_Set_vertex(mesh,Th.vertices[k].x,Th.vertices[k].y,
29 Th.vertices[k].z, Th.vertices[k].lab, k+1) != 1 ) {
30 exit(EXIT_FAILURE);
31 }
32 }
33
34 for (int k = 0; k < Th.nt; k++) {
35 const Tet &K(Th.elements[k]);
36 if ( MMG3D_Set_tetrahedron(mesh,Th.operator()(K[0])+1,Th.operator()(K[1])+1,
37 Th.operator()(K[2])+1,Th.operator()(K[3])+1,K.lab,k+1) != 1 ) {
38 exit(EXIT_FAILURE);
39 }
40 }
41
42 for (int k = 0; k < Th.nbe; k++) {
43 const Triangle3 &K(Th.be(k));
44 if ( MMG3D_Set_triangle(mesh,Th.operator()(K[0])+1,Th.operator()(K[1])+1,Th.operator()(K[2])+1,
45 K.lab,k+1) != 1 ) {
46 exit(EXIT_FAILURE);
47 }
48 }
49
50 return 0;
51 }
52
53 template<>
ffmesh_to_MMG5_pMesh(const MeshS & Th,MMG5_pMesh & mesh)54 int ffmesh_to_MMG5_pMesh<MeshS>(const MeshS &Th, MMG5_pMesh& mesh) {
55 int nVertices = Th.nv;
56 int nTetrahedra = 0;
57 int nPrisms = 0;
58 int nTriangles = Th.nt;
59 int nQuadrilaterals = 0;
60 int nEdges = 0;
61
62 if ( MMGS_Set_meshSize(mesh,nVertices,nTriangles,nEdges) != 1 ) {
63 exit(EXIT_FAILURE);
64 }
65
66 for (int k = 0; k < Th.nv; k++) {
67 if ( MMGS_Set_vertex(mesh,Th.vertices[k].x,Th.vertices[k].y,
68 Th.vertices[k].z, Th.vertices[k].lab, k+1) != 1 ) {
69 exit(EXIT_FAILURE);
70 }
71 }
72
73 for (int k = 0; k < Th.nt; k++) {
74 const TriangleS &K(Th.elements[k]);
75 if ( MMGS_Set_triangle(mesh,Th.operator()(K[0])+1,Th.operator()(K[1])+1,Th.operator()(K[2])+1,
76 K.lab,k+1) != 1 ) {
77 exit(EXIT_FAILURE);
78 }
79 }
80
81 return 0;
82 }
83
84 template<class ffmesh> int MMG5_pMesh_to_ffmesh(const MMG5_pMesh&, ffmesh *&);
85
86 template<>
MMG5_pMesh_to_ffmesh(const MMG5_pMesh & mesh,Mesh3 * & T_TH3)87 int MMG5_pMesh_to_ffmesh<Mesh3>(const MMG5_pMesh& mesh, Mesh3 *&T_TH3) {
88 int ier;
89
90 int nVertices = 0;
91 int nTetrahedra = 0;
92 int nTriangles = 0;
93 int nEdges = 0;
94
95 if ( MMG3D_Get_meshSize(mesh,&nVertices,&nTetrahedra,NULL,&nTriangles,NULL,
96 &nEdges) !=1 ) {
97 ier = MMG5_STRONGFAILURE;
98 }
99
100 Vertex3 *v = new Vertex3[nVertices];
101 Tet *t = new Tet[nTetrahedra];
102 Tet *tt = t;
103 Triangle3 *b = new Triangle3[nTriangles];
104 Triangle3 *bb = b;
105 int k;
106
107 int corner, required;
108
109 for (k = 0; k < nVertices; k++) {
110 if ( MMG3D_Get_vertex(mesh,&(v[k].x),&(v[k].y),&(v[k].z),
111 &(v[k].lab),&(corner),&(required)) != 1 ) {
112 cout << "Unable to get mesh vertex " << k << endl;
113 ier = MMG5_STRONGFAILURE;
114 }
115 }
116
117 for ( k=0; k<nTetrahedra; k++ ) {
118 int iv[4], lab;
119 if ( MMG3D_Get_tetrahedron(mesh,
120 &(iv[0]),&(iv[1]),
121 &(iv[2]),&(iv[3]),
122 &(lab),&(required)) != 1 ) {
123 cout << "Unable to get mesh tetra " << k << endl;
124 ier = MMG5_STRONGFAILURE;
125 }
126 for (int i=0; i<4; i++)
127 iv[i]--;
128 tt++->set(v, iv, lab);
129 }
130
131 for ( k=0; k<nTriangles; k++ ) {
132 int iv[3], lab;
133 if ( MMG3D_Get_triangle(mesh,
134 &(iv[0]),&(iv[1]),&(iv[2]),
135 &(lab),&(required)) != 1 ) {
136 cout << "Unable to get mesh triangle " << k << endl;
137 ier = MMG5_STRONGFAILURE;
138 }
139 for (int i=0; i<3; i++)
140 iv[i]--;
141 bb++->set(v, iv, lab);
142 }
143
144 T_TH3 = new Mesh3(nVertices, nTetrahedra, nTriangles, v, t, b, 1);
145
146 if (verbosity > 1) {
147 cout << "transformation maillage --> mesh3 " << endl;
148 cout << "vertices =" << nVertices << endl;
149 cout << "tetrahedrons =" << nTetrahedra << endl;
150 cout << "triangles =" << nTriangles << endl;
151 cout << "T_TH3" << T_TH3->nv << " " << T_TH3->nt << " " << T_TH3->nbe << endl;
152 }
153
154 return 0;
155 }
156
157 template<>
MMG5_pMesh_to_ffmesh(const MMG5_pMesh & mesh,MeshS * & T_TH3)158 int MMG5_pMesh_to_ffmesh<MeshS>(const MMG5_pMesh& mesh, MeshS *&T_TH3) {
159 int ier;
160
161 int nVertices = 0;
162 int nTriangles = 0;
163 int nEdges = 0;
164
165 if ( MMGS_Get_meshSize(mesh,&nVertices,&nTriangles,&nEdges) !=1 ) {
166 ier = MMG5_STRONGFAILURE;
167 }
168
169 Vertex3 *v = new Vertex3[nVertices];
170 TriangleS *t = new TriangleS[nTriangles];
171 TriangleS *tt = t;
172 BoundaryEdgeS *b = new BoundaryEdgeS[nEdges];
173 BoundaryEdgeS *bb = b;
174 int k;
175
176 int corner, required, ridge;
177
178 for (k = 0; k < nVertices; k++) {
179 if ( MMGS_Get_vertex(mesh,&(v[k].x),&(v[k].y),&(v[k].z),
180 &(v[k].lab),&(corner),&(required)) != 1 ) {
181 cout << "Unable to get mesh vertex " << k << endl;
182 ier = MMG5_STRONGFAILURE;
183 }
184 }
185
186 for ( k=0; k<nTriangles; k++ ) {
187 int iv[3], lab;
188 if ( MMGS_Get_triangle(mesh,
189 &(iv[0]),&(iv[1]),&(iv[2]),
190 &(lab),&(required)) != 1 ) {
191 cout << "Unable to get mesh triangle " << k << endl;
192 ier = MMG5_STRONGFAILURE;
193 }
194 for (int i=0; i<3; i++)
195 iv[i]--;
196 tt++->set(v, iv, lab);
197 }
198
199 for ( k=0; k<nEdges; k++ ) {
200 int iv[2], lab;
201 if ( MMGS_Get_edge(mesh,
202 &(iv[0]),&(iv[1]),
203 &(lab), &(ridge), &(required)) != 1 ) {
204 cout << "Unable to get mesh edge " << k << endl;
205 ier = MMG5_STRONGFAILURE;
206 }
207 for (int i=0; i<2; i++)
208 iv[i]--;
209 bb++->set(v, iv, lab);
210 }
211
212 T_TH3 = new MeshS(nVertices, nTriangles, nEdges, v, t, b, 1);
213
214 if (verbosity > 1) {
215 cout << "transformation maillage --> meshS " << endl;
216 cout << "vertices =" << nVertices << endl;
217 cout << "triangles =" << nTriangles << endl;
218 cout << "edges =" << nEdges << endl;
219 cout << "T_TH3" << T_TH3->nv << " " << T_TH3->nt << " " << T_TH3->nbe << endl;
220 }
221
222 return 0;
223 }
224
225 template<class ffmesh>
226 class mmg_Op : public E_F0mps {
227 public:
228 Expression eTh, xx, yy, zz;
229 static const int n_name_param = std::is_same<ffmesh,Mesh3>::value ? 27 : 19;
230 static basicAC_F0::name_and_type name_param[];
231 Expression nargs[n_name_param];
232
arg(int i,Stack stack,KN_<long> a) const233 KN_< long > arg(int i, Stack stack, KN_< long > a) const {
234 return nargs[i] ? GetAny< KN_< long > >((*nargs[i])(stack)) : a;
235 }
236
arg(int i,Stack stack,KN_<double> a) const237 KN_< double > arg(int i, Stack stack, KN_< double > a) const {
238 return nargs[i] ? GetAny< KN_< double > >((*nargs[i])(stack)) : a;
239 }
240
arg(int i,Stack stack,double a) const241 double arg(int i, Stack stack, double a) const {
242 return nargs[i] ? GetAny< double >((*nargs[i])(stack)) : a;
243 }
244
arg(int i,Stack stack,long a) const245 long arg(int i, Stack stack, long a) const {
246 return nargs[i] ? GetAny< long >((*nargs[i])(stack)) : a;
247 }
248
249 public:
mmg_Op(const basicAC_F0 & args,Expression tth)250 mmg_Op(const basicAC_F0 &args, Expression tth) : eTh(tth) {
251 if (verbosity > 1) {
252 cout << "mmg" << endl;
253 }
254
255 args.SetNameParam(n_name_param, name_param, nargs);
256
257 const E_Array *a1 = 0;
258 if (nargs[2]) {
259 a1 = dynamic_cast< const E_Array * >(nargs[2]);
260 }
261
262 if (a1) {
263 if (a1->size( ) != 3) {
264 CompileError("mmg3d(Th,displacement=[X,Y,Z],) ");
265 }
266
267 xx = to< double >((*a1)[0]);
268 yy = to< double >((*a1)[1]);
269 zz = to< double >((*a1)[2]);
270 }
271 }
272
273 AnyType operator( )(Stack stack) const;
274 };
275
276 template<>
277 basicAC_F0::name_and_type mmg_Op<Mesh3>::name_param[] = {
278 {"metric" , &typeid(KN< double > *)},
279 {"verbose" , &typeid(long)},/*!< [-1..10], Tune level of verbosity */
280 {"mem" , &typeid(long)},/*!< [n/-1], Set memory size to n Mbytes or keep the default value */
281 {"debug" , &typeid(bool)},/*!< [1/0], Turn on/off debug mode */
282 {"angle" , &typeid(bool)},/*!< [1/0], Turn on/off angle detection */
283 {"iso" , &typeid(bool)},/*!< [1/0], Level-set meshing */
284 {"nofem" , &typeid(bool)},/*!< [1/0], Generate a non finite element mesh */
285 {"opnbdy" , &typeid(bool)},/*!< [1/0], Preserve triangles at interface of 2 domains with same reference */
286 {"lag" , &typeid(long)},/*!< [-1/0/1/2], Lagrangian option */
287 {"optim" , &typeid(bool)},/*!< [1/0], Optimize mesh keeping its initial edge sizes */
288 {"optimLES" , &typeid(bool)},/*!< [1/0], Strong mesh optimization for Les computations */
289 {"noinsert" , &typeid(bool)},/*!< [1/0], Avoid/allow point insertion */
290 {"noswap" , &typeid(bool)},/*!< [1/0], Avoid/allow edge or face flipping */
291 {"nomove" , &typeid(bool)},/*!< [1/0], Avoid/allow point relocation */
292 {"nosurf" , &typeid(bool)},/*!< [1/0], Avoid/allow surface modifications */
293 //{"nreg" , &typeid(bool)},/*!< [0/1], Enable normal regularization */
294 {"renum" , &typeid(bool)},/*!< [1/0], Turn on/off point relocation with Scotch */
295 {"anisosize" , &typeid(bool)},/*!< [1/0], Turn on/off anisotropic metric creation when no metric is provided */
296 {"octree" , &typeid(long)},/*!< [n], Specify the max number of points per PROctree cell (DELAUNAY) */
297 {"angleDetection" , &typeid(double)},/*!< [val], Value for angle detection */
298 {"hmin" , &typeid(double)},/*!< [val], Minimal mesh size */
299 {"hmax" , &typeid(double)},/*!< [val], Maximal mesh size */
300 {"hsiz" , &typeid(double)},/*!< [val], Constant mesh size */
301 {"hausd" , &typeid(double)},/*!< [val], Control global Hausdorff distance (on all the boundary surfaces of the mesh) */
302 {"hgrad" , &typeid(double)},/*!< [val], Control gradation */
303 {"ls" , &typeid(double)},/*!< [val], Value of level-set */
304 //{"rmc" , &typeid(double)},/*!< [-1/val], Remove small connex componants in level-set mode */
305 {"requiredTriangle" , &typeid(KN<long>*)},/*!< [val], Value of level-set */
306 {"localParameter" , &typeid(KNM<double>*)}/*!< [val], Value of level-set */
307 };
308
309 template<>
310 basicAC_F0::name_and_type mmg_Op<MeshS>::name_param[] = {
311 {"metric" , &typeid(KN< double > *)},
312 {"verbose" , &typeid(long)},/*!< [-1..10], Tune level of verbosity */
313 {"mem" , &typeid(long)},/*!< [n/-1], Set memory size to n Mbytes or keep the default value */
314 {"debug" , &typeid(bool)},/*!< [1/0], Turn on/off debug mode */
315 {"angle" , &typeid(bool)},/*!< [1/0], Turn on/off angle detection */
316 {"iso" , &typeid(bool)},/*!< [1/0], Level-set meshing */
317 {"keepRef" , &typeid(bool)},/*!< [1/0], Preserve the initial domain references in level-set mode */
318 //{"optim" , &typeid(bool)},/*!< [1/0], Optimize mesh keeping its initial edge sizes */
319 {"noinsert" , &typeid(bool)},/*!< [1/0], Avoid/allow point insertion */
320 {"noswap" , &typeid(bool)},/*!< [1/0], Avoid/allow edge or face flipping */
321 {"nomove" , &typeid(bool)},/*!< [1/0], Avoid/allow point relocation */
322 {"nreg" , &typeid(bool)},/*!< [0/1], Disabled/enabled normal regularization */
323 {"renum" , &typeid(bool)},/*!< [1/0], Turn on/off point relocation with Scotch */
324 {"angleDetection" , &typeid(double)},/*!< [val], Value for angle detection */
325 {"hmin" , &typeid(double)},/*!< [val], Minimal mesh size */
326 {"hmax" , &typeid(double)},/*!< [val], Maximal mesh size */
327 {"hsiz" , &typeid(double)},/*!< [val], Constant mesh size */
328 {"hausd" , &typeid(double)},/*!< [val], Control global Hausdorff distance (on all the boundary surfaces of the mesh) */
329 {"hgrad" , &typeid(double)},/*!< [val], Control gradation */
330 {"ls" , &typeid(double)} /*!< [val], Value of level-set */
331 };
332
333 template<class ffmesh>
334 class mmg_ff : public OneOperator {
335 public:
mmg_ff()336 mmg_ff( ) : OneOperator(atype< const ffmesh* >( ), atype< const ffmesh* >( )) {}
337
code(const basicAC_F0 & args) const338 E_F0 *code(const basicAC_F0 &args) const { return new mmg_Op<ffmesh>(args, t[0]->CastTo(args[0])); }
339 };
340
341 template<>
operator ( )(Stack stack) const342 AnyType mmg_Op<Mesh3>::operator( )(Stack stack) const {
343 // initialisation
344 MeshPoint *mp(MeshPointStack(stack)), mps = *mp;
345 Mesh3 *pTh = GetAny< Mesh3 * >((*eTh)(stack));
346
347 ffassert(pTh);
348 Mesh3 &Th = *pTh;
349 int nv = Th.nv;
350 int nt = Th.nt;
351 int nbe = Th.nbe;
352
353 KN< double > *pmetric = 0;
354 KN< long > *prequiredTriangle = 0;
355 KNM< double > *plocalParameter = 0;
356
357 if (nargs[0]) {
358 pmetric = GetAny< KN< double > * >((*nargs[0])(stack));
359 }
360 if (nargs[25]) {
361 prequiredTriangle = GetAny< KN< long > * >((*nargs[25])(stack));
362 }
363 if (nargs[26]) {
364 plocalParameter = GetAny< KNM< double > * >((*nargs[26])(stack));
365 }
366
367 MMG5_pMesh mesh;
368 MMG5_pSol sol;
369 MMG5_pSol met;
370
371 mesh = nullptr;
372 sol = nullptr;
373 met = nullptr;
374
375 MMG3D_Init_mesh(MMG5_ARG_start,
376 MMG5_ARG_ppMesh,&mesh,MMG5_ARG_ppMet,&sol,
377 MMG5_ARG_end);
378
379 ffmesh_to_MMG5_pMesh(Th, mesh);
380
381 if (pmetric && pmetric->N( ) > 0) {
382 const KN< double > &metric = *pmetric;
383 if (metric.N( ) == Th.nv) {
384 if ( MMG3D_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Scalar) != 1 ) {
385 printf("Unable to allocate the metric array.\n");
386 exit(EXIT_FAILURE);
387 }
388 if ( MMG3D_Set_scalarSols(sol,metric) != 1 ) {
389 printf("Unable to set metric.\n");
390 exit(EXIT_FAILURE);
391 }
392 }
393 else {
394 if ( MMG3D_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Tensor) != 1 ) {
395 printf("Unable to allocate the metric array.\n");
396 exit(EXIT_FAILURE);
397 }
398 static const int perm[6] = {0, 1, 3, 2, 4, 5};
399 for (int k=0; k<Th.nv; k++) {
400 if ( MMG3D_Set_tensorSol(sol, metric[6*k+perm[0]], metric[6*k+perm[1]], metric[6*k+perm[2]],
401 metric[6*k+perm[3]], metric[6*k+perm[4]], metric[6*k+perm[5]], k+1) != 1 ) {
402 printf("Unable to set metric.\n");
403 exit(EXIT_FAILURE);
404 }
405 }
406 }
407 }
408 if (prequiredTriangle && prequiredTriangle->N( ) > 0) {
409 const KN< long > &requiredTriangle = *prequiredTriangle;
410 std::sort(requiredTriangle + 0, requiredTriangle + requiredTriangle.N());
411 int nt;
412 if ( MMG3D_Get_meshSize(mesh,NULL,NULL,NULL,&nt,NULL,NULL) !=1 ) {
413 exit(EXIT_FAILURE);
414 }
415 for (int k=1; k<=nt; k++) {
416 int ref, dummy;
417 if ( MMG3D_Get_triangle(mesh,&dummy,&dummy,&dummy,
418 &ref,NULL) != 1 ) {
419 exit(EXIT_FAILURE);
420 }
421 if (std::binary_search(requiredTriangle + 0, requiredTriangle + requiredTriangle.N(), ref)) {
422 if ( MMG3D_Set_requiredTriangle(mesh,k) != 1 ) {
423 exit(EXIT_FAILURE);
424 }
425 }
426 }
427 }
428 if (plocalParameter && plocalParameter->M( ) > 0) {
429 const KNM< double > &localParameter = *plocalParameter;
430 ffassert(localParameter.N() == 4);
431 if ( MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_numberOfLocalParam,localParameter.M()) != 1 ) {
432 exit(EXIT_FAILURE);
433 }
434 for(int j = 0; j < localParameter.M(); ++j) {
435 if ( MMG3D_Set_localParameter(mesh,sol,MMG5_Triangle,localParameter(0,j),localParameter(1,j),localParameter(2,j),localParameter(3,j)) != 1 ) {
436 exit(EXIT_FAILURE);
437 }
438 }
439 }
440
441 int i=1;
442 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_verbose, arg(i,stack,0L)); i++; /*!< [-1..10], Tune level of verbosity */
443 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_mem, arg(i,stack,0L)); i++; /*!< [n/-1], Set memory size to n Mbytes or keep the default value */
444 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_debug, arg(i,stack,0L)); i++; /*!< [1/0], Turn on/off debug mode */
445 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_angle, arg(i,stack,0L)); i++; /*!< [1/0], Turn on/off angle detection */
446 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_iso, arg(i,stack,0L)); i++; /*!< [1/0], Level-set meshing */
447 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_nofem, arg(i,stack,0L)); i++; /*!< [1/0], Generate a non finite element mesh */
448 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_opnbdy, arg(i,stack,0L)); i++; /*!< [1/0], Preserve triangles at interface of 2 domains with same reference */
449 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_lag, arg(i,stack,0L)); i++; /*!< [-1/0/1/2], Lagrangian option */
450 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_optim, arg(i,stack,0L)); i++; /*!< [1/0], Optimize mesh keeping its initial edge sizes */
451 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_optimLES, arg(i,stack,0L)); i++; /*!< [1/0], Strong mesh optimization for Les computations */
452 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_noinsert, arg(i,stack,0L)); i++; /*!< [1/0], Avoid/allow point insertion */
453 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_noswap, arg(i,stack,0L)); i++; /*!< [1/0], Avoid/allow edge or face flipping */
454 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_nomove, arg(i,stack,0L)); i++; /*!< [1/0], Avoid/allow point relocation */
455 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_nosurf, arg(i,stack,0L)); i++; /*!< [1/0], Avoid/allow surface modifications */
456 //if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_nreg, arg(i,stack,0L)); i++; /*!< [0/1], Enable normal regularization */
457 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_renum, arg(i,stack,0L)); i++; /*!< [1/0], Turn on/off point relocation with Scotch */
458 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_anisosize, arg(i,stack,0L)); i++; /*!< [1/0], Turn on/off anisotropic metric creation when no metric is provided */
459 if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_octree, arg(i,stack,0L)); i++; /*!< [n], Specify the max number of points per PROctree cell (DELAUNAY) */
460 if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_angleDetection,arg(i,stack,0.)); i++; /*!< [val], Value for angle detection */
461 if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hmin, arg(i,stack,0.)); i++; /*!< [val], Minimal mesh size */
462 if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hmax, arg(i,stack,0.)); i++; /*!< [val], Maximal mesh size */
463 if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hsiz, arg(i,stack,0.)); i++; /*!< [val], Constant mesh size */
464 if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hausd, arg(i,stack,0.)); i++; /*!< [val], Control global Hausdorff distance (on all the boundary surfaces of the mesh) */
465 if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hgrad, arg(i,stack,0.)); i++; /*!< [val], Control gradation */
466 if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_ls, arg(i,stack,0.)); i++; /*!< [val], Value of level-set */
467 //if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_rmc, arg(i,stack,0.)); i++; /*!< [-1/val], Remove small connex componants in level-set mode */
468
469 bool bls = MMG3D_Get_iparameter(mesh,MMG3D_IPARAM_iso);
470
471 int ier;
472 if (!bls)
473 ier = MMG3D_mmg3dlib(mesh,sol);
474 else
475 ier = MMG3D_mmg3dls(mesh,sol,met);
476
477 /*
478 if ( MMG5_saveMesh_centralized(mesh,"test.mesh") != 1 ) {
479 fprintf(stdout,"UNABLE TO SAVE MESH\n");
480 ier = MMG5_STRONGFAILURE;
481 }
482 MMG5_saveMet_centralized(mesh,"test.meshb");
483 */
484
485 Mesh3 *Th_T = nullptr;
486
487 MMG5_pMesh_to_ffmesh(mesh,Th_T);
488
489 MMG3D_Free_all(MMG5_ARG_start,
490 MMG5_ARG_ppMesh,&mesh,MMG5_ARG_ppMet,&sol,
491 MMG5_ARG_end);
492
493 Th_T->BuildGTree();
494
495 Add2StackOfPtr2FreeRC(stack, Th_T);
496 return Th_T;
497 }
498
499 template<>
operator ( )(Stack stack) const500 AnyType mmg_Op<MeshS>::operator( )(Stack stack) const {
501 // initialisation
502 MeshPoint *mp(MeshPointStack(stack)), mps = *mp;
503 MeshS *pTh = GetAny< MeshS * >((*eTh)(stack));
504
505 ffassert(pTh);
506 MeshS &Th = *pTh;
507 int nv = Th.nv;
508 int nt = Th.nt;
509 int nbe = Th.nbe;
510
511 KN< double > *pmetric = 0;
512
513 if (nargs[0]) {
514 pmetric = GetAny< KN< double > * >((*nargs[0])(stack));
515 }
516
517 MMG5_pMesh mesh;
518 MMG5_pSol sol;
519
520 mesh = nullptr;
521 sol = nullptr;
522
523 MMGS_Init_mesh(MMG5_ARG_start,
524 MMG5_ARG_ppMesh,&mesh,MMG5_ARG_ppMet,&sol,
525 MMG5_ARG_end);
526
527 ffmesh_to_MMG5_pMesh(Th, mesh);
528
529 if (pmetric && pmetric->N( ) > 0) {
530 const KN< double > &metric = *pmetric;
531 if (metric.N( ) == Th.nv) {
532 if ( MMGS_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Scalar) != 1 ) {
533 printf("Unable to allocate the metric array.\n");
534 exit(EXIT_FAILURE);
535 }
536 if ( MMGS_Set_scalarSols(sol,metric) != 1 ) {
537 printf("Unable to set metric.\n");
538 exit(EXIT_FAILURE);
539 }
540 }
541 else {
542 if ( MMGS_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Tensor) != 1 ) {
543 printf("Unable to allocate the metric array.\n");
544 exit(EXIT_FAILURE);
545 }
546 static const int perm[6] = {0, 1, 3, 2, 4, 5};
547 for (int k=0; k<Th.nv; k++) {
548 if ( MMGS_Set_tensorSol(sol, metric[6*k+perm[0]], metric[6*k+perm[1]], metric[6*k+perm[2]],
549 metric[6*k+perm[3]], metric[6*k+perm[4]], metric[6*k+perm[5]], k+1) != 1 ) {
550 printf("Unable to set metric.\n");
551 exit(EXIT_FAILURE);
552 }
553 }
554 }
555 }
556
557 int i=1;
558 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_verbose, arg(i,stack,0L)); i++; /*!< [-1..10], Tune level of verbosity */
559 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_mem, arg(i,stack,0L)); i++; /*!< [n/-1], Set memory size to n Mbytes or keep the default value */
560 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_debug, arg(i,stack,0L)); i++; /*!< [1/0], Turn on/off debug mode */
561 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_angle, arg(i,stack,0L)); i++; /*!< [1/0], Turn on/off angle detection */
562 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_iso, arg(i,stack,0L)); i++; /*!< [1/0], Level-set meshing */
563 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_keepRef, arg(i,stack,0L)); i++; /*!< [1/0], Preserve the initial domain references in level-set mode */
564 //if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_optim, arg(i,stack,0L)); i++; /*!< [1/0], Optimize mesh keeping its initial edge sizes */
565 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_noinsert, arg(i,stack,0L)); i++; /*!< [1/0], Avoid/allow point insertion */
566 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_noswap, arg(i,stack,0L)); i++; /*!< [1/0], Avoid/allow edge or face flipping */
567 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_nomove, arg(i,stack,0L)); i++; /*!< [1/0], Avoid/allow point relocation */
568 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_nreg, arg(i,stack,0L)); i++; /*!< [0/1], Disabled/enabled normal regularization */
569 if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_renum, arg(i,stack,0L)); i++; /*!< [1/0], Turn on/off point relocation with Scotch */
570 if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_angleDetection,arg(i,stack,0.)); i++; /*!< [val], Value for angle detection */
571 if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hmin, arg(i,stack,0.)); i++; /*!< [val], Minimal mesh size */
572 if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hmax, arg(i,stack,0.)); i++; /*!< [val], Maximal mesh size */
573 if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hsiz, arg(i,stack,0.)); i++; /*!< [val], Constant mesh size */
574 if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hausd, arg(i,stack,0.)); i++; /*!< [val], Control global Hausdorff distance (on all the boundary surfaces of the mesh) */
575 if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hgrad, arg(i,stack,0.)); i++; /*!< [val], Control gradation */
576 if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_ls, arg(i,stack,0.)); i++; /*!< [val], Value of level-set */
577
578 int ier = MMGS_mmgslib(mesh,sol);
579
580 /*
581 if ( MMG5_saveMesh_centralized(mesh,"test.mesh") != 1 ) {
582 fprintf(stdout,"UNABLE TO SAVE MESH\n");
583 ier = MMG5_STRONGFAILURE;
584 }
585 MMG5_saveMet_centralized(mesh,"test.meshb");
586 */
587
588 MeshS *Th_T = nullptr;
589
590 MMG5_pMesh_to_ffmesh(mesh,Th_T);
591
592 MMGS_Free_all(MMG5_ARG_start,
593 MMG5_ARG_ppMesh,&mesh,MMG5_ARG_ppMet,&sol,
594 MMG5_ARG_end);
595
596 Th_T->BuildGTree();
597
598 Add2StackOfPtr2FreeRC(stack, Th_T);
599 return Th_T;
600 }
601
Load_Init()602 static void Load_Init( ) {
603 if (verbosity) {
604 cout << " load: mmg " << endl;
605 }
606
607 Global.Add("mmg3d", "(", new mmg_ff<Mesh3>);
608 Global.Add("mmgs", "(", new mmg_ff<MeshS>);
609 }
610
611 LOADFUNC(Load_Init)
612