PolyFEM
Loading...
Searching...
No Matches
MMGRemesh.cpp
Go to the documentation of this file.
1// Original source from cellogram (https://github.com/cellogram/cellogram/blob/master/src/cellogram/remesh_adaptive.h)
2// Authors: Tobias Lendenmann, Teseo Schneider, Jérémie Dumas, Marco Tarini
3// License: MIT (https://github.com/cellogram/cellogram/blob/master/LICENSE)
4
5#ifdef POLYFEM_WITH_MMG
6
8#include "MMGRemesh.hpp"
9// #include <cellogram/MeshUtils.h>
10#include <algorithm>
11#include <iomanip>
12#include <cassert>
13#include <vector>
14#include <geogram/basic/attributes.h>
15#include <geogram/mesh/mesh.h>
16#include <geogram/mesh/mesh_io.h>
17#include <igl/boundary_facets.h>
18#include <mmg/libmmg.h>
19
20#include <polyfem/utils/GeometryUtils.hpp>
21#include <polyfem/utils/Logger.hpp>
23//
24// Wrapper for 3D remeshing comes from:
25// https://github.com/mxncr/mmgig
26//
27
28#ifdef WIN32
29typedef unsigned int uint;
30#endif // WIN32
31
32namespace polyfem::mesh
33{
34 namespace
35 {
36 // Flip triangle/tet if there signed volume is negative.
37 bool orient_simplices(const Eigen::MatrixXd &vertices, Eigen::MatrixXi &elements)
38 {
39 for (int e = 0; e < elements.rows(); ++e)
40 {
41 double vol = 0;
42 if (elements.cols() == 3)
43 {
44 vol = utils::triangle_area_2D(
45 vertices.row(elements(e, 0)),
46 vertices.row(elements(e, 1)),
47 vertices.row(elements(e, 2)));
48 }
49 else
50 {
51 assert(elements.cols() == 4);
52 vol = utils::tetrahedron_volume(
53 vertices.row(elements(e, 0)),
54 vertices.row(elements(e, 1)),
55 vertices.row(elements(e, 2)),
56 vertices.row(elements(e, 3)));
57 }
58
59 if (vol == 0)
60 {
61 logger().error("MMG produced a degenerate simplex {}", e);
62 return false;
63 }
64 if (vol < 0)
65 std::swap(elements(e, 1), elements(e, 2));
66 }
67 return true;
68 }
69
70 void to_geogram_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, GEO::Mesh &M)
71 {
72 M.clear();
73 // Setup vertices
74 M.vertices.create_vertices((int)V.rows());
75 for (int i = 0; i < (int)M.vertices.nb(); ++i)
76 {
77 GEO::vec3 &p = M.vertices.point(i);
78 p[0] = V(i, 0);
79 p[1] = V(i, 1);
80 p[2] = (V.cols() == 2 ? 0 : V(i, 2));
81 }
82 // Setup faces
83 if (F.cols() == 3)
84 {
85 M.facets.create_triangles((int)F.rows());
86 }
87 else if (F.cols() == 4)
88 {
89 M.facets.create_quads((int)F.rows());
90 }
91 else
92 {
93 throw std::runtime_error("Mesh faces not supported");
94 }
95 for (int c = 0; c < (int)M.facets.nb(); ++c)
96 {
97 for (int lv = 0; lv < F.cols(); ++lv)
98 {
99 M.facets.set_vertex(c, lv, F(c, lv));
100 }
101 }
102 M.facets.connect();
103 }
104
105 void to_geogram_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXi &T, GEO::Mesh &M)
106 {
107 to_geogram_mesh(V, F, M);
108 if (T.cols() == 4)
109 {
110 M.cells.create_tets((int)T.rows());
111 }
112 else if (T.rows() != 0)
113 {
114 throw std::runtime_error("Mesh cells not supported");
115 }
116 for (int c = 0; c < (int)M.cells.nb(); ++c)
117 {
118 for (int lv = 0; lv < T.cols(); ++lv)
119 {
120 M.cells.set_vertex(c, lv, T(c, lv));
121 }
122 }
123 M.cells.connect();
124 }
125
126 void from_geogram_mesh(const GEO::Mesh &M, Eigen::MatrixXd &V, Eigen::MatrixXi &F, Eigen::MatrixXi &T)
127 {
128 V.resize(M.vertices.nb(), 3);
129 for (int i = 0; i < (int)M.vertices.nb(); ++i)
130 {
131 GEO::vec3 p = M.vertices.point(i);
132 V.row(i) << p[0], p[1], p[2];
133 }
134 assert(M.facets.are_simplices());
135 F.resize(M.facets.nb(), 3);
136 for (int c = 0; c < (int)M.facets.nb(); ++c)
137 {
138 for (int lv = 0; lv < 3; ++lv)
139 {
140 F(c, lv) = M.facets.vertex(c, lv);
141 }
142 }
143 assert(M.cells.are_simplices());
144 T.resize(M.cells.nb(), 4);
145 for (int c = 0; c < (int)M.cells.nb(); ++c)
146 {
147 for (int lv = 0; lv < 4; ++lv)
148 {
149 T(c, lv) = M.cells.vertex(c, lv);
150 }
151 }
152 }
153
154 bool mmg_to_geo(const MMG5_pMesh mmg, GEO::Mesh &M)
155 {
156 logger().trace("converting MMG5_pMesh to GEO::Mesh ...");
157 /* Notes:
158 * - indexing seems to start at 1 in MMG */
159
160 assert(mmg->dim == 2 || mmg->dim == 3);
161 M.clear();
162 M.vertices.create_vertices((uint)mmg->np);
163 M.edges.create_edges((uint)mmg->na);
164 M.facets.create_triangles((uint)mmg->nt);
165 M.cells.create_tets((uint)mmg->ne);
166
167 for (uint v = 0; v < M.vertices.nb(); ++v)
168 {
169 for (uint d = 0; d < (uint)mmg->dim; ++d)
170 {
171 M.vertices.point_ptr(v)[d] = mmg->point[v + 1].c[d];
172 }
173 // Geogram mesh is 3D, init z coordinate to zero for 2D MMG mesh.
174 if (mmg->dim == 2)
175 M.vertices.point_ptr(v)[2] = 0.;
176 }
177 for (uint e = 0; e < M.edges.nb(); ++e)
178 {
179 M.edges.set_vertex(e, 0, (uint)mmg->edge[e + 1].a - 1);
180 M.edges.set_vertex(e, 1, (uint)mmg->edge[e + 1].b - 1);
181 }
182 for (uint t = 0; t < M.facets.nb(); ++t)
183 {
184 M.facets.set_vertex(t, 0, (uint)mmg->tria[t + 1].v[0] - 1);
185 M.facets.set_vertex(t, 1, (uint)mmg->tria[t + 1].v[1] - 1);
186 M.facets.set_vertex(t, 2, (uint)mmg->tria[t + 1].v[2] - 1);
187 }
188 for (uint c = 0; mmg->dim == 3 && c < M.cells.nb(); ++c)
189 {
190 M.cells.set_vertex(c, 0, (uint)mmg->tetra[c + 1].v[0] - 1);
191 M.cells.set_vertex(c, 1, (uint)mmg->tetra[c + 1].v[1] - 1);
192 M.cells.set_vertex(c, 2, (uint)mmg->tetra[c + 1].v[2] - 1);
193 M.cells.set_vertex(c, 3, (uint)mmg->tetra[c + 1].v[3] - 1);
194 }
195 M.facets.connect();
196 M.cells.connect();
197
198 return true;
199 }
200
201 bool geo_to_mmg(const GEO::Mesh &M, MMG5_pMesh &mmg, MMG5_pSol &sol, bool volume_mesh = true)
202 {
203 logger().trace("converting GEO::M to MMG5_pMesh ...");
204 assert(M.vertices.dimension() == 3);
205 if (M.facets.nb() > 0)
206 assert(M.facets.are_simplices());
207 if (M.cells.nb() > 0)
208 assert(M.cells.are_simplices());
209
210 if (volume_mesh)
211 {
212 MMG3D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
213 }
214 else
215 {
216 MMGS_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
217 }
218
219 if (volume_mesh && MMG3D_Set_meshSize(mmg, (int)M.vertices.nb(), (int)M.cells.nb(), 0, /* nb prisms */
220 (int)M.facets.nb(), 0, /* nb quad */
221 (int)M.edges.nb() /* nb edges */
222 )
223 != 1)
224 {
225 logger().error("failed to MMG3D_Set_meshSize");
226 return false;
227 }
228 else if (!volume_mesh && MMGS_Set_meshSize(mmg, (int)M.vertices.nb(), (int)M.facets.nb(), (int)M.edges.nb() /* nb edges */
229 )
230 != 1)
231 {
232 logger().error("failed to MMGS_Set_meshSize");
233 return false;
234 }
235
236 for (uint v = 0; v < (uint)mmg->np; ++v)
237 {
238 for (uint d = 0; d < M.vertices.dimension(); ++d)
239 {
240 mmg->point[v + 1].c[d] = M.vertices.point_ptr(v)[d];
241 }
242 }
243 for (uint e = 0; e < (uint)mmg->na; ++e)
244 {
245 mmg->edge[e + 1].a = (int)M.edges.vertex(e, 0) + 1;
246 mmg->edge[e + 1].b = (int)M.edges.vertex(e, 1) + 1;
247 }
248 for (uint t = 0; t < (uint)mmg->nt; ++t)
249 {
250 mmg->tria[t + 1].v[0] = (int)M.facets.vertex(t, 0) + 1;
251 mmg->tria[t + 1].v[1] = (int)M.facets.vertex(t, 1) + 1;
252 mmg->tria[t + 1].v[2] = (int)M.facets.vertex(t, 2) + 1;
253 }
254 if (volume_mesh)
255 {
256 for (uint c = 0; c < (uint)mmg->ne; ++c)
257 {
258 mmg->tetra[c + 1].v[0] = (int)M.cells.vertex(c, 0) + 1;
259 mmg->tetra[c + 1].v[1] = (int)M.cells.vertex(c, 1) + 1;
260 mmg->tetra[c + 1].v[2] = (int)M.cells.vertex(c, 2) + 1;
261 mmg->tetra[c + 1].v[3] = (int)M.cells.vertex(c, 3) + 1;
262 }
263 }
264
265 if (volume_mesh && MMG3D_Set_solSize(mmg, sol, MMG5_Vertex, (int)M.vertices.nb(), MMG5_Scalar) != 1)
266 {
267 logger().error("failed to MMG3D_Set_solSize");
268 return false;
269 }
270 else if (!volume_mesh && MMGS_Set_solSize(mmg, sol, MMG5_Vertex, (int)M.vertices.nb(), MMG5_Scalar) != 1)
271 {
272 logger().error("failed to MMGS_Set_solSize");
273 return false;
274 }
275 for (uint v = 0; v < M.vertices.nb(); ++v)
276 {
277 sol->m[v + 1] = 1.;
278 }
279 if (volume_mesh && MMG3D_Chk_meshData(mmg, sol) != 1)
280 {
281 logger().error("error in mmg: inconsistent mesh and sol");
282 return false;
283 }
284 else if (!volume_mesh && MMGS_Chk_meshData(mmg, sol) != 1)
285 {
286 logger().error("error in mmg: inconsistent mesh and sol");
287 return false;
288 }
289
290 if (volume_mesh)
291 {
292 MMG3D_Set_handGivenMesh(mmg); /* because we don't use the API functions */
293 }
294
295 return true;
296 }
297
298 bool geo_to_mmg2d(const GEO::Mesh &M, MMG5_pMesh &mmg, MMG5_pSol &sol)
299 {
300 assert(M.vertices.dimension() == 3);
301 assert(M.facets.are_simplices());
302
303 MMG2D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
304 if (!MMG2D_Set_meshSize(mmg, M.vertices.nb(), M.facets.nb(), 0, M.edges.nb()))
305 {
306 logger().error("mmg2d_remesh: failed to allocate the MMG mesh");
307 return false;
308 }
309
310 for (int v = 0; v < mmg->np; ++v)
311 {
312 mmg->point[v + 1].c[0] = M.vertices.point_ptr(v)[0];
313 mmg->point[v + 1].c[1] = M.vertices.point_ptr(v)[1];
314 }
315 for (int e = 0; e < mmg->na; ++e)
316 {
317 // MMG use 1 based indexing.
318 mmg->edge[e + 1].a = M.edges.vertex(e, 0) + 1;
319 mmg->edge[e + 1].b = M.edges.vertex(e, 1) + 1;
320 }
321 for (int t = 0; t < mmg->nt; ++t)
322 {
323 // MMG use 1 based indexing.
324 mmg->tria[t + 1].v[0] = M.facets.vertex(t, 0) + 1;
325 mmg->tria[t + 1].v[1] = M.facets.vertex(t, 1) + 1;
326 mmg->tria[t + 1].v[2] = M.facets.vertex(t, 2) + 1;
327 }
328
329 if (!MMG2D_Set_solSize(mmg, sol, MMG5_Vertex, M.vertices.nb(), MMG5_Scalar))
330 {
331 logger().error("failed to MMG2D_Set_solSize");
332 return false;
333 }
334
335 for (int v = 0; v < M.vertices.nb(); ++v)
336 sol->m[v + 1] = 1.;
337
338 if (!MMG2D_Chk_meshData(mmg, sol))
339 {
340 logger().error("error in mmg: inconsistent mesh and sol");
341 return false;
342 }
343 return true;
344 }
345
346 void mmg2d_free(MMG5_pMesh mmg, MMG5_pSol sol)
347 {
348 MMG2D_Free_all(MMG5_ARG_start,
349 MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
350 }
351
352 void mmg3d_free(MMG5_pMesh mmg, MMG5_pSol sol)
353 {
354 MMG3D_Free_all(MMG5_ARG_start,
355 MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
356 }
357
358 void mmgs_free(MMG5_pMesh mmg, MMG5_pSol sol)
359 {
360 MMGS_Free_all(MMG5_ARG_start,
361 MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
362 }
363
364 bool mmg_wrapper_test_geo2mmg2geo(const GEO::Mesh &M_in, GEO::Mesh &M_out)
365 {
366 MMG5_pMesh mmg = nullptr;
367 MMG5_pSol sol = nullptr;
368 bool ok = geo_to_mmg(M_in, mmg, sol);
369 if (!ok)
370 return false;
371 ok = mmg_to_geo(mmg, M_out);
372 mmg3d_free(mmg, sol);
373 return ok;
374 }
375
376 bool mmg2d_tri_remesh(
377 const GEO::Mesh &M,
378 const std::vector<int> *pinned_vertices,
379 GEO::Mesh &M_out,
380 const MmgOptions &opt)
381 {
382 MMG5_pMesh mesh = nullptr;
383 MMG5_pSol met = nullptr;
384 bool ok = geo_to_mmg2d(M, mesh, met);
385 if (!ok)
386 {
387 logger().error("mmg2d_remesh: failed to convert mesh to MMG5_pMesh");
388 mmg2d_free(mesh, met);
389 return false;
390 }
391
392 // Pin vertices and edges.
393 std::vector<bool> is_pinned(M.vertices.nb(), false);
394 if (pinned_vertices != nullptr)
395 {
396 for (const int v : *pinned_vertices)
397 {
398 assert(v >= 0 && v < M.vertices.nb());
399 is_pinned[v] = true;
400 MMG2D_Set_requiredVertex(mesh, v + 1);
401 }
402 }
403 for (int e = 0; e < M.edges.nb(); ++e)
404 {
405 if (is_pinned[M.edges.vertex(e, 0)]
406 && is_pinned[M.edges.vertex(e, 1)])
407 {
408 MMG2D_Set_requiredEdge(mesh, e + 1);
409 }
410 }
411
412 /* Set remeshing options */
413 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_angleDetection, opt.angle_value);
414 if (opt.enable_anisotropy)
415 {
416 MMG2D_Set_solSize(mesh, met, MMG5_Vertex, 0, MMG5_Tensor);
417 }
418 const bool with_metric = opt.metric_attribute != "no_metric";
419 if (opt.hsiz == 0. || with_metric)
420 {
421 if (with_metric || !opt.optim)
422 {
423 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hmin, opt.hmin);
424 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hmax, opt.hmax);
425 }
426 else
427 {
428 met->np = 0;
429 }
430 }
431 else
432 {
433 met->np = 0;
434 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hsiz, opt.hsiz);
435 }
436 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hausd, opt.hausd);
437 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hgrad, opt.hgrad);
438 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_angle, int(opt.angle_detection));
439 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_noswap, int(opt.noswap));
440 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_noinsert, int(opt.noinsert));
441 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_nomove, int(opt.nomove));
442 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_nosurf, int(opt.nosurf));
443 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_optim, int(opt.optim));
444 if (with_metric)
445 {
446 if (!M.vertices.attributes().is_defined(opt.metric_attribute))
447 {
448 logger().error("mmg2D_remesh: {} is not a vertex attribute, cancel", opt.metric_attribute);
449 mmg2d_free(mesh, met);
450 return false;
451 }
452 GEO::Attribute<double> h_local(M.vertices.attributes(), opt.metric_attribute);
453 for (int v = 0; v < M.vertices.nb(); ++v)
454 {
455 met->m[v + 1] = h_local[v];
456 }
457 }
458
459 int ier = MMG2D_mmg2dlib(mesh, met);
460 if (ier != MMG5_SUCCESS)
461 {
462 logger().error("mmg2d_remesh: failed to remesh");
463 mmg2d_free(mesh, met);
464 return false;
465 }
466
467 ok = mmg_to_geo(mesh, M_out);
468
469 mmg2d_free(mesh, met);
470 return ok;
471 }
472
473 bool mmgs_tri_remesh(const GEO::Mesh &M, GEO::Mesh &M_out, const MmgOptions &opt)
474 {
475 MMG5_pMesh mesh = nullptr;
476 MMG5_pSol met = nullptr;
477 bool ok = geo_to_mmg(M, mesh, met, false);
478 if (!ok)
479 {
480 logger().error("mmgs_remesh: failed to convert mesh to MMG5_pMesh");
481 mmgs_free(mesh, met);
482 return false;
483 }
484
485 /* Set remeshing options */
486 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_angleDetection, opt.angle_value);
487 if (opt.enable_anisotropy)
488 {
489 MMGS_Set_solSize(mesh, met, MMG5_Vertex, 0, MMG5_Tensor);
490 }
491 if (opt.hsiz == 0. || opt.metric_attribute != "no_metric")
492 {
493 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hmin, opt.hmin);
494 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hmax, opt.hmax);
495 }
496 else
497 {
498 met->np = 0;
499 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hsiz, opt.hsiz);
500 }
501 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hausd, opt.hausd);
502 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hgrad, opt.hgrad);
503 MMGS_Set_iparameter(mesh, met, MMGS_IPARAM_angle, int(opt.angle_detection));
504 MMGS_Set_iparameter(mesh, met, MMGS_IPARAM_noswap, int(opt.noswap));
505 MMGS_Set_iparameter(mesh, met, MMGS_IPARAM_noinsert, int(opt.noinsert));
506 MMGS_Set_iparameter(mesh, met, MMGS_IPARAM_nomove, int(opt.nomove));
507 if (opt.metric_attribute != "no_metric")
508 {
509 if (!M.vertices.attributes().is_defined(opt.metric_attribute))
510 {
511 logger().error("mmgs_remesh: {} is not a vertex attribute, cancel", opt.metric_attribute);
512 return false;
513 }
514 GEO::Attribute<double> h_local(M.vertices.attributes(), opt.metric_attribute);
515 for (uint v = 0; v < M.vertices.nb(); ++v)
516 {
517 met->m[v + 1] = h_local[v];
518 }
519 }
520
521 int ier = MMGS_mmgslib(mesh, met);
522 if (ier != MMG5_SUCCESS)
523 {
524 logger().error("mmgs_remesh: failed to remesh");
525 mmgs_free(mesh, met);
526 return false;
527 }
528
529 ok = mmg_to_geo(mesh, M_out);
530
531 mmgs_free(mesh, met);
532 return ok;
533 }
534
535 bool mmg3d_tet_remesh(
536 const GEO::Mesh &M,
537 const std::vector<int> *pinned_vertices,
538 GEO::Mesh &M_out,
539 const MmgOptions &opt)
540 {
541 MMG5_pMesh mesh = nullptr;
542 MMG5_pSol met = nullptr;
543 bool ok = geo_to_mmg(M, mesh, met, true);
544 if (!ok)
545 {
546 logger().error("mmg3d_remesh: failed to convert mesh to MMG5_pMesh");
547 mmg3d_free(mesh, met);
548 return false;
549 }
550
551 // Pin vertices and triangles.
552 std::vector<bool> is_pinned(M.vertices.nb(), false);
553 if (pinned_vertices != nullptr)
554 {
555 for (const int v : *pinned_vertices)
556 {
557 assert(v >= 0 && v < M.vertices.nb());
558 is_pinned[v] = true;
559 MMG3D_Set_requiredVertex(mesh, v + 1);
560 }
561 }
562 for (int f = 0; f < M.facets.nb(); ++f)
563 {
564 if (is_pinned[M.facets.vertex(f, 0)]
565 && is_pinned[M.facets.vertex(f, 1)]
566 && is_pinned[M.facets.vertex(f, 2)])
567 {
568 MMG3D_Set_requiredTriangle(mesh, f + 1);
569 }
570 }
571
572 /* Set remeshing options */
573 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_angleDetection, opt.angle_value);
574 if (opt.enable_anisotropy)
575 {
576 MMG3D_Set_solSize(mesh, met, MMG5_Vertex, 0, MMG5_Tensor);
577 }
578 const bool with_metric = opt.metric_attribute != "no_metric";
579 if (opt.hsiz == 0. || with_metric)
580 {
581 if (with_metric || !opt.optim)
582 {
583 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hmin, opt.hmin);
584 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hmax, opt.hmax);
585 }
586 else
587 {
588 met->np = 0;
589 }
590 }
591 else
592 {
593 met->np = 0;
594 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hsiz, opt.hsiz);
595 }
596 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hausd, opt.hausd);
597 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hgrad, opt.hgrad);
598 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_angle, int(opt.angle_detection));
599 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_noswap, int(opt.noswap));
600 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_noinsert, int(opt.noinsert));
601 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_nomove, int(opt.nomove));
602 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_nosurf, int(opt.nosurf));
603 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_opnbdy, int(opt.opnbdy));
604 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_optim, int(opt.optim));
605 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_optimLES, int(opt.optimLES));
606 if (with_metric)
607 {
608 if (!M.vertices.attributes().is_defined(opt.metric_attribute))
609 {
610 logger().error("mmg3D_remesh: {} is not a vertex attribute, cancel", opt.metric_attribute);
611 mmg3d_free(mesh, met);
612 return false;
613 }
614 GEO::Attribute<double> h_local(M.vertices.attributes(), opt.metric_attribute);
615 for (int v = 0; v < M.vertices.nb(); ++v)
616 {
617 met->m[v + 1] = h_local[v];
618 }
619 }
620
621 int ier = MMG3D_mmg3dlib(mesh, met);
622 if (ier != MMG5_SUCCESS)
623 {
624 logger().error("mmg3d_remesh: failed to remesh");
625 mmg3d_free(mesh, met);
626 return false;
627 }
628
629 ok = mmg_to_geo(mesh, M_out);
630
631 mmg3d_free(mesh, met);
632 return ok;
633 }
634
635 bool mmg3d_extract_iso(const GEO::Mesh &M, GEO::Mesh &M_out, const MmgOptions &opt)
636 {
637 if (!opt.level_set || opt.ls_attribute == "no_ls" || !M.vertices.attributes().is_defined(opt.ls_attribute))
638 {
639 logger().error("mmg3D_iso: {} is not a vertex attribute, cancel", opt.ls_attribute);
640 return false;
641 }
642 if (opt.angle_detection)
643 {
644 logger().warn("mmg3D_iso: angle_detection shoud probably be disabled because level set functions are smooth");
645 }
646
647 MMG5_pMesh mesh = nullptr;
648 MMG5_pSol met = nullptr;
649 bool ok = geo_to_mmg(M, mesh, met, true);
650 if (!ok)
651 {
652 logger().error("mmg3d_remesh: failed to convert mesh to MMG5_pMesh");
653 mmg3d_free(mesh, met);
654 return false;
655 }
656 GEO::Attribute<double> ls(M.vertices.attributes(), opt.ls_attribute);
657 for (uint v = 0; v < M.vertices.nb(); ++v)
658 {
659 met->m[v + 1] = ls[v];
660 }
661
662 /* Flag border for future deletion */
663 // std::vector<bool> on_border(M.vertices.nb(), false);
664 // for (index_t t = 0; t < M.cells.nb(); ++t) {
665 // for (index_t lf = 0; lf < M.cells.nb_facets(t); ++lf) {
666 // if (M.cells.adjacent(t,lf) != GEO::NO_CELL) continue;
667 // for (index_t lv = 0; lv < M.cells.facet_nb_vertices(t,lf); ++lv) {
668 // on_border[M.cells.facet_vertex(t,lf,lv)] = true;
669 // }
670 // }
671 // }
672
673 /* Set remeshing options */
674 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_iso, 1);
675 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_ls, opt.ls_value);
676 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_angleDetection, opt.angle_value);
677 if (opt.hsiz == 0.)
678 {
679 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hmin, opt.hmin);
680 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hmax, opt.hmax);
681 }
682 else
683 {
684 logger().error("mmg3d_iso: should not use hsiz parameter for level set mode");
685 mmg3d_free(mesh, met);
686 return false;
687 }
688 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hausd, opt.hausd);
689 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hgrad, opt.hgrad);
690 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_angle, int(opt.angle_detection));
691 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_noswap, int(opt.noswap));
692 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_noinsert, 1);
693 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_nomove, 1);
694 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_nosurf, 1);
695
696 // TODO: Check this is correct
697 // Used to be
698 // int ier = MMG3D_mmg3dls(mesh, met);
699 int ier = MMG3D_mmg3dls(mesh, met, nullptr);
700 if (ier != MMG5_SUCCESS)
701 {
702 logger().error("mmg3d_iso: failed to remesh isovalue");
703 mmg3d_free(mesh, met);
704 return false;
705 }
706
707 /* Convert back */
708 ok = mmg_to_geo(mesh, M_out);
709 GEO::Attribute<double> ls_out(M_out.vertices.attributes(), opt.ls_attribute);
710 for (uint v = 0; v < M_out.vertices.nb(); ++v)
711 {
712 ls_out[v] = met->m[v + 1];
713 }
714 /* Extract only the border */
715 // M_out.cells.clear(false,false);
716 // M_out.vertices.remove_isolated();
717 // GEO::vector<index_t> to_del(M_out.facets.nb(), 0);
718 // for (index_t f = 0; f < M_out.facets.nb(); ++f) {
719 // double d = 0;
720 // bool f_on_border = true;
721 // for (index_t lv = 0; lv < M_out.facets.nb_vertices(f); ++lv) {
722 // d = geo_max(d,std::abs(ls_out[M_out.facets.vertex(f,0)] - opt.ls_value));
723 // if (M_out.facets.vertex(f,lv) < M.vertices.nb()) {
724 // if (!on_border[M_out.facets.vertex(f,lv)]) f_on_border = false;
725 // } else {
726 // f_on_border = false;
727 // }
728 // }
729 // // if (d > 1.1 * opt.hmin) {
730 // // to_del[f] = 1;
731 // // }
732 // if (f_on_border) to_del[f] = 1;
733 // }
734 // M_out.facets.delete_elements(to_del, true);
735
736 mmg3d_free(mesh, met);
737 return ok;
738 }
739
740 } // anonymous namespace
741
743
744 bool remesh_2d(
745 const Eigen::MatrixXd &V,
746 const Eigen::MatrixXi &F,
747 Eigen::MatrixXd &OV,
748 Eigen::MatrixXi &OF,
749 MmgOptions opt,
750 const std::vector<int> *pinned_vertices)
751 {
752 assert(V.cols() == 2);
753 assert(F.cols() == 3);
754 GEO::Mesh M, M_out;
755 to_geogram_mesh(V, F, M);
756 if (!mmg2d_tri_remesh(M, pinned_vertices, M_out, opt))
757 return false;
758 Eigen::MatrixXi OT;
759 from_geogram_mesh(M_out, OV, OF, OT);
760 // Geogram mesh is 3d only, truncate unused third column.
761 OV.conservativeResize(OV.rows(), 2);
762 // MMG does not preserve orientation. Ensure singed area is positive.
763 return orient_simplices(OV, OF);
764 }
765
766 void remesh_adaptive_2d(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::VectorXd &S,
767 Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, MmgOptions opt)
768 {
769 assert(V.cols() == 2 || V.cols() == 3);
770 assert(V.rows() == S.size());
771 assert(F.cols() == 3);
772 GEO::Mesh M, M_out;
773 to_geogram_mesh(V, F, M);
774
775 opt.metric_attribute = "scalar";
776 GEO::Attribute<double> scalar(M.vertices.attributes(), opt.metric_attribute);
777 for (int v = 0; v < M.vertices.nb(); ++v)
778 scalar[v] = S(v);
779
780 if (!mmg2d_tri_remesh(M, nullptr, M_out, opt))
781 return;
782 Eigen::MatrixXi OT;
783 from_geogram_mesh(M_out, OV, OF, OT);
784 OV.conservativeResize(OV.rows(), 2);
785 orient_simplices(OV, OF);
786 }
787
788 bool remesh_3d(
789 const Eigen::MatrixXd &V,
790 const Eigen::MatrixXi &T,
791 Eigen::MatrixXd &OV,
792 Eigen::MatrixXi &OF,
793 Eigen::MatrixXi &OT,
794 MmgOptions opt,
795 const std::vector<int> *pinned_vertices)
796 {
797 assert(V.cols() == 3);
798 assert(T.cols() == 4);
799 GEO::Mesh M, M_out;
800 Eigen::MatrixXi boundary;
801 igl::boundary_facets(T, boundary);
802 to_geogram_mesh(V, boundary, T, M);
803 if (!mmg3d_tet_remesh(M, pinned_vertices, M_out, opt))
804 return false;
805 from_geogram_mesh(M_out, OV, OF, OT);
806 // MMG does not preserve orientation. Ensure singed area is positive.
807 return orient_simplices(OV, OT);
808 }
809
810 void remesh_adaptive_3d(const Eigen::MatrixXd &V, const Eigen::MatrixXi &T, const Eigen::VectorXd &S,
811 Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::MatrixXi &OT, MmgOptions opt)
812 {
813 assert(V.cols() == 3);
814 assert(V.rows() == S.size());
815 GEO::Mesh M, M_out;
816 to_geogram_mesh(V, Eigen::MatrixXi(0, 3), T, M);
817
818 // Remeshing options
819 opt.metric_attribute = "scalar";
820
821 GEO::Attribute<double> scalar(M.vertices.attributes(), opt.metric_attribute);
822 for (int v = 0; v < M.vertices.nb(); ++v)
823 {
824 scalar[v] = S(v);
825 }
826
827 // Remesh volume
828 if (!mmg3d_tet_remesh(M, nullptr, M_out, opt))
829 return;
830
831 // Convert output
832 from_geogram_mesh(M_out, OV, OF, OT);
833 orient_simplices(OV, OT);
834 }
835
836} // namespace polyfem::mesh
837
838#endif
V
int V
Definition AdjointNLProblem.cpp:48
eigs.M
M
Definition eigs.py:94
elasticity_rhs.f
f
Definition elasticity_rhs.py:173
generate_rotation_mtx.e
list e
Definition generate_rotation_mtx.py:49
generate_rotation_mtx.T
T
Definition generate_rotation_mtx.py:43
p_bases.vertices
list vertices
Definition p_bases.py:238
polyfem::mesh::to_geogram_mesh
void to_geogram_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, GEO::Mesh &M)
Converts a triangle mesh to a Geogram mesh.
Definition MeshUtils.cpp:578
polyfem::mesh::from_geogram_mesh
void from_geogram_mesh(const GEO::Mesh &M, Eigen::MatrixXd &V, Eigen::MatrixXi &F, Eigen::MatrixXi &T)
Extract simplices from a Geogram mesh.
Definition MeshUtils.cpp:642
polyfem::utils::tetrahedron_volume
double tetrahedron_volume(const Eigen::Vector3d &a, const Eigen::Vector3d &b, const Eigen::Vector3d &c, const Eigen::Vector3d &d)
Compute the signed volume of a tetrahedron defined by four points.
Definition GeometryUtils.cpp:54
polyfem::utils::triangle_area_2D
double triangle_area_2D(const Eigen::Vector2d &a, const Eigen::Vector2d &b, const Eigen::Vector2d &c)
Compute the signed area of a 2D triangle defined by three points.
Definition GeometryUtils.cpp:25
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44