Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MeshUtils.cpp
Go to the documentation of this file.
1
2#include "MeshUtils.hpp"
3
10
11#include <unordered_set>
12
13#include <igl/PI.h>
14#include <igl/read_triangle_mesh.h>
15
16#include <geogram/mesh/mesh_io.h>
17#include <geogram/mesh/mesh_geometry.h>
18#include <geogram/basic/geometry.h>
19#include <geogram/mesh/mesh_preprocessing.h>
20#include <geogram/mesh/mesh_topology.h>
21#include <geogram/mesh/mesh_geometry.h>
22#include <geogram/mesh/mesh_repair.h>
23#include <geogram/mesh/mesh_AABB.h>
24#include <geogram/voronoi/CVT.h>
25#include <geogram/basic/logger.h>
27
28using namespace polyfem::io;
29using namespace polyfem::utils;
30
31bool polyfem::mesh::is_planar(const GEO::Mesh &M, const double tol)
32{
33 if (M.vertices.dimension() == 2)
34 return true;
35
36 assert(M.vertices.dimension() == 3);
37 GEO::vec3 min_corner, max_corner;
38 GEO::get_bbox(M, &min_corner[0], &max_corner[0]);
39 const double diff = (max_corner[2] - min_corner[2]);
40
41 return fabs(diff) < tol;
42}
43
44GEO::vec3 polyfem::mesh::mesh_vertex(const GEO::Mesh &M, GEO::index_t v)
45{
46 using GEO::index_t;
47 GEO::vec3 p(0, 0, 0);
48 for (index_t d = 0; d < std::min(3u, (index_t)M.vertices.dimension()); ++d)
49 {
50 if (M.vertices.double_precision())
51 {
52 p[d] = M.vertices.point_ptr(v)[d];
53 }
54 else
55 {
56 p[d] = M.vertices.single_precision_point_ptr(v)[d];
57 }
58 }
59 return p;
60}
61
62// -----------------------------------------------------------------------------
63
64GEO::vec3 polyfem::mesh::facet_barycenter(const GEO::Mesh &M, GEO::index_t f)
65{
66 using GEO::index_t;
67 GEO::vec3 p(0, 0, 0);
68 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
69 {
70 p += polyfem::mesh::mesh_vertex(M, M.facets.vertex(f, lv));
71 }
72 return p / M.facets.nb_vertices(f);
73}
74
75// -----------------------------------------------------------------------------
76
77GEO::index_t polyfem::mesh::mesh_create_vertex(GEO::Mesh &M, const GEO::vec3 &p)
78{
79 using GEO::index_t;
80 auto v = M.vertices.create_vertex();
81 for (index_t d = 0; d < std::min(3u, (index_t)M.vertices.dimension()); ++d)
82 {
83 if (M.vertices.double_precision())
84 {
85 M.vertices.point_ptr(v)[d] = p[d];
86 }
87 else
88 {
89 M.vertices.single_precision_point_ptr(v)[d] = (float)p[d];
90 }
91 }
92 return v;
93}
94
96
97void polyfem::mesh::compute_element_tags(const GEO::Mesh &M, std::vector<ElementType> &element_tags)
98{
99 using GEO::index_t;
100
101 std::vector<ElementType> old_tags = element_tags;
102
103 element_tags.resize(M.facets.nb());
104
105 // Step 0: Compute boundary vertices as true boundary + vertices incident to a polygon
106 GEO::Attribute<bool> is_boundary_vertex(M.vertices.attributes(), "boundary_vertex");
107 std::vector<bool> is_interface_vertex(M.vertices.nb(), false);
108 {
109 for (index_t f = 0; f < M.facets.nb(); ++f)
110 {
111 if (M.facets.nb_vertices(f) != 4 || (!old_tags.empty() && old_tags[f] == ElementType::INTERIOR_POLYTOPE))
112 {
113 // Vertices incident to polygonal facets (triangles or > 4 vertices) are marked as interface
114 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
115 {
116 is_interface_vertex[M.facets.vertex(f, lv)] = true;
117 }
118 }
119 }
120 }
121
122 // Step 1: Determine which vertices are regular or not
123 //
124 // Interior vertices are regular if they are incident to exactly 4 quads
125 // Boundary vertices are regular if they are incident to at most 2 quads, and no other facets
126 std::vector<int> degree(M.vertices.nb(), 0);
127 std::vector<bool> is_regular_vertex(M.vertices.nb());
128 for (index_t f = 0; f < M.facets.nb(); ++f)
129 {
130 if (M.facets.nb_vertices(f) == 4)
131 {
132 // Only count incident quads for the degree
133 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
134 {
135 index_t v = M.facets.vertex(f, lv);
136 degree[v]++;
137 }
138 }
139 }
140 for (index_t v = 0; v < M.vertices.nb(); ++v)
141 {
142 // assert(degree[v] > 0); // We assume there are no isolated vertices here
143 if (is_boundary_vertex[v] || is_interface_vertex[v])
144 {
145 is_regular_vertex[v] = (degree[v] <= 2);
146 }
147 else
148 {
149 is_regular_vertex[v] = (degree[v] == 4);
150 }
151 }
152
153 // Step 2: Iterate over the facets and determine the type
154 for (index_t f = 0; f < M.facets.nb(); ++f)
155 {
156 assert(M.facets.nb_vertices(f) > 2);
157 if (!old_tags.empty() && old_tags[f] == ElementType::INTERIOR_POLYTOPE)
158 continue;
159
160 if (M.facets.nb_vertices(f) == 4)
161 {
162 // Quad facet
163
164 // a) Determine if it is on the mesh boundary
165 bool is_boundary_facet = false;
166 bool is_interface_facet = false;
167 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
168 {
169 if (is_boundary_vertex[M.facets.vertex(f, lv)])
170 {
171 is_boundary_facet = true;
172 }
173 if (is_interface_vertex[M.facets.vertex(f, lv)])
174 {
175 is_interface_facet = true;
176 }
177 }
178
179 // b) Determine if it is regular or not
180 if (is_boundary_facet || is_interface_facet)
181 {
182 // A boundary quad is regular iff all its vertices are incident to at most 2 other quads
183 // We assume that non-boundary vertices of a boundary quads are always regular
184 bool is_singular = false;
185 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
186 {
187 index_t v = M.facets.vertex(f, lv);
188 if (is_boundary_vertex[v] || is_interface_vertex[v])
189 {
190 if (!is_regular_vertex[v])
191 {
192 is_singular = true;
193 break;
194 }
195 }
196 else
197 {
198 if (!is_regular_vertex[v])
199 {
200 element_tags[f] = ElementType::UNDEFINED;
201 break;
202 }
203 }
204 }
205
206 if (is_interface_facet)
207 {
208 element_tags[f] = ElementType::INTERFACE_CUBE;
209 }
210 else if (is_singular)
211 {
212 element_tags[f] = ElementType::SIMPLE_SINGULAR_BOUNDARY_CUBE;
213 }
214 else
215 {
216 element_tags[f] = ElementType::REGULAR_BOUNDARY_CUBE;
217 }
218 }
219 else
220 {
221 // An interior quad is regular if all its vertices are singular
222 int nb_singulars = 0;
223 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
224 {
225 if (!is_regular_vertex[M.facets.vertex(f, lv)])
226 {
227 ++nb_singulars;
228 }
229 }
230
231 if (nb_singulars == 0)
232 {
233 element_tags[f] = ElementType::REGULAR_INTERIOR_CUBE;
234 }
235 else if (nb_singulars == 1)
236 {
237 element_tags[f] = ElementType::SIMPLE_SINGULAR_INTERIOR_CUBE;
238 }
239 else
240 {
241 element_tags[f] = ElementType::MULTI_SINGULAR_INTERIOR_CUBE;
242 }
243 }
244 }
245 else
246 {
247 // Polygonal facet
248
249 // Note: In this function, we consider triangles as polygonal facets
250 ElementType tag = ElementType::INTERIOR_POLYTOPE;
251 GEO::Attribute<bool> boundary_vertices(M.vertices.attributes(), "boundary_vertex");
252 for (index_t lv = 0; lv < M.facets.nb_vertices(f); ++lv)
253 {
254 if (boundary_vertices[M.facets.vertex(f, lv)])
255 {
256 tag = ElementType::BOUNDARY_POLYTOPE;
257 break;
258 }
259 }
260
261 element_tags[f] = tag;
262 }
263 }
264
265 // TODO what happens at the neighs?
266 // Override for simplices
267 for (index_t f = 0; f < M.facets.nb(); ++f)
268 {
269 if (M.facets.nb_vertices(f) == 3)
270 {
271 element_tags[f] = ElementType::SIMPLEX;
272 }
273 }
274}
275
277
278namespace
279{
280
281 // Signed area of a polygonal facet
282 double signed_area(const GEO::Mesh &M, GEO::index_t f)
283 {
284 using namespace GEO;
285 double result = 0;
286 index_t v0 = M.facet_corners.vertex(M.facets.corners_begin(f));
287 const vec3 &p0 = Geom::mesh_vertex(M, v0);
288 for (index_t c =
289 M.facets.corners_begin(f) + 1;
290 c + 1 < M.facets.corners_end(f); c++)
291 {
292 index_t v1 = M.facet_corners.vertex(c);
293 const vec3 &p1 = polyfem::mesh::mesh_vertex(M, v1);
294 index_t v2 = M.facet_corners.vertex(c + 1);
295 const vec3 &p2 = polyfem::mesh::mesh_vertex(M, v2);
296 result += Geom::triangle_signed_area(vec2(&p0[0]), vec2(&p1[0]), vec2(&p2[0]));
297 }
298 return result;
299 }
300
301} // anonymous namespace
302
304{
305 using namespace GEO;
306 vector<index_t> component;
307 index_t nb_components = get_connected_components(M, component);
308 vector<double> comp_signed_volume(nb_components, 0.0);
309 for (index_t f = 0; f < M.facets.nb(); ++f)
310 {
311 comp_signed_volume[component[f]] += signed_area(M, f);
312 }
313 for (index_t f = 0; f < M.facets.nb(); ++f)
314 {
315 if (comp_signed_volume[component[f]] < 0.0)
316 {
317 M.facets.flip(f);
318 }
319 }
320}
321
323
324void polyfem::mesh::reorder_mesh(Eigen::MatrixXd &V, Eigen::MatrixXi &F, const Eigen::VectorXi &C, Eigen::VectorXi &R)
325{
326 assert(V.rows() == C.size());
327 int num_colors = C.maxCoeff() + 1;
328 Eigen::VectorXi count(num_colors);
329 count.setZero();
330 for (int i = 0; i < C.size(); ++i)
331 {
332 ++count[C(i)];
333 }
334 R.resize(num_colors + 1);
335 R(0) = 0;
336 for (int c = 0; c < num_colors; ++c)
337 {
338 R(c + 1) = R(c) + count(c);
339 }
340 count.setZero();
341 Eigen::VectorXi remap(C.size());
342 for (int i = 0; i < C.size(); ++i)
343 {
344 remap[i] = R(C(i)) + count[C(i)];
345 ++count[C(i)];
346 }
347 // Remap vertices
348 Eigen::MatrixXd NV(V.rows(), V.cols());
349 for (int v = 0; v < V.rows(); ++v)
350 {
351 NV.row(remap(v)) = V.row(v);
352 }
353 V = NV;
354 // Remap face indices
355 for (int f = 0; f < F.rows(); ++f)
356 {
357 for (int lv = 0; lv < F.cols(); ++lv)
358 {
359 F(f, lv) = remap(F(f, lv));
360 }
361 }
362}
363
365
366namespace
367{
368
369 void compute_unsigned_distance_field(const GEO::Mesh &M,
370 const GEO::MeshFacetsAABB &aabb_tree, const Eigen::MatrixXd &P, Eigen::VectorXd &D)
371 {
372 assert(P.cols() == 3);
373 D.resize(P.rows());
374#pragma omp parallel for
375 for (int i = 0; i < P.rows(); ++i)
376 {
377 GEO::vec3 pos(P(i, 0), P(i, 1), P(i, 2));
378 double sq_dist = aabb_tree.squared_distance(pos);
379 D(i) = sq_dist;
380 }
381 }
382
383 // calculate twice signed area of triangle (0,0)-(x1,y1)-(x2,y2)
384 // return an SOS-determined sign (-1, +1, or 0 only if it's a truly degenerate triangle)
385 int orientation(
386 double x1, double y1, double x2, double y2, double &twice_signed_area)
387 {
388 twice_signed_area = y1 * x2 - x1 * y2;
389 if (twice_signed_area > 0)
390 return 1;
391 else if (twice_signed_area < 0)
392 return -1;
393 else if (y2 > y1)
394 return 1;
395 else if (y2 < y1)
396 return -1;
397 else if (x1 > x2)
398 return 1;
399 else if (x1 < x2)
400 return -1;
401 else
402 return 0; // only true when x1==x2 and y1==y2
403 }
404
405 // robust test of (x0,y0) in the triangle (x1,y1)-(x2,y2)-(x3,y3)
406 // if true is returned, the barycentric coordinates are set in a,b,c.
407 //
408 // Note: This function comes from SDFGen by Christopher Batty.
409 // https://github.com/christopherbatty/SDFGen/blob/master/makelevelset3.cpp
410 bool point_in_triangle_2d(
411 double x0, double y0, double x1, double y1,
412 double x2, double y2, double x3, double y3,
413 double &a, double &b, double &c)
414 {
415 x1 -= x0;
416 x2 -= x0;
417 x3 -= x0;
418 y1 -= y0;
419 y2 -= y0;
420 y3 -= y0;
421 int signa = orientation(x2, y2, x3, y3, a);
422 if (signa == 0)
423 return false;
424 int signb = orientation(x3, y3, x1, y1, b);
425 if (signb != signa)
426 return false;
427 int signc = orientation(x1, y1, x2, y2, c);
428 if (signc != signa)
429 return false;
430 double sum = a + b + c;
431 geo_assert(sum != 0); // if the SOS signs match and are nonzero, there's no way all of a, b, and c are zero.
432 a /= sum;
433 b /= sum;
434 c /= sum;
435 return true;
436 }
437
438 // -----------------------------------------------------------------------------
439
440 // \brief Computes the (approximate) orientation predicate in 2d.
441 // \details Computes the sign of the (approximate) signed volume of
442 // the triangle p0, p1, p2
443 // \param[in] p0 first vertex of the triangle
444 // \param[in] p1 second vertex of the triangle
445 // \param[in] p2 third vertex of the triangle
446 // \retval POSITIVE if the triangle is oriented positively
447 // \retval ZERO if the triangle is flat
448 // \retval NEGATIVE if the triangle is oriented negatively
449 // \todo check whether orientation is inverted as compared to
450 // Shewchuk's version.
451 // Taken from geogram/src/lib/geogram/delaunay/delaunay_2d.cpp
452 inline GEO::Sign orient_2d_inexact(GEO::vec2 p0, GEO::vec2 p1, GEO::vec2 p2)
453 {
454 double a11 = p1[0] - p0[0];
455 double a12 = p1[1] - p0[1];
456
457 double a21 = p2[0] - p0[0];
458 double a22 = p2[1] - p0[1];
459
460 double Delta = GEO::det2x2(
461 a11, a12,
462 a21, a22);
463
464 return GEO::geo_sgn(Delta);
465 }
466
467 // -----------------------------------------------------------------------------
468
479 template <int X = 0, int Y = 1, int Z = 2>
480 int intersect_ray_z(const GEO::Mesh &M, GEO::index_t f, const GEO::vec3 &q, double &z)
481 {
482 using namespace GEO;
483
484 index_t c = M.facets.corners_begin(f);
485 const vec3 &p1 = Geom::mesh_vertex(M, M.facet_corners.vertex(c++));
486 const vec3 &p2 = Geom::mesh_vertex(M, M.facet_corners.vertex(c++));
487 const vec3 &p3 = Geom::mesh_vertex(M, M.facet_corners.vertex(c));
488
489 double u, v, w;
490 if (point_in_triangle_2d(
491 q[X], q[Y], p1[X], p1[Y], p2[X], p2[Y], p3[X], p3[Y], u, v, w))
492 {
493 z = u * p1[Z] + v * p2[Z] + w * p3[Z];
494 auto sign = orient_2d_inexact(vec2(p1[X], p1[Y]), vec2(p2[X], p2[Y]), vec2(p3[X], p3[Y]));
495 switch (sign)
496 {
497 case GEO::POSITIVE:
498 return 1;
499 case GEO::NEGATIVE:
500 return -1;
501 case GEO::ZERO:
502 default:
503 return 0;
504 }
505 }
506
507 return 0;
508 }
509
510 // -----------------------------------------------------------------------------
511
512 void compute_sign(const GEO::Mesh &M, const GEO::MeshFacetsAABB &aabb_tree,
513 const Eigen::MatrixXd &P, Eigen::VectorXd &D)
514 {
515 assert(P.cols() == 3);
516 assert(D.size() == P.rows());
517
518 GEO::vec3 min_corner, max_corner;
519 GEO::get_bbox(M, &min_corner[0], &max_corner[0]);
520
521#pragma omp parallel for
522 for (int k = 0; k < P.rows(); ++k)
523 {
524 GEO::vec3 center(P(k, 0), P(k, 1), P(k, 2));
525
526 GEO::Box box;
527 box.xyz_min[0] = box.xyz_max[0] = center[0];
528 box.xyz_min[1] = box.xyz_max[1] = center[1];
529 box.xyz_min[2] = min_corner[2];
530 box.xyz_max[2] = max_corner[2];
531
532 std::vector<std::pair<double, int>> inter;
533 auto action = [&M, &inter, &center](GEO::index_t f) {
534 double z;
535 if (int s = intersect_ray_z(M, f, center, z))
536 {
537 inter.emplace_back(z, s);
538 }
539 };
540 aabb_tree.compute_bbox_facet_bbox_intersections(box, action);
541 std::sort(inter.begin(), inter.end());
542
543 std::vector<double> reduced;
544 for (int i = 0, s = 0; i < (int)inter.size(); ++i)
545 {
546 const int ds = inter[i].second;
547 s += ds;
548 if ((s == -1 && ds < 0) || (s == 0 && ds > 0))
549 {
550 reduced.push_back(inter[i].first);
551 }
552 }
553
554 int num_before = 0;
555 for (double z : reduced)
556 {
557 if (z < center[2])
558 {
559 ++num_before;
560 }
561 }
562 if (num_before % 2 == 1)
563 {
564 // Point is inside
565 D(k) *= -1.0;
566 }
567 }
568 }
569
570} // anonymous namespace
571
573
574void polyfem::mesh::to_geogram_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, GEO::Mesh &M)
575{
576 M.clear();
577 // Setup vertices
578 M.vertices.create_vertices((int)V.rows());
579 for (int i = 0; i < (int)M.vertices.nb(); ++i)
580 {
581 GEO::vec3 &p = M.vertices.point(i);
582 p[0] = V(i, 0);
583 p[1] = V(i, 1);
584 p[2] = V.cols() >= 3 ? V(i, 2) : 0;
585 }
586 // Setup faces
587 if (F.cols() == 3)
588 {
589 M.facets.create_triangles((int)F.rows());
590 }
591 else if (F.cols() == 4)
592 {
593 M.facets.create_quads((int)F.rows());
594 }
595 else
596 {
597 throw std::runtime_error("Mesh format not supported");
598 }
599 for (int c = 0; c < (int)M.facets.nb(); ++c)
600 {
601 for (int lv = 0; lv < F.cols(); ++lv)
602 {
603 M.facets.set_vertex(c, lv, F(c, lv));
604 }
605 }
606}
607
608// void polyfem::mesh::to_geogram_mesh_3d(const Eigen::MatrixXd &V, const Eigen::MatrixXi &C, GEO::Mesh &M) {
609// M.clear();
610// // Setup vertices
611// M.vertices.create_vertices((int) V.rows());
612// assert(V.cols() == 3);
613// for (int i = 0; i < (int) M.vertices.nb(); ++i) {
614// GEO::vec3 &p = M.vertices.point(i);
615// p[0] = V(i, 0);
616// p[1] = V(i, 1);
617// p[2] = V(i, 2);
618// }
619
620// if(C.cols() == 4)
621// M.cells.create_tets((int) C.rows());
622// else if(C.cols() == 8)
623// M.cells.create_hexes((int) C.rows());
624// else
625// assert(false);
626
627// for (int c = 0; c < (int) M.cells.nb(); ++c) {
628// for (int lv = 0; lv < C.cols(); ++lv) {
629// M.cells.set_vertex(c, lv, C(c, lv));
630// }
631// }
632// M.cells.connect();
633// // GEO::mesh_reorient(M);
634// }
635
636// -----------------------------------------------------------------------------
637
638void polyfem::mesh::from_geogram_mesh(const GEO::Mesh &M, Eigen::MatrixXd &V, Eigen::MatrixXi &F, Eigen::MatrixXi &T)
639{
640 V.resize(M.vertices.nb(), 3);
641 for (int i = 0; i < (int)M.vertices.nb(); ++i)
642 {
643 GEO::vec3 p = M.vertices.point(i);
644 V.row(i) << p[0], p[1], p[2];
645 }
646 assert(M.facets.are_simplices());
647 F.resize(M.facets.nb(), 3);
648 for (int c = 0; c < (int)M.facets.nb(); ++c)
649 {
650 for (int lv = 0; lv < 3; ++lv)
651 {
652 F(c, lv) = M.facets.vertex(c, lv);
653 }
654 }
655 assert(M.cells.are_simplices());
656 T.resize(M.cells.nb(), 4);
657 for (int c = 0; c < (int)M.cells.nb(); ++c)
658 {
659 for (int lv = 0; lv < 4; ++lv)
660 {
661 T(c, lv) = M.cells.vertex(c, lv);
662 }
663 }
664}
665
666// -----------------------------------------------------------------------------
667
668void polyfem::mesh::signed_squared_distances(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
669 const Eigen::MatrixXd &P, Eigen::VectorXd &D)
670{
671 GEO::Mesh M;
672 to_geogram_mesh(V, F, M);
673 GEO::MeshFacetsAABB aabb_tree(M);
674 compute_unsigned_distance_field(M, aabb_tree, P, D);
675 compute_sign(M, aabb_tree, P, D);
676}
677
678// -----------------------------------------------------------------------------
679
680double polyfem::mesh::signed_volume(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F)
681{
682 assert(F.cols() == 3);
683 assert(V.cols() == 3);
684 std::array<Eigen::RowVector3d, 4> t;
685 t[3] = Eigen::RowVector3d::Zero(V.cols());
686 double volume_total = 0;
687 for (int f = 0; f < F.rows(); ++f)
688 {
689 for (int lv = 0; lv < F.cols(); ++lv)
690 {
691 t[lv] = V.row(F(f, lv));
692 }
693 double vol = GEO::Geom::tetra_signed_volume(t[0].data(), t[1].data(), t[2].data(), t[3].data());
694 volume_total += vol;
695 }
696 return -volume_total;
697}
698
699// -----------------------------------------------------------------------------
700
701void polyfem::mesh::orient_closed_surface(const Eigen::MatrixXd &V, Eigen::MatrixXi &F, bool positive)
702{
703 if ((positive ? 1 : -1) * signed_volume(V, F) < 0)
704 {
705 for (int f = 0; f < F.rows(); ++f)
706 {
707 F.row(f) = F.row(f).reverse().eval();
708 }
709 }
710}
711
712// -----------------------------------------------------------------------------
713
714void polyfem::mesh::extract_polyhedra(const Mesh3D &mesh, std::vector<std::unique_ptr<GEO::Mesh>> &polys, bool triangulated)
715{
716 std::vector<int> vertex_g2l(mesh.n_vertices() + mesh.n_faces(), -1);
717 std::vector<int> vertex_l2g;
718 for (int c = 0; c < mesh.n_cells(); ++c)
719 {
720 if (!mesh.is_polytope(c))
721 {
722 continue;
723 }
724 auto poly = std::make_unique<GEO::Mesh>();
725 int nv = mesh.n_cell_vertices(c);
726 int nf = mesh.n_cell_faces(c);
727 poly->vertices.create_vertices((triangulated ? nv + nf : nv) + 1);
728 vertex_l2g.clear();
729 vertex_l2g.reserve(nv);
730 for (int lf = 0; lf < nf; ++lf)
731 {
732 GEO::vector<GEO::index_t> facet_vertices;
733 auto index = mesh.get_index_from_element(c, lf, 0);
734 for (int lv = 0; lv < mesh.n_face_vertices(index.face); ++lv)
735 {
736 Eigen::RowVector3d p = mesh.point(index.vertex);
737 if (vertex_g2l[index.vertex] < 0)
738 {
739 vertex_g2l[index.vertex] = vertex_l2g.size();
740 vertex_l2g.push_back(index.vertex);
741 }
742 int v1 = vertex_g2l[index.vertex];
743 facet_vertices.push_back(v1);
744 poly->vertices.point(v1) = GEO::vec3(p.data());
745 index = mesh.next_around_face(index);
746 }
747 if (triangulated)
748 {
749 GEO::vec3 p(0, 0, 0);
750 for (GEO::index_t lv = 0; lv < facet_vertices.size(); ++lv)
751 {
752 p += poly->vertices.point(facet_vertices[lv]);
753 }
754 p /= facet_vertices.size();
755 int v0 = vertex_l2g.size();
756 vertex_l2g.push_back(0);
757 poly->vertices.point(v0) = p;
758 for (GEO::index_t lv = 0; lv < facet_vertices.size(); ++lv)
759 {
760 int v1 = facet_vertices[lv];
761 int v2 = facet_vertices[(lv + 1) % facet_vertices.size()];
762 poly->facets.create_triangle(v0, v1, v2);
763 }
764 }
765 else
766 {
767 poly->facets.create_polygon(facet_vertices);
768 }
769 }
770 {
771 Eigen::RowVector3d p = mesh.kernel(c);
772 poly->vertices.point(nv) = GEO::vec3(p.data());
773 }
774 assert(vertex_l2g.size() == size_t(triangulated ? nv + nf : nv));
775
776 for (int v : vertex_l2g)
777 {
778 vertex_g2l[v] = -1;
779 }
780
781 poly->facets.compute_borders();
782 poly->facets.connect();
783
784 polys.emplace_back(std::move(poly));
785 }
786}
787
788// -----------------------------------------------------------------------------
789
790// In Geogram, local vertices of a hex are numbered as follows:
791//
792// v5────v7
793// ╱┆ ╱│
794// v1─┼──v3 │
795// │v4┄┄┄┼v6
796// │╱ │╱
797// v0────v2
798//
799// However, `get_ordered_vertices_from_hex()` retrieves the local vertices in
800// this order:
801//
802// v7────v6
803// ╱┆ ╱│
804// v4─┼──v5 │
805// │v3┄┄┄┼v2
806// │╱ │╱
807// v0────v1
808//
809
810void polyfem::mesh::to_geogram_mesh(const Mesh3D &mesh, GEO::Mesh &M)
811{
812 M.clear();
813 // Convert vertices
814 M.vertices.create_vertices((int)mesh.n_vertices());
815 for (int i = 0; i < (int)M.vertices.nb(); ++i)
816 {
817 auto pt = mesh.point(i);
818 GEO::vec3 &p = M.vertices.point(i);
819 p[0] = pt[0];
820 p[1] = pt[1];
821 p[2] = pt[2];
822 }
823 // Convert faces
824 for (int f = 0, lf = 0; f < mesh.n_faces(); ++f)
825 {
826 if (mesh.is_boundary_face(f))
827 {
828 int nv = mesh.n_face_vertices(f);
829 M.facets.create_polygon(nv);
830 for (int lv = 0; lv < nv; ++lv)
831 {
832 M.facets.set_vertex(lf, lv, mesh.face_vertex(f, lv));
833 }
834 ++lf;
835 }
836 }
837 // Convert cells
838 typedef std::array<int, 8> Vector8i;
839 Vector8i g2p = {{0, 4, 1, 5, 3, 7, 2, 6}};
840 for (int c = 0; c < mesh.n_cells(); ++c)
841 {
842 if (mesh.is_cube(c))
843 {
844 Vector8i lvp = mesh.get_ordered_vertices_from_hex(c);
845 Vector8i lvg;
846 for (size_t k = 0; k < 8; ++k)
847 {
848 lvg[k] = lvp[g2p[k]];
849 }
850 std::reverse(lvg.begin(), lvg.end());
851 M.cells.create_hex(
852 lvg[0], lvg[1], lvg[2], lvg[3],
853 lvg[4], lvg[5], lvg[6], lvg[7]);
854 }
855 else
856 {
857 // TODO: Support conversion of tets as well!
858 }
859 }
860 M.facets.connect();
861 M.cells.connect();
862 // M.cells.compute_borders();
863 GEO::mesh_reorient(M);
864}
865
867
868void polyfem::mesh::tertrahedralize_star_shaped_surface(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
869 const Eigen::RowVector3d &kernel, Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::MatrixXi &OT)
870{
871 assert(V.cols() == 3);
872 OV.resize(V.rows() + 1, V.cols());
873 OV.topRows(V.rows()) = V;
874 OV.bottomRows(1) = kernel;
875 OF = F;
876 OT.resize(OF.rows(), 4);
877 OT.col(0).setConstant(V.rows());
878 OT.rightCols(3) = F;
879}
880
882
883void polyfem::mesh::sample_surface(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, int num_samples,
884 Eigen::MatrixXd &P, Eigen::MatrixXd *N, int num_lloyd, int num_newton)
885{
886 assert(num_samples > 3);
887 GEO::Mesh M;
888 to_geogram_mesh(V, F, M);
889 GEO::CentroidalVoronoiTesselation CVT(&M);
890 // GEO::mesh_save(M, "foo.obj");
891 bool was_quiet = GEO::Logger::instance()->is_quiet();
892 GEO::Logger::instance()->set_quiet(true);
893 CVT.compute_initial_sampling(num_samples);
894 GEO::Logger::instance()->set_quiet(was_quiet);
895
896 if (num_lloyd > 0)
897 {
898 CVT.Lloyd_iterations(num_lloyd);
899 }
900
901 if (num_newton > 0)
902 {
903 CVT.Newton_iterations(num_newton);
904 }
905
906 P.resize(3, num_samples);
907 std::copy_n(CVT.embedding(0), 3 * num_samples, P.data());
908 P.transposeInPlace();
909
910 if (N)
911 {
912 GEO::MeshFacetsAABB aabb(M);
913 N->resizeLike(P);
914 for (int i = 0; i < num_samples; ++i)
915 {
916 GEO::vec3 p(P(i, 0), P(i, 1), P(i, 2));
917 GEO::vec3 nearest_point;
918 double sq_dist;
919 auto f = aabb.nearest_facet(p, nearest_point, sq_dist);
920 GEO::vec3 n = normalize(GEO::Geom::mesh_facet_normal(M, f));
921 N->row(i) << n[0], n[1], n[2];
922 }
923 }
924}
925
927
928namespace
929{
930
931 bool approx_aligned(const double *a_, const double *b_, const double *p_, const double *q_, double tol = 1e-6)
932 {
933 using namespace GEO;
934 vec3 a(a_), b(b_), p(p_), q(q_);
935 double da = std::sqrt(Geom::point_segment_squared_distance(a, p, q));
936 double db = std::sqrt(Geom::point_segment_squared_distance(b, p, q));
937 double cos_theta = Geom::cos_angle(b - a, p - q);
938 return (da < tol && db < tol && std::abs(std::abs(cos_theta) - 1.0) < tol);
939 }
940
941} // anonymous namespace
942
943// -----------------------------------------------------------------------------
944
945void polyfem::mesh::extract_parent_edges(const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IE,
946 const Eigen::MatrixXd &BV, const Eigen::MatrixXi &BE, Eigen::MatrixXi &OE)
947{
948 assert(IV.cols() == 2 || IV.cols() == 3);
949 assert(BV.cols() == 2 || BV.cols() == 3);
950 typedef std::pair<int, int> Edge;
951 std::vector<Edge> selected;
952 for (int e1 = 0; e1 < IE.rows(); ++e1)
953 {
954 Eigen::RowVector3d a;
955 a.setZero();
956 a.head(IV.cols()) = IV.row(IE(e1, 0));
957 Eigen::RowVector3d b;
958 b.setZero();
959 b.head(IV.cols()) = IV.row(IE(e1, 1));
960 for (int e2 = 0; e2 < BE.rows(); ++e2)
961 {
962 Eigen::RowVector3d p;
963 p.setZero();
964 p.head(BV.cols()) = BV.row(BE(e2, 0));
965 Eigen::RowVector3d q;
966 q.setZero();
967 q.head(BV.cols()) = BV.row(BE(e2, 1));
968 if (approx_aligned(a.data(), b.data(), p.data(), q.data()))
969 {
970 selected.emplace_back(IE(e1, 0), IE(e1, 1));
971 break;
972 }
973 }
974 }
975
976 OE.resize(selected.size(), 2);
977 for (int e = 0; e < OE.rows(); ++e)
978 {
979 OE.row(e) << selected[e].first, selected[e].second;
980 }
981}
982
984
986 const Eigen::MatrixXd &vertices,
987 const Eigen::MatrixXi &codim_edges,
988 const Eigen::MatrixXi &faces,
989 Eigen::VectorXi &codim_vertices)
990{
991 std::vector<bool> is_vertex_codim(vertices.rows(), true);
992 for (int i = 0; i < codim_edges.rows(); i++)
993 {
994 for (int j = 0; j < codim_edges.cols(); j++)
995 {
996 is_vertex_codim[codim_edges(i, j)] = false;
997 }
998 }
999 for (int i = 0; i < faces.rows(); i++)
1000 {
1001 for (int j = 0; j < faces.cols(); j++)
1002 {
1003 is_vertex_codim[faces(i, j)] = false;
1004 }
1005 }
1006 const auto n_codim_vertices = std::count(is_vertex_codim.begin(), is_vertex_codim.end(), true);
1007 codim_vertices.resize(n_codim_vertices);
1008 for (int i = 0, ci = 0; i < vertices.rows(); i++)
1009 {
1010 if (is_vertex_codim[i])
1011 {
1012 codim_vertices[ci++] = i;
1013 }
1014 }
1015}
1016
1018 const Eigen::MatrixXi &tets,
1019 Eigen::MatrixXi &faces)
1020{
1021 std::unordered_set<Eigen::Vector3i, HashMatrix> tri_to_tet(4 * tets.rows());
1022 for (int i = 0; i < tets.rows(); i++)
1023 {
1024 tri_to_tet.emplace(tets(i, 0), tets(i, 2), tets(i, 1));
1025 tri_to_tet.emplace(tets(i, 0), tets(i, 3), tets(i, 2));
1026 tri_to_tet.emplace(tets(i, 0), tets(i, 1), tets(i, 3));
1027 tri_to_tet.emplace(tets(i, 1), tets(i, 2), tets(i, 3));
1028 }
1029
1030 std::vector<Eigen::RowVector3i> faces_vector;
1031 for (const auto &tri : tri_to_tet)
1032 {
1033 // find dual triangle with reversed indices:
1034 bool is_surface_triangle =
1035 tri_to_tet.find(Eigen::Vector3i(tri[2], tri[1], tri[0])) == tri_to_tet.end()
1036 && tri_to_tet.find(Eigen::Vector3i(tri[1], tri[0], tri[2])) == tri_to_tet.end()
1037 && tri_to_tet.find(Eigen::Vector3i(tri[0], tri[2], tri[1])) == tri_to_tet.end();
1038 if (is_surface_triangle)
1039 {
1040 faces_vector.emplace_back(tri[0], tri[1], tri[2]);
1041 }
1042 }
1043
1044 faces.resize(faces_vector.size(), 3);
1045 for (int i = 0; i < faces.rows(); i++)
1046 {
1047 faces.row(i) = faces_vector[i];
1048 }
1049}
1050
1052 const Eigen::MatrixXd &vertices,
1053 const Eigen::MatrixXi &tets,
1054 Eigen::MatrixXd &surface_vertices,
1055 Eigen::MatrixXi &tris)
1056{
1057 Eigen::MatrixXi full_tris;
1058 find_triangle_surface_from_tets(tets, full_tris);
1059
1060 std::unordered_map<int, int> full_to_surface;
1061 std::vector<size_t> surface_to_full;
1062 for (int i = 0; i < full_tris.rows(); i++)
1063 {
1064 for (int j = 0; j < full_tris.cols(); j++)
1065 {
1066 if (full_to_surface.find(full_tris(i, j)) == full_to_surface.end())
1067 {
1068 full_to_surface[full_tris(i, j)] = surface_to_full.size();
1069 surface_to_full.push_back(full_tris(i, j));
1070 }
1071 }
1072 }
1073
1074 surface_vertices.resize(surface_to_full.size(), 3);
1075 for (int i = 0; i < surface_to_full.size(); i++)
1076 {
1077 surface_vertices.row(i) = vertices.row(surface_to_full[i]);
1078 }
1079
1080 tris.resize(full_tris.rows(), full_tris.cols());
1081 for (int i = 0; i < tris.rows(); i++)
1082 {
1083 for (int j = 0; j < tris.cols(); j++)
1084 {
1085 tris(i, j) = full_to_surface[full_tris(i, j)];
1086 }
1087 }
1088}
1089
1091 const std::string &mesh_path,
1092 Eigen::MatrixXd &vertices,
1093 Eigen::VectorXi &codim_vertices,
1094 Eigen::MatrixXi &codim_edges,
1095 Eigen::MatrixXi &faces)
1096{
1097 vertices.resize(0, 0);
1098 codim_vertices.resize(0);
1099 codim_edges.resize(0, 0);
1100 faces.resize(0, 0);
1101
1102 std::string lowername = mesh_path;
1103 std::transform(
1104 lowername.begin(), lowername.end(), lowername.begin(), ::tolower);
1105
1106 if (StringUtils::endswith(lowername, ".msh"))
1107 {
1108 Eigen::MatrixXi cells;
1109 std::vector<std::vector<int>> elements;
1110 std::vector<std::vector<double>> weights;
1111 std::vector<int> body_ids;
1112 if (!MshReader::load(mesh_path, vertices, cells, elements, weights, body_ids))
1113 {
1114 logger().error("Unable to load mesh: {}", mesh_path);
1115 return false;
1116 }
1117
1118 if (cells.cols() == 1)
1119 codim_vertices = cells;
1120 else if (cells.cols() == 2)
1121 codim_edges = cells;
1122 else if (cells.cols() == 3)
1123 faces = cells;
1124 else if (cells.cols() == 4)
1125 {
1126 if (vertices.cols() == 2)
1127 {
1128 logger().error("read_surface_mesh not implemented for 2D quad meshes");
1129 return false;
1130 }
1131 else
1132 {
1133 // TODO: how to distinguish between 3D tet mesh and 3D surface quad mesh?
1134 assert(vertices.cols() == 3);
1135 Eigen::MatrixXd surface_vertices;
1136 extract_triangle_surface_from_tets(vertices, cells, surface_vertices, faces);
1137 vertices = surface_vertices;
1138 }
1139 }
1140 else
1141 {
1142 logger().error("read_surface_mesh not implemented for hexahedral and polygonal/polyhedral meshes");
1143 return false;
1144 }
1145 }
1146 else if (StringUtils::endswith(lowername, ".obj")) // Use specialized OBJ reader function with polyline support
1147 {
1148 if (!OBJReader::read(mesh_path, vertices, codim_edges, faces))
1149 {
1150 logger().error("Unable to load mesh: {}", mesh_path);
1151 return false;
1152 }
1153 }
1154 else if (!igl::read_triangle_mesh(mesh_path, vertices, faces))
1155 {
1156 GEO::Mesh mesh;
1157 if (!GEO::mesh_load(mesh_path, mesh))
1158 {
1159 logger().error("Unable to load mesh: {}", mesh_path);
1160 return false;
1161 }
1162
1163 int dim = is_planar(mesh) ? 2 : 3;
1164 vertices.resize(mesh.vertices.nb(), dim);
1165 for (int vi = 0; vi < mesh.vertices.nb(); vi++)
1166 {
1167 const auto &v = mesh.vertices.point(vi);
1168 for (int vj = 0; vj < dim; vj++)
1169 {
1170 vertices(vi, vj) = v[vj];
1171 }
1172 }
1173
1174 // TODO: Check that this works even for a volumetric mesh
1175 assert(mesh.facets.nb());
1176 int face_cols = mesh.facets.nb_vertices(0);
1177 faces.resize(mesh.facets.nb(), face_cols);
1178 for (int fi = 0; fi < mesh.facets.nb(); fi++)
1179 {
1180 assert(face_cols == mesh.facets.nb_vertices(fi));
1181 for (int fj = 0; fj < mesh.facets.nb_vertices(fi); fj++)
1182 {
1183 faces(fi, fj) = mesh.facets.vertex(fi, fj);
1184 }
1185 }
1186 }
1187
1188 find_codim_vertices(vertices, codim_edges, faces, codim_vertices);
1189
1190 return true;
1191}
1192
1193int polyfem::mesh::count_faces(const int dim, const Eigen::MatrixXi &cells)
1194{
1195 std::unordered_set<std::vector<int>, HashVector> boundaries;
1196
1197 auto insert = [&](std::vector<int> v) {
1198 std::sort(v.begin(), v.end());
1199 boundaries.insert(v);
1200 };
1201
1202 for (int i = 0; i < cells.rows(); i++)
1203 {
1204 const auto &cell = cells.row(i);
1205 if (cells.cols() == 3) // triangle
1206 {
1207 insert({{cell(0), cell(1)}});
1208 insert({{cell(1), cell(2)}});
1209 insert({{cell(2), cell(0)}});
1210 }
1211 else if (cells.cols() == 4 && dim == 2) // quadralateral
1212 {
1213 insert({{cell(0), cell(1)}});
1214 insert({{cell(1), cell(2)}});
1215 insert({{cell(2), cell(3)}});
1216 insert({{cell(3), cell(0)}});
1217 }
1218 else if (cells.cols() == 4 && dim == 3) // tetrahedron
1219 {
1220 insert({{cell(0), cell(2), cell(1)}});
1221 insert({{cell(0), cell(3), cell(2)}});
1222 insert({{cell(0), cell(1), cell(3)}});
1223 insert({{cell(1), cell(2), cell(3)}});
1224 }
1225 else if (cells.cols() == 8) // hexahedron
1226 {
1227 insert({{cell(0), cell(1), cell(2), cell(3)}});
1228 insert({{cell(1), cell(5), cell(6), cell(3)}});
1229 insert({{cell(5), cell(4), cell(7), cell(6)}});
1230 insert({{cell(0), cell(4), cell(7), cell(3)}});
1231 insert({{cell(0), cell(4), cell(5), cell(1)}});
1232 insert({{cell(2), cell(6), cell(7), cell(3)}});
1233 }
1234 else
1235 {
1236 throw std::runtime_error("count_boundary_elements not implemented for polygons");
1237 }
1238 }
1239
1240 return boundaries.size();
1241}
1242
1244{
1245 using namespace GEO;
1246 typedef std::pair<index_t, index_t> Edge;
1247
1248 if (M.edges.nb() > 0)
1249 return;
1250
1251 M.facets.connect();
1252 M.cells.connect();
1253 if (M.cells.nb() != 0 && M.facets.nb() == 0)
1254 {
1255 M.cells.compute_borders();
1256 }
1257
1258 // Compute a list of all the edges, and store edge index as a corner attribute
1259 std::vector<std::pair<Edge, index_t>> e2c; // edge to corner id
1260 for (index_t f = 0; f < M.facets.nb(); ++f)
1261 {
1262 for (index_t c = M.facets.corners_begin(f); c < M.facets.corners_end(f); ++c)
1263 {
1264 index_t v = M.facet_corners.vertex(c);
1265 index_t c2 = M.facets.next_corner_around_facet(f, c);
1266 index_t v2 = M.facet_corners.vertex(c2);
1267 e2c.emplace_back(std::make_pair(std::min(v, v2), std::max(v, v2)), c);
1268 }
1269 }
1270 std::sort(e2c.begin(), e2c.end());
1271
1272 // Assign unique id to edges
1273 M.edges.clear();
1274 Edge prev_e(-1, -1);
1275 for (const auto &kv : e2c)
1276 {
1277 Edge e = kv.first;
1278 if (e != prev_e)
1279 {
1280 M.edges.create_edge(e.first, e.second);
1281 prev_e = e;
1282 }
1283 }
1284}
int V
QuadratureVector da
Definition Assembler.cpp:23
void find_codim_vertices(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &codim_edges, const Eigen::MatrixXi &faces, Eigen::VectorXi &codim_vertices)
void find_triangle_surface_from_tets(const Eigen::MatrixXi &tets, Eigen::MatrixXi &faces)
std::vector< Eigen::VectorXi > faces
int z
static bool load(const std::string &path, Eigen::MatrixXd &vertices, Eigen::MatrixXi &cells, std::vector< std::vector< int > > &elements, std::vector< std::vector< double > > &weights, std::vector< int > &body_ids)
Definition MshReader.cpp:29
static bool read(const std::string obj_file_name, std::vector< std::vector< double > > &V, std::vector< std::vector< double > > &TC, std::vector< std::vector< double > > &N, std::vector< std::vector< int > > &F, std::vector< std::vector< int > > &FTC, std::vector< std::vector< int > > &FN, std::vector< std::vector< int > > &L)
Read a mesh from an ascii obj file.
Definition OBJReader.cpp:32
virtual Navigation3D::Index get_index_from_element(int hi, int lf, int lv) const =0
virtual RowVectorNd kernel(const int cell_id) const =0
std::array< int, 8 > get_ordered_vertices_from_hex(const int element_index) const
Definition Mesh3D.cpp:322
virtual int n_cell_faces(const int c_id) const =0
virtual Navigation3D::Index next_around_face(Navigation3D::Index idx) const =0
virtual int n_vertices() const =0
number of vertices
bool is_polytope(const int el_id) const
checks if element is polygon compatible
Definition Mesh.cpp:365
bool is_cube(const int el_id) const
checks if element is cube compatible
Definition Mesh.cpp:352
virtual RowVectorNd point(const int global_index) const =0
point coordinates
virtual bool is_boundary_face(const int face_global_id) const =0
is face boundary
virtual int n_cells() const =0
number of cells
virtual int n_faces() const =0
number of faces
virtual int n_face_vertices(const int f_id) const =0
number of vertices of a face
virtual int n_cell_vertices(const int c_id) const =0
number of vertices of a cell
virtual int face_vertex(const int f_id, const int lv_id) const =0
id of the face vertex
M
Definition eigs.py:94
Eigen::ArrayXd P(const int m, const int p, const Eigen::ArrayXd &z)
Definition p_n_bases.cpp:42
bool is_planar(const GEO::Mesh &M, const double tol=1e-5)
Determine if the given mesh is planar (2D or tiny z-range).
Definition MeshUtils.cpp:31
void sample_surface(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, int num_samples, Eigen::MatrixXd &P, Eigen::MatrixXd *N=nullptr, int num_lloyd=10, int num_newton=10)
Samples points on a surface.
void orient_closed_surface(const Eigen::MatrixXd &V, Eigen::MatrixXi &F, bool positive=true)
Orient a triangulated surface to have positive volume.
double signed_volume(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F)
Compute the signed volume of a surface mesh.
void reorder_mesh(Eigen::MatrixXd &V, Eigen::MatrixXi &F, const Eigen::VectorXi &C, Eigen::VectorXi &R)
Reorder vertices of a mesh using color tags, so that vertices are ordered by increasing colors.
void orient_normals_2d(GEO::Mesh &M)
Orient facets of a 2D mesh so that each connected component has positive volume.
void extract_parent_edges(const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IE, const Eigen::MatrixXd &BV, const Eigen::MatrixXi &BE, Eigen::MatrixXi &OE)
Extract a set of edges that are overlap with a set given set of parent edges, using vertices position...
ElementType
Type of Element, check [Poly-Spline Finite Element Method] for a complete description.
Definition Mesh.hpp:23
void extract_triangle_surface_from_tets(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &tets, Eigen::MatrixXd &surface_vertices, Eigen::MatrixXi &tris)
Extract triangular surface from a tetmesh.
GEO::vec3 mesh_vertex(const GEO::Mesh &M, GEO::index_t v)
Retrieve a 3D vector with the position of a given vertex.
Definition MeshUtils.cpp:44
GEO::index_t mesh_create_vertex(GEO::Mesh &M, const GEO::vec3 &p)
Definition MeshUtils.cpp:77
void to_geogram_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, GEO::Mesh &M)
Converts a triangle mesh to a Geogram mesh.
void signed_squared_distances(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXd &P, Eigen::VectorXd &D)
Computes the signed squared distance from a list of points to a triangle mesh.
GEO::vec3 facet_barycenter(const GEO::Mesh &M, GEO::index_t f)
Definition MeshUtils.cpp:64
void from_geogram_mesh(const GEO::Mesh &M, Eigen::MatrixXd &V, Eigen::MatrixXi &F, Eigen::MatrixXi &T)
Extract simplices from a Geogram mesh.
void tertrahedralize_star_shaped_surface(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::RowVector3d &kernel, Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::MatrixXi &OT)
Tetrahedralize a star-shaped mesh, with a given point in its kernel.
void generate_edges(GEO::Mesh &M)
assing edges to M
bool read_surface_mesh(const std::string &mesh_path, Eigen::MatrixXd &vertices, Eigen::VectorXi &codim_vertices, Eigen::MatrixXi &codim_edges, Eigen::MatrixXi &faces)
read a surface mesh
void extract_polyhedra(const Mesh3D &mesh, std::vector< std::unique_ptr< GEO::Mesh > > &polys, bool triangulated=false)
Extract polyhedra from a 3D volumetric mesh.
int count_faces(const int dim, const Eigen::MatrixXi &cells)
Count the number of boundary elements (triangles for tetmesh and edges for triangle mesh)
void compute_element_tags(const GEO::Mesh &M, std::vector< ElementType > &element_tags)
Compute the type of each facet in a surface mesh.
Definition MeshUtils.cpp:97
bool endswith(const std::string &str, const std::string &suffix)
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44