Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Refinement.cpp
Go to the documentation of this file.
1
2#include "Refinement.hpp"
3#include "Navigation.hpp"
5#include "PolygonUtils.hpp"
7
8#include <iostream>
9#include <vector>
10#include <array>
11#include <numeric>
13
14using namespace polyfem;
15
17 const Eigen::MatrixXi &Q, Eigen::MatrixXi &edge_index,
18 std::vector<std::vector<int>> &adj,
19 std::vector<std::pair<int, int>> *pairs_of_edges,
20 std::vector<std::pair<int, int>> *pairs_of_quads,
21 Eigen::MatrixXi *quad_index)
22{
23 assert(Q.cols() == 4); // Assumes quad mesh
24 typedef std::pair<int, int> Edge;
25
26 std::vector<std::pair<int, int>> edges;
27
28 edge_index.resizeLike(Q); // (f, lv) -> edge id
29 adj.clear();
30 edges.clear();
31 if (pairs_of_quads)
32 {
33 pairs_of_quads->clear();
34 }
35 if (quad_index)
36 {
37 quad_index->setConstant(Q.rows(), Q.cols(), -1);
38 }
39
40 // Number mesh edges & build dual edge-graph
41 std::vector<std::pair<Edge, std::pair<int, int>>> accu;
42 accu.reserve((size_t)Q.size());
43 for (int f = 0; f < Q.rows(); ++f)
44 {
45 for (int lv = 0; lv < 4; ++lv)
46 {
47 int x = Q(f, lv);
48 int y = Q(f, (lv + 1) % 4);
49 accu.emplace_back(Edge(std::min(x, y), std::max(x, y)), std::make_pair(f, lv));
50 }
51 }
52 std::sort(accu.begin(), accu.end());
53 edges.reserve(accu.size() / 2);
54 edge_index.setConstant(-1);
55 int id = -1;
56 int prev_quad = -1;
57 int prev_lv = -1;
58 Edge prev_side(-1, -1);
59 for (const auto &kv : accu)
60 {
61 // kv = (v1, v2), (f, lv)
62 Edge e = kv.first;
63 int f = kv.second.first;
64 int lv = kv.second.second;
65 // printf("(%d,%d) -> %d, %d\n", e.first, e.second, f, lv);
66 if (e != prev_side)
67 {
68 ++id;
69 edges.emplace_back(e);
70 }
71 else
72 {
73 if (pairs_of_quads)
74 {
75 if (quad_index)
76 {
77 (*quad_index)(f, lv) = prev_quad;
78 (*quad_index)(prev_quad, prev_lv) = f;
79 }
80 pairs_of_quads->emplace_back(prev_quad, f);
81 }
82 }
83 edge_index(f, lv) = id;
84 prev_side = e;
85 prev_quad = f;
86 prev_lv = lv;
87 }
88 adj.resize(edges.size());
89 for (int f = 0; f < Q.rows(); ++f)
90 {
91 for (int lv = 0; lv < 4; ++lv)
92 {
93 int e1 = edge_index(f, lv);
94 int e2 = edge_index(f, (lv + 2) % 4);
95 assert(e1 != -1);
96 if (e1 < e2)
97 {
98 adj[e1].push_back(e2);
99 adj[e2].push_back(e1);
100 }
101 }
102 }
103
104 if (pairs_of_edges)
105 {
106 std::swap(*pairs_of_edges, edges);
107 }
108}
109
111
112namespace
113{
114
115 enum class PatternType
116 {
117 NOT_PERIODIC, // Opposite borders do not match
118 NOT_SYMMETRIC, // At least one border is not symmetric along its bisector
119 SIMPLE_PERIODIC, // Opposite borders match individually
120 DOUBLE_PERIODIC, // All four borders match between each others
121 };
122
123 // -----------------------------------------------------------------------------
124
125 // Determine the type of periodicity of the given pattern. This function rely
126 // only on the input vertex positions, but does not check if the underlying
127 // topology is consistent.
128 //
129 // @param[in] V #V x 3 matrix of vertex positions
130 // @param[in] F #F x 3 matrix of triangle indexes
131 // @param[out] border_vertices List of vertices along each side
132 //
133 // @return { Type of periodicity of the pattern. }
134 //
135 PatternType compute_pattern_type(
136 const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
137 std::array<Eigen::VectorXi, 4> &border_vertices)
138 {
139 Eigen::Vector2d lower = V.colwise().minCoeff().head<2>();
140 Eigen::Vector2d upper = V.colwise().maxCoeff().head<2>();
141
142 std::array<Eigen::MatrixXd, 2> border;
143
144 // Compute vertices along the border. Sides are numbered as follows:
145 //
146 // 2
147 // v3 ----- v2
148 // | |
149 // 3 | | 1
150 // | |
151 // v0 ----- v1
152 // 0
153 // ↑y
154 // └─→x
155 for (int d = 0; d < 2; ++d)
156 {
157 std::vector<std::pair<double, int>> low, upp;
158 for (int i = 0; i < V.rows(); ++i)
159 {
160 if (V(i, 1 - d) == lower[1 - d])
161 {
162 low.emplace_back(V(i, d), i);
163 }
164 if (V(i, 1 - d) == upper[1 - d])
165 {
166 upp.emplace_back(V(i, d), i);
167 }
168 }
169 std::sort(low.begin(), low.end());
170 std::sort(upp.begin(), upp.end());
171 border_vertices[d ? 3 : 0].resize(low.size());
172 for (size_t i = 0; i < low.size(); ++i)
173 {
174 border_vertices[d ? 3 : 0][i] = low[i].second;
175 }
176 border_vertices[d ? 1 : 2].resize(upp.size());
177 for (size_t i = 0; i < upp.size(); ++i)
178 {
179 border_vertices[d ? 1 : 2][i] = upp[i].second;
180 }
181 }
182 for (int i : {2, 3})
183 {
184 border_vertices[i].reverseInPlace();
185 }
186
187 // Check if borders have the same topology and geometry
188 for (int d = 0; d < 2; ++d)
189 {
190 std::vector<double> low, upp;
191 for (int i = 0; i < V.rows(); ++i)
192 {
193 if (V(i, d) == lower[d])
194 {
195 low.push_back(V(i, 1 - d));
196 }
197 if (V(i, d) == upper[d])
198 {
199 upp.push_back(V(i, 1 - d));
200 }
201 }
202 std::sort(low.begin(), low.end());
203 std::sort(upp.begin(), upp.end());
204 if (low.size() != upp.size())
205 {
206 return PatternType::NOT_PERIODIC;
207 }
208 size_t n = low.size();
209 border[d].resize(n, 2);
210 for (size_t i = 0; i < n; ++i)
211 {
212 border[d](i, 0) = low[i];
213 border[d](i, 1) = upp[i];
214 }
215 if (!border[d].col(0).isApprox(border[d].col(1)))
216 {
217 return PatternType::NOT_PERIODIC;
218 }
219 }
220
221 // Check if borders are symmetric
222 for (int d = 0; d < 2; ++d)
223 {
224 if (!(border[d].col(0).array() - lower[d]).isApprox(upper[d] - border[d].col(0).array().reverse()))
225 {
226 return PatternType::NOT_SYMMETRIC;
227 }
228 }
229
230 // Check if horizontal and vertical borders are matching
231 if (border[0].size() == border[1].size())
232 {
233 if (!border[0].col(0).isApprox(border[1].col(0)))
234 {
235 logger().warn("Pattern boundaries have the same number of vertices, but their position differ slighly.");
236 Eigen::MatrixXd X(border[0].size(), 2);
237 X.col(0) = border[0].col(0);
238 X.col(1) = border[1].col(0);
239 }
240 return PatternType::DOUBLE_PERIODIC;
241 }
242 else
243 {
244 return PatternType::SIMPLE_PERIODIC;
245 }
246
247 return PatternType::SIMPLE_PERIODIC;
248 }
249
250 // -----------------------------------------------------------------------------
251
252 Eigen::VectorXi vertex_degree(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F)
253 {
254 Eigen::VectorXi count = Eigen::VectorXi::Zero(V.rows());
255
256 for (unsigned i = 0; i < F.rows(); ++i)
257 {
258 for (unsigned j = 0; j < F.cols(); ++j)
259 {
260 // Avoid duplicate edges
261 if (F(i, j) < F(i, (j + 1) % F.cols()))
262 {
263 count(F(i, j)) += 1;
264 count(F(i, (j + 1) % F.cols())) += 1;
265 }
266 }
267 }
268
269 return count;
270 }
271
272} // anonymous namespace
273
275
277 const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IQ,
278 const Eigen::MatrixXd &PV, const Eigen::MatrixXi &PF,
279 Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::VectorXi *SF,
280 EvalParametersFunc evalFunc, GetAdjacentLocalEdge getAdjLocalEdge)
281{
282 // List of vertices along each border (from lv to (lv+1)%4)
283 std::array<Eigen::VectorXi, 4> border_vertices;
284
285 auto pattern = compute_pattern_type(PV, PF, border_vertices);
286 auto valence = vertex_degree(IV, IQ);
287 // auto border = igl::is_border_vertex(IV, IQ);
288
289 // Check input validity
290 switch (pattern)
291 {
292 case PatternType::NOT_PERIODIC:
293 logger().error("Pattern not periodic");
294 return false;
295 case PatternType::NOT_SYMMETRIC:
296 logger().error("Pattern not symmetric");
297 return false;
298 case PatternType::SIMPLE_PERIODIC:
299 logger().error("Pattern simple periodic");
300 return false;
301 case PatternType::DOUBLE_PERIODIC:
302 // Not implemented
303 break;
304 default:
305 logger().error("Unknown patter type");
306 return false;
307 }
308
309 // Normalized coordinates (between 0 and 1 for barycentric coordinates)
310 auto lower = PV.colwise().minCoeff();
311 auto upper = PV.colwise().maxCoeff();
312 Eigen::MatrixXd PVN = (PV.rowwise() - lower).array().rowwise() / (upper - lower).cwiseMax(1e-5).array();
313
314 // If the eval function is undefined, don't do any remapping
315 if (!evalFunc)
316 {
317 evalFunc = [&](const Eigen::MatrixXd &uv, Eigen::MatrixXd &mapped, int q) {
318 const auto &u = uv.col(0).array();
319 const auto &v = uv.col(1).array();
320 Eigen::RowVector3d a = IV.row(IQ(q, 0));
321 Eigen::RowVector3d b = IV.row(IQ(q, 1));
322 Eigen::RowVector3d c = IV.row(IQ(q, 2));
323 Eigen::RowVector3d d = IV.row(IQ(q, 3));
324 mapped = ((1 - u) * (1 - v)).matrix() * a + (u * (1 - v)).matrix() * b + (u * v).matrix() * c + ((1 - u) * v).matrix() * d;
325 };
326 }
327
328 // Mapped uv samples
329 Eigen::MatrixXd mapped;
330
331 // Instantiate (duplicating vertices)
332 Eigen::MatrixXd P0 = PVN.row(0);
333 evalFunc(P0, mapped, 0);
334 Eigen::MatrixXd V(IQ.rows() * PV.rows(), mapped.cols());
335 Eigen::MatrixXi F(IQ.rows() * PF.rows(), PF.cols());
336 if (SF)
337 {
338 SF->resize(V.rows());
339 }
340
341 // Remapped vertex id (after duplicate removal)
342 Eigen::VectorXi remap(V.rows());
343 remap.setConstant(-1);
344
345 for (int q = 0; q < IQ.rows(); ++q)
346 {
347 evalFunc(PVN, mapped, q);
348 int v0 = (int)PVN.rows() * q;
349 V.middleRows(v0, mapped.rows()) = mapped;
350 int f0 = (int)PF.rows() * q;
351 F.middleRows(f0, PF.rows()) = PF.array() + v0;
352 if (SF)
353 {
354 SF->segment(v0, mapped.rows()).setConstant(q);
355 }
356 }
357
358 auto min_max = [](int x, int y) {
359 return std::make_pair(std::min(x, y), std::max(x, y));
360 };
361
362 if (getAdjLocalEdge)
363 {
364 // Adjacency function has already been provided
365 int cnt = 0;
366 for (int q1 = 0; q1 < IQ.rows(); ++q1)
367 {
368 const int v1 = (int)PVN.rows() * q1;
369 for (int lv1 = 0; lv1 < 4; ++lv1)
370 {
371 const auto res = getAdjLocalEdge(q1, lv1);
372 const int q2 = std::get<0>(res);
373 const int v2 = (int)PVN.rows() * q2;
374 const int lv2 = std::get<1>(res);
375 const bool rev = std::get<2>(res);
376 if (q2 > q1)
377 {
378 Eigen::VectorXi side1 = border_vertices[lv1].array() + v1;
379 Eigen::VectorXi side2 = border_vertices[lv2].array() + v2;
380 if (rev)
381 {
382 side2.reverseInPlace();
383 }
384 for (int ii = 0; ii < (int)side1.size(); ++ii)
385 {
386 const int x1 = side1[ii];
387 const int x2 = side2[ii];
388 if (remap(x1) < 0)
389 {
390 remap(x1) = cnt++;
391 }
392 remap(x2) = remap(x1);
393 }
394 }
395 }
396 }
397 for (int v = 0; v < V.rows(); ++v)
398 {
399 if (remap(v) < 0)
400 {
401 remap(v) = cnt++;
402 }
403 }
404 }
405 else
406 {
407 // Compute adjacency info on the quad mesh
408 Eigen::MatrixXi edge_index;
409 std::vector<std::vector<int>> adj;
410 std::vector<std::pair<int, int>> pairs_of_quads;
411 edge_adjacency_graph(IQ, edge_index, adj, nullptr, &pairs_of_quads);
412 // Stitch vertices from adjacent quads
413 for (const auto &qq : pairs_of_quads)
414 {
415 int q1 = qq.first;
416 int q2 = qq.second;
417 if (q2 > q1)
418 {
419 std::swap(q1, q2);
420 }
421 const int v1 = (int)PVN.rows() * q1;
422 const int v2 = (int)PVN.rows() * q2;
423 int lv1 = 0;
424 int lv2 = 0;
425 for (; lv1 < 4; ++lv1)
426 {
427 int x1 = IQ(q1, lv1);
428 int y1 = IQ(q1, (lv1 + 1) % 4);
429 for (lv2 = 0; lv2 < 4; ++lv2)
430 {
431 int x2 = IQ(q2, lv2);
432 int y2 = IQ(q2, (lv2 + 1) % 4);
433 if (min_max(x1, y1) == min_max(x2, y2))
434 {
435 int e = edge_index(q1, lv1);
436 // tfm::printf("quads: (%s, %s) and (%s, %s)\n", q1, lv1, q2, lv2);
437 // tfm::printf("└ edge: %s-%s (id: %e)\n", x1, y1, e);
438 assert(edge_index(q2, lv2) == e);
439 Eigen::VectorXi side1 = border_vertices[lv1].array() + v1;
440 Eigen::VectorXi side2 = border_vertices[lv2].array() + v2;
441 if (x1 > y1)
442 {
443 side1.reverseInPlace();
444 }
445 if (x2 > y2)
446 {
447 side2.reverseInPlace();
448 }
449 for (int ii = 0; ii < (int)side1.size(); ++ii)
450 {
451 remap(side2[ii]) = side1[ii];
452 }
453 lv1 = 4;
454 lv2 = 4;
455 }
456 }
457 }
458 }
459 }
460
461 // OV = V;
462 // OF = F;
463 // return true;
464
465 // Remap vertices according to 'remap'
466 int num_remapped = remap.maxCoeff() + 1;
467 OV.resize(num_remapped, V.cols());
468 for (int v = 0; v < V.rows(); ++v)
469 {
470 OV.row(remap(v)) = V.row(v);
471 }
472 for (int f = 0; f < F.rows(); ++f)
473 {
474 for (int lv = 0; lv < F.cols(); ++lv)
475 {
476 int ov = F(f, lv);
477 int nv = remap(ov);
478 F(f, lv) = nv;
479 }
480 }
481 OF = F;
482 // Eigen::VectorXi I;
483 // igl::remove_unreferenced(V, F, OV, OF, I);
484
485 // Remap tags on vertices
486 if (SF)
487 {
488 Eigen::VectorXi tmp = *SF;
489 SF->resize(OV.rows());
490 SF->setZero();
491 for (int v = 0; v < V.rows(); ++v)
492 {
493 (*SF)(remap(v)) = tmp(v);
494 }
495 }
496
497 return true;
498}
499
501
503 const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IF,
504 Eigen::MatrixXd &OV, Eigen::MatrixXi &OF)
505{
506 Eigen::MatrixXd PV(9, 3);
507 Eigen::MatrixXi PF(4, 4);
508 PV << -1, 0, 0,
509 0, 0, 0,
510 0, 1, 0,
511 -1, 1, 0,
512 1, 1, 0,
513 1, 0, 0,
514 1, -1, 0,
515 0, -1, 0,
516 -1, -1, 0;
517 PF << 1, 2, 3, 4,
518 3, 2, 6, 5,
519 6, 2, 8, 7,
520 2, 1, 9, 8;
521 PF = PF.array() - 1;
522 bool res = instantiate_pattern(IV, IF, PV, PF, OV, OF);
523 assert(res);
524}
525
527
528void polyfem::mesh::Polygons::polar_split(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector<std::vector<int>> &OF, double t)
529{
530 assert(IV.cols() == 2 || IV.cols() == 3);
531 Eigen::RowVector3d bary;
532 if (is_star_shaped(IV, bary))
533 {
534 if (t == 0.0)
535 {
536 // Special case 1: collapse the remaining polygon into a single vertex
537 OV.resize(IV.rows() + 1, IV.cols());
538 OV.topRows(IV.rows()) = IV;
539 OV.bottomRows(1) = bary.head(IV.cols());
540 int n = (int)IV.rows();
541 for (int i = 0; i < IV.rows(); ++i)
542 {
543 OF.push_back({i, (i + 1) % n, n});
544 }
545 }
546 else if (t < 0.0)
547 {
548 throw std::invalid_argument("Value t should be >= 0.0.");
549 }
550 else if (t >= 1.0)
551 {
552 throw std::invalid_argument("Value t >= 1.0 would create degenerate elements.");
553 }
554 else
555 {
556 // Create 1-ring around the central point
557 OV.resize(2 * IV.rows(), IV.cols());
558 OV.topRows(IV.rows()) = IV;
559 OV.bottomRows(IV.rows()) = (t * IV).rowwise() + (1.0 - t) * bary.head(IV.cols());
560 int n = (int)IV.rows();
561 OF.push_back({});
562 for (int i = 0; i < IV.rows(); ++i)
563 {
564 OF.front().push_back(i + n);
565 OF.push_back({i, (i + 1) % n, (i + 1) % n + n, i + n});
566 }
567 }
568 }
569 else
570 {
571 throw std::invalid_argument("Non star-shaped input polygon.");
572 }
573}
574
575// -----------------------------------------------------------------------------
576
577void polyfem::mesh::Polygons::catmul_clark_split(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector<std::vector<int>> &OF)
578{
579 assert(IV.cols() == 2 || IV.cols() == 3);
580 assert(IV.rows() % 2 == 0);
581 Eigen::RowVector3d bary;
582 if (is_star_shaped(IV, bary))
583 {
584 // Create 1-ring around the central point
585 OV.resize(IV.rows() + 1, IV.cols());
586 OV.topRows(IV.rows()) = IV;
587 OV.bottomRows(1) = bary.head(IV.cols());
588 int n = (int)IV.rows();
589 for (int i = 1; i < IV.rows(); i += 2)
590 {
591 OF.push_back({i, (i + 1) % n, (i + 2) % n, n});
592 }
593 }
594 else
595 {
596 throw std::invalid_argument("Non star-shaped input polygon.");
597 }
598}
599
600// -----------------------------------------------------------------------------
601
602void polyfem::mesh::Polygons::no_split(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector<std::vector<int>> &OF)
603{
604 OV = IV;
605 OF.clear();
606 OF.push_back({});
607 OF.front().resize(IV.rows());
608 std::iota(OF.front().begin(), OF.front().end(), 0);
609}
610
612
613void polyfem::mesh::refine_polygonal_mesh(const GEO::Mesh &M_in, GEO::Mesh &M_out, Polygons::SplitFunction split_func)
614{
615 using GEO::index_t;
616 using Navigation::Index;
617
618 // Step 0: Clear output mesh, and fill it with M_in's vertices
619 assert(&M_in != &M_out);
620 M_out.copy(M_in);
621 M_out.edges.clear();
622 M_out.facets.clear();
623 GEO::Attribute<GEO::index_t> c2e(M_in.facet_corners.attributes(), "edge_id");
624
625 // Step 1: Iterate over facets and refine triangles and quads
626 std::vector<int> edge_to_midpoint(M_in.edges.nb(), -1);
627 for (index_t f = 0; f < M_in.facets.nb(); ++f)
628 {
629 index_t nv = M_in.facets.nb_vertices(f);
630 assert(nv > 2);
631 Index idx = Navigation::get_index_from_face(M_in, c2e, f, 0);
632 assert(Navigation::switch_vertex(M_in, idx).vertex == (int)M_in.facets.vertex(f, 1));
633
634 // Create mid-edge vertices
635 for (index_t lv = 0; lv < M_in.facets.nb_vertices(f); ++lv, idx = Navigation::next_around_face(M_in, c2e, idx))
636 {
637 assert(idx.vertex == (int)M_in.facets.vertex(f, lv));
638 if (edge_to_midpoint[idx.edge] == -1)
639 {
640 GEO::vec3 coords = 0.5 * (mesh_vertex(M_in, idx.vertex) + mesh_vertex(M_in, Navigation::switch_vertex(M_in, idx).vertex));
641 edge_to_midpoint[idx.edge] = mesh_create_vertex(M_out, coords);
642 assert(edge_to_midpoint[idx.edge] + 1 == (int)M_out.vertices.nb());
643 }
644 }
645
646 if (/*nv == 3 ||*/ nv == 4)
647 {
648 // Create mid-face vertex
649 index_t vf = mesh_create_vertex(M_out, facet_barycenter(M_in, f));
650 assert(vf + 1 == M_out.vertices.nb());
651
652 // Create quads
653 for (index_t lv = 0; lv < M_in.facets.nb_vertices(f); ++lv)
654 {
655 idx = Navigation::get_index_from_face(M_in, c2e, f, lv);
656 assert(Navigation::switch_vertex(M_in, idx).vertex == (int)M_in.facets.vertex(f, (lv + 1) % M_in.facets.nb_vertices(f)));
657 int v1 = idx.vertex;
658 int v12 = edge_to_midpoint[idx.edge];
659 int v01 = edge_to_midpoint[Navigation::switch_edge(M_in, c2e, idx).edge];
660 assert(v12 != -1 && v01 != -1);
661 M_out.facets.create_quad(v1, v12, vf, v01);
662 }
663 }
664 }
665
666 // Step 2: Create polygonal faces following vertices around holes
667 for (index_t f = 0; f < M_in.facets.nb(); ++f)
668 {
669 if (M_in.facets.nb_vertices(f) <= 4)
670 {
671 continue;
672 }
673 GEO::vector<index_t> hole;
674 {
675 auto idx = Navigation::get_index_from_face(M_in, c2e, f, 0);
676 for (index_t lv = 0; lv < M_in.facets.nb_vertices(f); ++lv)
677 {
678 hole.push_back(idx.vertex);
679 hole.push_back(edge_to_midpoint[idx.edge]);
680 idx = Navigation::next_around_face(M_in, c2e, idx);
681 }
682 }
683
684 if (M_in.vertices.dimension() != 2)
685 {
686 // std::cerr << "WARNING: Input mesh has dimension > 2, but polygonal facets will be split considering their XY coordinates only." << std::endl;
687 }
688
689 // Subdivide the hole using polar refinement
690 index_t n = hole.size();
691 Eigen::MatrixXd P(n, 2), V;
692 std::vector<std::vector<int>> F;
693 for (index_t k = 0; k < n; ++k)
694 {
695 GEO::vec3 pk = polyfem::mesh::mesh_vertex(M_out, hole[k]);
696 P.row(k) << pk[0], pk[1];
697 }
698 split_func(P, V, F);
699 assert(V.rows() >= n);
700 std::vector<int> remap(V.rows() - n);
701 for (index_t k = n, lk = 0; k < V.rows(); ++k, ++lk)
702 {
703 GEO::vec3 qk = GEO::vec3(V(k, 0), V(k, 1), 0);
704 remap[lk] = mesh_create_vertex(M_out, qk);
705 }
706 for (const auto &poly : F)
707 {
708 GEO::vector<index_t> vertices;
709 for (int vk : poly)
710 {
711 if (vk < (int)n)
712 {
713 vertices.push_back(hole[vk]);
714 }
715 else
716 {
717 vertices.push_back(remap[vk - n]);
718 }
719 }
720 M_out.facets.create_polygon(vertices);
721 }
722 }
723}
724
726
727void polyfem::mesh::refine_triangle_mesh(const GEO::Mesh &M_in, GEO::Mesh &M_out)
728{
729 using GEO::index_t;
730 using Navigation::Index;
731
732 // Step 1: Clear output mesh, and fill it with M_in's vertices
733 assert(&M_in != &M_out);
734 M_out.copy(M_in);
735 M_out.edges.clear();
736 M_out.facets.clear();
737 GEO::Attribute<GEO::index_t> c2e(M_in.facet_corners.attributes(), "edge_id");
738
739 // Step 2: Iterate over facets and refine triangles
740 std::vector<int> edge_to_midpoint(M_in.edges.nb(), -1);
741 for (index_t f = 0; f < M_in.facets.nb(); ++f)
742 {
743 index_t nv = M_in.facets.nb_vertices(f);
744 assert(nv == 3);
745 Index idx = Navigation::get_index_from_face(M_in, c2e, f, 0);
746 assert(Navigation::switch_vertex(M_in, idx).vertex == (int)M_in.facets.vertex(f, 1));
747
748 // Create mid-edge vertices
749 for (index_t lv = 0; lv < M_in.facets.nb_vertices(f); ++lv, idx = Navigation::next_around_face(M_in, c2e, idx))
750 {
751 assert(idx.vertex == (int)M_in.facets.vertex(f, lv));
752 if (edge_to_midpoint[idx.edge] == -1)
753 {
754 GEO::vec3 coords = 0.5 * (mesh_vertex(M_in, idx.vertex) + mesh_vertex(M_in, Navigation::switch_vertex(M_in, idx).vertex));
755 edge_to_midpoint[idx.edge] = mesh_create_vertex(M_out, coords);
756 assert(edge_to_midpoint[idx.edge] + 1 == (int)M_out.vertices.nb());
757 }
758 }
759
760 // Create triangles
761 std::array<index_t, 3> e2v;
762 for (index_t lv = 0; lv < M_in.facets.nb_vertices(f); ++lv)
763 {
764 idx = Navigation::get_index_from_face(M_in, c2e, f, lv);
765 assert(Navigation::switch_vertex(M_in, idx).vertex == (int)M_in.facets.vertex(f, (lv + 1) % M_in.facets.nb_vertices(f)));
766 int v1 = idx.vertex;
767 int v12 = edge_to_midpoint[idx.edge];
768 int v01 = edge_to_midpoint[Navigation::switch_edge(M_in, c2e, idx).edge];
769 e2v[lv] = v12;
770 assert(v12 != -1 && v01 != -1);
771 M_out.facets.create_triangle(v1, v12, v01);
772 }
773 M_out.facets.create_triangle(e2v[0], e2v[1], e2v[2]);
774 }
775}
vector< list< int > > adj
int V
int y
int x
Index get_index_from_face(const GEO::Mesh &M, const GEO::Attribute< GEO::index_t > &c2e, int f, int lv)
Index switch_edge(const GEO::Mesh &M, const GEO::Attribute< GEO::index_t > &c2e, Index idx)
Index next_around_face(const GEO::Mesh &M, const GEO::Attribute< GEO::index_t > &c2e, Index idx)
Index switch_vertex(const GEO::Mesh &M, Index idx)
void polar_split(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector< std::vector< int > > &OF, double t=0.5)
Split a polygon using polar refinement.
void catmul_clark_split(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector< std::vector< int > > &OF)
Split a polygon using polar refinement.
std::function< void(Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector< std::vector< int > > &OF)> SplitFunction
void no_split(const Eigen::MatrixXd &IV, Eigen::MatrixXd &OV, std::vector< std::vector< int > > &OF)
Don't split polygons.
void refine_triangle_mesh(const GEO::Mesh &M_in, GEO::Mesh &M_out)
Refine a triangle mesh.
void refine_quad_mesh(const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IF, Eigen::MatrixXd &OV, Eigen::MatrixXi &OF)
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 edge_adjacency_graph(const Eigen::MatrixXi &Q, Eigen::MatrixXi &edge_index, std::vector< std::vector< int > > &adj, std::vector< std::pair< int, int > > *pairs_of_edges=nullptr, std::vector< std::pair< int, int > > *pairs_of_quads=nullptr, Eigen::MatrixXi *quad_index=nullptr)
GEO::vec3 facet_barycenter(const GEO::Mesh &M, GEO::index_t f)
Definition MeshUtils.cpp:64
bool instantiate_pattern(const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IF, const Eigen::MatrixXd &PV, const Eigen::MatrixXi &PF, Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::VectorXi *SF=nullptr, EvalParametersFunc evalFunc=nullptr, GetAdjacentLocalEdge getAdjLocalEdge=nullptr)
std::function< void(const Eigen::MatrixXd &, Eigen::MatrixXd &, int)> EvalParametersFunc
void refine_polygonal_mesh(const GEO::Mesh &M_in, GEO::Mesh &M_out, Polygons::SplitFunction split_func)
Refine a polygonal mesh.
std::function< std::tuple< int, int, bool >(int, int)> GetAdjacentLocalEdge
bool is_star_shaped(const Eigen::MatrixXd &IV, Eigen::RowVector3d &bary)
Determine whether a polygon is star-shaped or not.
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44