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 <iomanip>
11#include <cassert>
12#include <geogram/basic/attributes.h>
13#include <geogram/mesh/mesh.h>
14#include <geogram/mesh/mesh_io.h>
15#include <mmg/libmmg.h>
16
17#include <polyfem/utils/Logger.hpp>
19//
20// Wrapper for 3D remeshing comes from:
21// https://github.com/mxncr/mmgig
22//
23
24#ifdef WIN32
25typedef unsigned int uint;
26#endif // WIN32
27
28namespace polyfem::mesh
29{
30 namespace
31 {
32 void to_geogram_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, GEO::Mesh &M)
33 {
34 M.clear();
35 // Setup vertices
36 M.vertices.create_vertices((int)V.rows());
37 for (int i = 0; i < (int)M.vertices.nb(); ++i)
38 {
39 GEO::vec3 &p = M.vertices.point(i);
40 p[0] = V(i, 0);
41 p[1] = V(i, 1);
42 p[2] = (V.cols() == 2 ? 0 : V(i, 2));
43 }
44 // Setup faces
45 if (F.cols() == 3)
46 {
47 M.facets.create_triangles((int)F.rows());
48 }
49 else if (F.cols() == 4)
50 {
51 M.facets.create_quads((int)F.rows());
52 }
53 else
54 {
55 throw std::runtime_error("Mesh faces not supported");
56 }
57 for (int c = 0; c < (int)M.facets.nb(); ++c)
58 {
59 for (int lv = 0; lv < F.cols(); ++lv)
60 {
61 M.facets.set_vertex(c, lv, F(c, lv));
62 }
63 }
64 M.facets.connect();
65 }
66
67 void to_geogram_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXi &T, GEO::Mesh &M)
68 {
69 to_geogram_mesh(V, F, M);
70 if (T.cols() == 4)
71 {
72 M.cells.create_tets((int)T.rows());
73 }
74 else if (T.rows() != 0)
75 {
76 throw std::runtime_error("Mesh cells not supported");
77 }
78 for (int c = 0; c < (int)M.cells.nb(); ++c)
79 {
80 for (int lv = 0; lv < T.cols(); ++lv)
81 {
82 M.cells.set_vertex(c, lv, T(c, lv));
83 }
84 }
85 M.cells.connect();
86 }
87
88 void from_geogram_mesh(const GEO::Mesh &M, Eigen::MatrixXd &V, Eigen::MatrixXi &F, Eigen::MatrixXi &T)
89 {
90 V.resize(M.vertices.nb(), 3);
91 for (int i = 0; i < (int)M.vertices.nb(); ++i)
92 {
93 GEO::vec3 p = M.vertices.point(i);
94 V.row(i) << p[0], p[1], p[2];
95 }
96 assert(M.facets.are_simplices());
97 F.resize(M.facets.nb(), 3);
98 for (int c = 0; c < (int)M.facets.nb(); ++c)
99 {
100 for (int lv = 0; lv < 3; ++lv)
101 {
102 F(c, lv) = M.facets.vertex(c, lv);
103 }
104 }
105 assert(M.cells.are_simplices());
106 T.resize(M.cells.nb(), 4);
107 for (int c = 0; c < (int)M.cells.nb(); ++c)
108 {
109 for (int lv = 0; lv < 4; ++lv)
110 {
111 T(c, lv) = M.cells.vertex(c, lv);
112 }
113 }
114 }
115
116 bool mmg_to_geo(const MMG5_pMesh mmg, GEO::Mesh &M)
117 {
118 logger().trace("converting MMG5_pMesh to GEO::Mesh ...");
119 /* Notes:
120 * - indexing seems to start at 1 in MMG */
121
122 assert(mmg->dim == 3);
123 M.clear();
124 M.vertices.create_vertices((uint)mmg->np);
125 M.edges.create_edges((uint)mmg->na);
126 M.facets.create_triangles((uint)mmg->nt);
127 M.cells.create_tets((uint)mmg->ne);
128
129 for (uint v = 0; v < M.vertices.nb(); ++v)
130 {
131 for (uint d = 0; d < (uint)mmg->dim; ++d)
132 {
133 M.vertices.point_ptr(v)[d] = mmg->point[v + 1].c[d];
134 }
135 }
136 for (uint e = 0; e < M.edges.nb(); ++e)
137 {
138 M.edges.set_vertex(e, 0, (uint)mmg->edge[e + 1].a - 1);
139 M.edges.set_vertex(e, 1, (uint)mmg->edge[e + 1].b - 1);
140 }
141 for (uint t = 0; t < M.facets.nb(); ++t)
142 {
143 M.facets.set_vertex(t, 0, (uint)mmg->tria[t + 1].v[0] - 1);
144 M.facets.set_vertex(t, 1, (uint)mmg->tria[t + 1].v[1] - 1);
145 M.facets.set_vertex(t, 2, (uint)mmg->tria[t + 1].v[2] - 1);
146 }
147 for (uint c = 0; c < M.cells.nb(); ++c)
148 {
149 M.cells.set_vertex(c, 0, (uint)mmg->tetra[c + 1].v[0] - 1);
150 M.cells.set_vertex(c, 1, (uint)mmg->tetra[c + 1].v[1] - 1);
151 M.cells.set_vertex(c, 2, (uint)mmg->tetra[c + 1].v[2] - 1);
152 M.cells.set_vertex(c, 3, (uint)mmg->tetra[c + 1].v[3] - 1);
153 }
154 M.facets.connect();
155 M.cells.connect();
156
157 return true;
158 }
159
160 bool geo_to_mmg(const GEO::Mesh &M, MMG5_pMesh &mmg, MMG5_pSol &sol, bool volume_mesh = true)
161 {
162 logger().trace("converting GEO::M to MMG5_pMesh ...");
163 assert(M.vertices.dimension() == 3);
164 if (M.facets.nb() > 0)
165 assert(M.facets.are_simplices());
166 if (M.cells.nb() > 0)
167 assert(M.cells.are_simplices());
168
169 if (volume_mesh)
170 {
171 MMG3D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
172 }
173 else
174 {
175 MMGS_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
176 }
177
178 if (volume_mesh && MMG3D_Set_meshSize(mmg, (int)M.vertices.nb(), (int)M.cells.nb(), 0, /* nb prisms */
179 (int)M.facets.nb(), 0, /* nb quad */
180 (int)M.edges.nb() /* nb edges */
181 )
182 != 1)
183 {
184 logger().error("failed to MMG3D_Set_meshSize");
185 return false;
186 }
187 else if (!volume_mesh && MMGS_Set_meshSize(mmg, (int)M.vertices.nb(), (int)M.facets.nb(), (int)M.edges.nb() /* nb edges */
188 )
189 != 1)
190 {
191 logger().error("failed to MMGS_Set_meshSize");
192 return false;
193 }
194
195 for (uint v = 0; v < (uint)mmg->np; ++v)
196 {
197 for (uint d = 0; d < M.vertices.dimension(); ++d)
198 {
199 mmg->point[v + 1].c[d] = M.vertices.point_ptr(v)[d];
200 }
201 }
202 for (uint e = 0; e < (uint)mmg->na; ++e)
203 {
204 mmg->edge[e + 1].a = (int)M.edges.vertex(e, 0) + 1;
205 mmg->edge[e + 1].b = (int)M.edges.vertex(e, 1) + 1;
206 }
207 for (uint t = 0; t < (uint)mmg->nt; ++t)
208 {
209 mmg->tria[t + 1].v[0] = (int)M.facets.vertex(t, 0) + 1;
210 mmg->tria[t + 1].v[1] = (int)M.facets.vertex(t, 1) + 1;
211 mmg->tria[t + 1].v[2] = (int)M.facets.vertex(t, 2) + 1;
212 }
213 if (volume_mesh)
214 {
215 for (uint c = 0; c < (uint)mmg->ne; ++c)
216 {
217 mmg->tetra[c + 1].v[0] = (int)M.cells.vertex(c, 0) + 1;
218 mmg->tetra[c + 1].v[1] = (int)M.cells.vertex(c, 1) + 1;
219 mmg->tetra[c + 1].v[2] = (int)M.cells.vertex(c, 2) + 1;
220 mmg->tetra[c + 1].v[3] = (int)M.cells.vertex(c, 3) + 1;
221 }
222 }
223
224 if (volume_mesh && MMG3D_Set_solSize(mmg, sol, MMG5_Vertex, (int)M.vertices.nb(), MMG5_Scalar) != 1)
225 {
226 logger().error("failed to MMG3D_Set_solSize");
227 return false;
228 }
229 else if (!volume_mesh && MMGS_Set_solSize(mmg, sol, MMG5_Vertex, (int)M.vertices.nb(), MMG5_Scalar) != 1)
230 {
231 logger().error("failed to MMGS_Set_solSize");
232 return false;
233 }
234 for (uint v = 0; v < M.vertices.nb(); ++v)
235 {
236 sol->m[v + 1] = 1.;
237 }
238 if (volume_mesh && MMG3D_Chk_meshData(mmg, sol) != 1)
239 {
240 logger().error("error in mmg: inconsistent mesh and sol");
241 return false;
242 }
243 else if (!volume_mesh && MMGS_Chk_meshData(mmg, sol) != 1)
244 {
245 logger().error("error in mmg: inconsistent mesh and sol");
246 return false;
247 }
248
249 if (volume_mesh)
250 {
251 MMG3D_Set_handGivenMesh(mmg); /* because we don't use the API functions */
252 }
253
254 return true;
255 }
256
257 void mmg2d_free(MMG5_pMesh mmg, MMG5_pSol sol)
258 {
259 MMG2D_Free_all(MMG5_ARG_start,
260 MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
261 }
262
263 void mmg3d_free(MMG5_pMesh mmg, MMG5_pSol sol)
264 {
265 MMG3D_Free_all(MMG5_ARG_start,
266 MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
267 }
268
269 void mmgs_free(MMG5_pMesh mmg, MMG5_pSol sol)
270 {
271 MMGS_Free_all(MMG5_ARG_start,
272 MMG5_ARG_ppMesh, &mmg, MMG5_ARG_ppMet, &sol, MMG5_ARG_end);
273 }
274
275 bool mmg_wrapper_test_geo2mmg2geo(const GEO::Mesh &M_in, GEO::Mesh &M_out)
276 {
277 MMG5_pMesh mmg = nullptr;
278 MMG5_pSol sol = nullptr;
279 bool ok = geo_to_mmg(M_in, mmg, sol);
280 if (!ok)
281 return false;
282 ok = mmg_to_geo(mmg, M_out);
283 mmg3d_free(mmg, sol);
284 return ok;
285 }
286
287 bool mmg2d_tri_remesh(const GEO::Mesh &M, GEO::Mesh &M_out, const MmgOptions &opt)
288 {
289 MMG5_pMesh mesh = nullptr;
290 MMG5_pSol met = nullptr;
291 bool ok = geo_to_mmg(M, mesh, met, false);
292 if (!ok)
293 {
294 logger().error("mmg2d_remesh: failed to convert mesh to MMG5_pMesh");
295 mmg2d_free(mesh, met);
296 return false;
297 }
298
299 /* Set remeshing options */
300 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_angleDetection, opt.angle_value);
301 if (opt.enable_anisotropy)
302 {
303 MMG2D_Set_solSize(mesh, met, MMG5_Vertex, 0, MMG5_Tensor);
304 }
305 if (opt.hsiz == 0. || opt.metric_attribute != "no_metric")
306 {
307 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hmin, opt.hmin);
308 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hmax, opt.hmax);
309 }
310 else
311 {
312 met->np = 0;
313 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hsiz, opt.hsiz);
314 }
315 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hausd, opt.hausd);
316 MMG2D_Set_dparameter(mesh, met, MMG2D_DPARAM_hgrad, opt.hgrad);
317 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_angle, int(opt.angle_detection));
318 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_noswap, int(opt.noswap));
319 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_noinsert, int(opt.noinsert));
320 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_nomove, int(opt.nomove));
321 MMG2D_Set_iparameter(mesh, met, MMG2D_IPARAM_nosurf, int(opt.nosurf));
322 if (opt.metric_attribute != "no_metric")
323 {
324 if (!M.vertices.attributes().is_defined(opt.metric_attribute))
325 {
326 logger().error("mmg2D_remesh: {} is not a vertex attribute, cancel", opt.metric_attribute);
327 return false;
328 }
329 GEO::Attribute<double> h_local(M.vertices.attributes(), opt.metric_attribute);
330 for (uint v = 0; v < M.vertices.nb(); ++v)
331 {
332 met->m[v + 1] = h_local[v];
333 }
334 }
335
336 int ier = MMG2D_mmg2dlib(mesh, met);
337 if (ier != MMG5_SUCCESS)
338 {
339 logger().error("mmg2d_remesh: failed to remesh");
340 mmg2d_free(mesh, met);
341 return false;
342 }
343
344 ok = mmg_to_geo(mesh, M_out);
345
346 mmg2d_free(mesh, met);
347 return ok;
348 }
349
350 bool mmgs_tri_remesh(const GEO::Mesh &M, GEO::Mesh &M_out, const MmgOptions &opt)
351 {
352 MMG5_pMesh mesh = nullptr;
353 MMG5_pSol met = nullptr;
354 bool ok = geo_to_mmg(M, mesh, met, false);
355 if (!ok)
356 {
357 logger().error("mmgs_remesh: failed to convert mesh to MMG5_pMesh");
358 mmgs_free(mesh, met);
359 return false;
360 }
361
362 /* Set remeshing options */
363 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_angleDetection, opt.angle_value);
364 if (opt.enable_anisotropy)
365 {
366 MMGS_Set_solSize(mesh, met, MMG5_Vertex, 0, MMG5_Tensor);
367 }
368 if (opt.hsiz == 0. || opt.metric_attribute != "no_metric")
369 {
370 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hmin, opt.hmin);
371 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hmax, opt.hmax);
372 }
373 else
374 {
375 met->np = 0;
376 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hsiz, opt.hsiz);
377 }
378 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hausd, opt.hausd);
379 MMGS_Set_dparameter(mesh, met, MMGS_DPARAM_hgrad, opt.hgrad);
380 MMGS_Set_iparameter(mesh, met, MMGS_IPARAM_angle, int(opt.angle_detection));
381 MMGS_Set_iparameter(mesh, met, MMGS_IPARAM_noswap, int(opt.noswap));
382 MMGS_Set_iparameter(mesh, met, MMGS_IPARAM_noinsert, int(opt.noinsert));
383 MMGS_Set_iparameter(mesh, met, MMGS_IPARAM_nomove, int(opt.nomove));
384 if (opt.metric_attribute != "no_metric")
385 {
386 if (!M.vertices.attributes().is_defined(opt.metric_attribute))
387 {
388 logger().error("mmgs_remesh: {} is not a vertex attribute, cancel", opt.metric_attribute);
389 return false;
390 }
391 GEO::Attribute<double> h_local(M.vertices.attributes(), opt.metric_attribute);
392 for (uint v = 0; v < M.vertices.nb(); ++v)
393 {
394 met->m[v + 1] = h_local[v];
395 }
396 }
397
398 int ier = MMGS_mmgslib(mesh, met);
399 if (ier != MMG5_SUCCESS)
400 {
401 logger().error("mmgs_remesh: failed to remesh");
402 mmgs_free(mesh, met);
403 return false;
404 }
405
406 ok = mmg_to_geo(mesh, M_out);
407
408 mmgs_free(mesh, met);
409 return ok;
410 }
411
412 bool mmg3d_tet_remesh(const GEO::Mesh &M, GEO::Mesh &M_out, const MmgOptions &opt)
413 {
414 MMG5_pMesh mesh = nullptr;
415 MMG5_pSol met = nullptr;
416 bool ok = geo_to_mmg(M, mesh, met, true);
417 if (!ok)
418 {
419 logger().error("mmg3d_remesh: failed to convert mesh to MMG5_pMesh");
420 mmg3d_free(mesh, met);
421 return false;
422 }
423
424 /* Set remeshing options */
425 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_angleDetection, opt.angle_value);
426 if (opt.enable_anisotropy)
427 {
428 MMG3D_Set_solSize(mesh, met, MMG5_Vertex, 0, MMG5_Tensor);
429 }
430 if (opt.hsiz == 0. || opt.metric_attribute != "no_metric")
431 {
432 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hmin, opt.hmin);
433 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hmax, opt.hmax);
434 }
435 else
436 {
437 met->np = 0;
438 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hsiz, opt.hsiz);
439 }
440 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hausd, opt.hausd);
441 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hgrad, opt.hgrad);
442 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_angle, int(opt.angle_detection));
443 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_noswap, int(opt.noswap));
444 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_noinsert, int(opt.noinsert));
445 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_nomove, int(opt.nomove));
446 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_nosurf, int(opt.nosurf));
447 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_opnbdy, int(opt.opnbdy));
448 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_optim, int(opt.optim));
449 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_optimLES, int(opt.optimLES));
450 if (opt.metric_attribute != "no_metric")
451 {
452 if (!M.vertices.attributes().is_defined(opt.metric_attribute))
453 {
454 logger().error("mmg3D_remesh: {} is not a vertex attribute, cancel", opt.metric_attribute);
455 return false;
456 }
457 GEO::Attribute<double> h_local(M.vertices.attributes(), opt.metric_attribute);
458 for (uint v = 0; v < M.vertices.nb(); ++v)
459 {
460 met->m[v + 1] = h_local[v];
461 }
462 }
463
464 int ier = MMG3D_mmg3dlib(mesh, met);
465 if (ier != MMG5_SUCCESS)
466 {
467 logger().error("mmg3d_remesh: failed to remesh");
468 mmg3d_free(mesh, met);
469 return false;
470 }
471
472 ok = mmg_to_geo(mesh, M_out);
473
474 mmg3d_free(mesh, met);
475 return ok;
476 }
477
478 bool mmg3d_extract_iso(const GEO::Mesh &M, GEO::Mesh &M_out, const MmgOptions &opt)
479 {
480 if (!opt.level_set || opt.ls_attribute == "no_ls" || !M.vertices.attributes().is_defined(opt.ls_attribute))
481 {
482 logger().error("mmg3D_iso: {} is not a vertex attribute, cancel", opt.ls_attribute);
483 return false;
484 }
485 if (opt.angle_detection)
486 {
487 logger().warn("mmg3D_iso: angle_detection shoud probably be disabled because level set functions are smooth");
488 }
489
490 MMG5_pMesh mesh = nullptr;
491 MMG5_pSol met = nullptr;
492 bool ok = geo_to_mmg(M, mesh, met, true);
493 if (!ok)
494 {
495 logger().error("mmg3d_remesh: failed to convert mesh to MMG5_pMesh");
496 mmg3d_free(mesh, met);
497 return false;
498 }
499 GEO::Attribute<double> ls(M.vertices.attributes(), opt.ls_attribute);
500 for (uint v = 0; v < M.vertices.nb(); ++v)
501 {
502 met->m[v + 1] = ls[v];
503 }
504
505 /* Flag border for future deletion */
506 // std::vector<bool> on_border(M.vertices.nb(), false);
507 // for (index_t t = 0; t < M.cells.nb(); ++t) {
508 // for (index_t lf = 0; lf < M.cells.nb_facets(t); ++lf) {
509 // if (M.cells.adjacent(t,lf) != GEO::NO_CELL) continue;
510 // for (index_t lv = 0; lv < M.cells.facet_nb_vertices(t,lf); ++lv) {
511 // on_border[M.cells.facet_vertex(t,lf,lv)] = true;
512 // }
513 // }
514 // }
515
516 /* Set remeshing options */
517 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_iso, 1);
518 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_ls, opt.ls_value);
519 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_angleDetection, opt.angle_value);
520 if (opt.hsiz == 0.)
521 {
522 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hmin, opt.hmin);
523 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hmax, opt.hmax);
524 }
525 else
526 {
527 logger().error("mmg3d_iso: should not use hsiz parameter for level set mode");
528 mmg3d_free(mesh, met);
529 return false;
530 }
531 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hausd, opt.hausd);
532 MMG3D_Set_dparameter(mesh, met, MMG3D_DPARAM_hgrad, opt.hgrad);
533 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_angle, int(opt.angle_detection));
534 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_noswap, int(opt.noswap));
535 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_noinsert, 1);
536 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_nomove, 1);
537 MMG3D_Set_iparameter(mesh, met, MMG3D_IPARAM_nosurf, 1);
538
539 // TODO: Check this is correct
540 // Used to be
541 // int ier = MMG3D_mmg3dls(mesh, met);
542 int ier = MMG3D_mmg3dls(mesh, met, nullptr);
543 if (ier != MMG5_SUCCESS)
544 {
545 logger().error("mmg3d_iso: failed to remesh isovalue");
546 mmg3d_free(mesh, met);
547 return false;
548 }
549
550 /* Convert back */
551 ok = mmg_to_geo(mesh, M_out);
552 GEO::Attribute<double> ls_out(M_out.vertices.attributes(), opt.ls_attribute);
553 for (uint v = 0; v < M_out.vertices.nb(); ++v)
554 {
555 ls_out[v] = met->m[v + 1];
556 }
557 /* Extract only the border */
558 // M_out.cells.clear(false,false);
559 // M_out.vertices.remove_isolated();
560 // GEO::vector<index_t> to_del(M_out.facets.nb(), 0);
561 // for (index_t f = 0; f < M_out.facets.nb(); ++f) {
562 // double d = 0;
563 // bool f_on_border = true;
564 // for (index_t lv = 0; lv < M_out.facets.nb_vertices(f); ++lv) {
565 // d = geo_max(d,std::abs(ls_out[M_out.facets.vertex(f,0)] - opt.ls_value));
566 // if (M_out.facets.vertex(f,lv) < M.vertices.nb()) {
567 // if (!on_border[M_out.facets.vertex(f,lv)]) f_on_border = false;
568 // } else {
569 // f_on_border = false;
570 // }
571 // }
572 // // if (d > 1.1 * opt.hmin) {
573 // // to_del[f] = 1;
574 // // }
575 // if (f_on_border) to_del[f] = 1;
576 // }
577 // M_out.facets.delete_elements(to_del, true);
578
579 mmg3d_free(mesh, met);
580 return ok;
581 }
582
583 } // anonymous namespace
584
586
587 void remesh_adaptive_2d(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::VectorXd &S,
588 Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, MmgOptions opt)
589 {
590 assert(V.cols() == 2 || V.cols() == 3);
591 assert(V.rows() == S.size());
592 GEO::Mesh M, M_out;
593 to_geogram_mesh(V, F, M);
594
595 // Remeshing options
596 opt.metric_attribute = "scalar";
597
598 GEO::Attribute<double> scalar(M.vertices.attributes(), opt.metric_attribute);
599 for (int v = 0; v < M.vertices.nb(); ++v)
600 {
601 scalar[v] = S(v);
602 }
603
604 // Remesh surface
605 mmg2d_tri_remesh(M, M_out, opt);
606
607 // Convert output
608 Eigen::MatrixXi OT;
609 from_geogram_mesh(M_out, OV, OF, OT);
610
611 // Promised a 2D mesh, so drop z-coordinate
612 OV.conservativeResize(OV.rows(), 2);
613 }
614
615 void remesh_adaptive_3d(const Eigen::MatrixXd &V, const Eigen::MatrixXi &T, const Eigen::VectorXd &S,
616 Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::MatrixXi &OT, MmgOptions opt)
617 {
618 assert(V.cols() == 3);
619 assert(V.rows() == S.size());
620 GEO::Mesh M, M_out;
621 to_geogram_mesh(V, Eigen::MatrixXi(0, 3), T, M);
622
623 // Remeshing options
624 opt.metric_attribute = "scalar";
625
626 GEO::Attribute<double> scalar(M.vertices.attributes(), opt.metric_attribute);
627 for (int v = 0; v < M.vertices.nb(); ++v)
628 {
629 scalar[v] = S(v);
630 }
631
632 // Remesh volume
633 mmg3d_tet_remesh(M, M_out, opt);
634
635 // Convert output
636 from_geogram_mesh(M_out, OV, OF, OT);
637 }
638
639} // namespace polyfem::mesh
640
641#endif
V
int V
Definition AdjointNLProblem.cpp:46
eigs.M
M
Definition eigs.py:94
generate_rotation_mtx.e
list e
Definition generate_rotation_mtx.py:49
generate_rotation_mtx.T
T
Definition generate_rotation_mtx.py:43
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:575
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:639
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:42