PolyFEM
Loading...
Searching...
No Matches
RefElementSampler.cpp
Go to the documentation of this file.
1#include "RefElementSampler.hpp"
2#include <polyfem/mesh/MeshUtils.hpp>
3
4#include <polyfem/quadrature/PolyhedronQuadrature.hpp>
5
6#include <igl/predicates/ear_clipping.h>
7
8#include <igl/edges.h>
9#include <igl/boundary_loop.h>
10#include <igl/doublearea.h>
11#include <igl/upsample.h>
12
13#ifdef POLYFEM_WITH_TRIANGLE
14#include <igl/triangle/triangulate.h>
15#endif
16
17#include <cassert>
18#include <cmath>
19
20namespace polyfem
21{
22 using namespace mesh;
23 namespace utils
24 {
25
34 void regular_2d_grid(const int n, bool tri, Eigen::MatrixXd &V, Eigen::MatrixXi &F)
35 {
36
37 V.resize(n * n, 2);
38 F.resize((n - 1) * (n - 1) * 2, 3);
39 const double delta = 1. / (n - 1.);
40 std::vector<int> map(n * n, -1);
41
42 int index = 0;
43 for (int i = 0; i < n; ++i)
44 {
45 for (int j = 0; j < n; ++j)
46 {
47 if (tri && i + j >= n)
48 continue;
49 map[i + j * n] = index;
50 V.row(index) << i * delta, j * delta;
51 ++index;
52 }
53 }
54
55 V.conservativeResize(index, 2);
56
57 std::array<int, 3> tmp;
58
59 index = 0;
60 for (int i = 0; i < n - 1; ++i)
61 {
62 for (int j = 0; j < n - 1; ++j)
63 {
64 tmp = {{map[i + j * n], map[i + 1 + j * n], map[i + (j + 1) * n]}};
65 if (tmp[0] >= 0 && tmp[1] >= 0 && tmp[2] >= 0)
66 {
67 F.row(index) << tmp[0], tmp[1], tmp[2];
68 ++index;
69 }
70
71 tmp = {{map[i + 1 + j * n], map[i + 1 + (j + 1) * n], map[i + (j + 1) * n]}};
72 if (tmp[0] >= 0 && tmp[1] >= 0 && tmp[2] >= 0)
73 {
74 F.row(index) << tmp[0], tmp[1], tmp[2];
75 ++index;
76 }
77 }
78 }
79
80 F.conservativeResize(index, 3);
81 }
82
83 namespace
84 {
85 constexpr double PYRAMID_APEX_EPS = 1e-8;
86
87 void avoid_pyramid_apex(Eigen::MatrixXd &V)
88 {
89 assert(V.cols() == 3);
90 for (int i = 0; i < V.rows(); ++i)
91 {
92 if (std::abs(V(i, 2) - 1.0) < PYRAMID_APEX_EPS)
93 V(i, 2) = 1.0 - PYRAMID_APEX_EPS;
94 }
95 }
96
97 void add_tet(const std::array<int, 4> &tmp, const Eigen::MatrixXd &V, int &index, Eigen::MatrixXi &T)
98 {
99 if (tmp[0] >= 0 && tmp[1] >= 0 && tmp[2] >= 0 && tmp[3] >= 0)
100 {
101 const Eigen::Vector3d e0 = V.row(tmp[1]) - V.row(tmp[0]);
102 const Eigen::Vector3d e1 = V.row(tmp[2]) - V.row(tmp[0]);
103 const Eigen::Vector3d e2 = V.row(tmp[3]) - V.row(tmp[0]);
104 double vol = (e0.cross(e1)).dot(e2);
105 if (std::abs(vol) < 1e-14)
106 return;
107 if (vol < 0)
108 T.row(index) << tmp[0], tmp[1], tmp[2], tmp[3];
109 else
110 T.row(index) << tmp[0], tmp[1], tmp[3], tmp[2];
111#ifndef NDEBUG
112 const Eigen::Vector3d ed0 = V.row(T(index, 1)) - V.row(T(index, 0));
113 const Eigen::Vector3d ed1 = V.row(T(index, 2)) - V.row(T(index, 0));
114 const Eigen::Vector3d ed2 = V.row(T(index, 3)) - V.row(T(index, 0));
115 assert((ed0.cross(ed1)).dot(ed2) < 0);
116#endif
117 ++index;
118 }
119 }
120
121 } // anonymous namespace
122
134 void regular_3d_grid(const int nn,
135 bool tet,
136 bool prism,
137 bool pyramid,
138 Eigen::MatrixXd &V,
139 Eigen::MatrixXi &F,
140 Eigen::MatrixXi &T)
141 {
142 const int n = nn;
143 const double delta = 1. / (n - 1.);
144
145 if (pyramid)
146 {
147 T.resize((n - 1) * (n - 1) * (n - 1) * 6, 4);
148 V.resize(n * n * (n - 1) + 1, 3);
149 std::vector<int> map(n * n * n, -1);
150
151 int index = 0;
152 int apex_index = -1;
153 for (int i = 0; i < n; ++i)
154 {
155 for (int j = 0; j < n; ++j)
156 {
157 for (int k = 0; k < n; ++k)
158 {
159 const int grid_index = (i + j * n) * n + k;
160 if (k == n - 1)
161 {
162 if (apex_index < 0)
163 {
164 apex_index = index++;
165 V.row(apex_index) << 0, 0, 1;
166 }
167 map[grid_index] = apex_index;
168 continue;
169 }
170
171 const double z = k * delta;
172 const double scale = 1.0 - z;
173 map[grid_index] = index;
174 V.row(index) << i * delta * scale, j * delta * scale, z;
175 ++index;
176 }
177 }
178 }
179 V.conservativeResize(index, 3);
180 avoid_pyramid_apex(V);
181
182 std::array<int, 8> indices;
183 std::array<int, 4> tmp;
184 index = 0;
185 for (int i = 0; i < n - 1; ++i)
186 {
187 for (int j = 0; j < n - 1; ++j)
188 {
189 for (int k = 0; k < n - 1; ++k)
190 {
191 indices = {{(i + j * n) * n + k,
192 (i + 1 + j * n) * n + k,
193 (i + 1 + (j + 1) * n) * n + k,
194 (i + (j + 1) * n) * n + k,
195
196 (i + j * n) * n + k + 1,
197 (i + 1 + j * n) * n + k + 1,
198 (i + 1 + (j + 1) * n) * n + k + 1,
199 (i + (j + 1) * n) * n + k + 1}};
200
201 tmp = {{map[indices[1 - 1]], map[indices[2 - 1]], map[indices[4 - 1]], map[indices[5 - 1]]}};
202 add_tet(tmp, V, index, T);
203
204 tmp = {{map[indices[6 - 1]], map[indices[3 - 1]], map[indices[7 - 1]], map[indices[8 - 1]]}};
205 add_tet(tmp, V, index, T);
206
207 tmp = {{map[indices[5 - 1]], map[indices[2 - 1]], map[indices[6 - 1]], map[indices[4 - 1]]}};
208 add_tet(tmp, V, index, T);
209
210 tmp = {{map[indices[5 - 1]], map[indices[4 - 1]], map[indices[8 - 1]], map[indices[6 - 1]]}};
211 add_tet(tmp, V, index, T);
212
213 tmp = {{map[indices[4 - 1]], map[indices[2 - 1]], map[indices[6 - 1]], map[indices[3 - 1]]}};
214 add_tet(tmp, V, index, T);
215
216 tmp = {{map[indices[3 - 1]], map[indices[4 - 1]], map[indices[8 - 1]], map[indices[6 - 1]]}};
217 add_tet(tmp, V, index, T);
218 }
219 }
220 }
221
222 T.conservativeResize(index, 4);
223
224 F.resize(4 * index, 3);
225
226 F.block(0, 0, index, 1) = T.col(1);
227 F.block(0, 1, index, 1) = T.col(0);
228 F.block(0, 2, index, 1) = T.col(2);
229
230 F.block(index, 0, index, 1) = T.col(0);
231 F.block(index, 1, index, 1) = T.col(1);
232 F.block(index, 2, index, 1) = T.col(3);
233
234 F.block(2 * index, 0, index, 1) = T.col(1);
235 F.block(2 * index, 1, index, 1) = T.col(2);
236 F.block(2 * index, 2, index, 1) = T.col(3);
237
238 F.block(3 * index, 0, index, 1) = T.col(2);
239 F.block(3 * index, 1, index, 1) = T.col(0);
240 F.block(3 * index, 2, index, 1) = T.col(3);
241
242 return;
243 }
244
245 T.resize((n - 1) * (n - 1) * (n - 1) * 6, 4);
246 V.resize(n * n * n, 3);
247 std::vector<int> map(n * n * n, -1);
248
249 int index = 0;
250 for (int i = 0; i < n; ++i)
251 {
252 for (int j = 0; j < n; ++j)
253 {
254 for (int k = 0; k < n; ++k)
255 {
256 if (tet && i + j + k >= n)
257 continue;
258 if (prism && i + j >= n)
259 continue;
260 // TODO
261 map[(i + j * n) * n + k] = index;
262 V.row(index) << i * delta, j * delta, k * delta;
263 ++index;
264 }
265 }
266 }
267 V.conservativeResize(index, 3);
268
269 std::array<int, 8> indices;
270 std::array<int, 4> tmp;
271 index = 0;
272 for (int i = 0; i < n - 1; ++i)
273 {
274 for (int j = 0; j < n - 1; ++j)
275 {
276 for (int k = 0; k < n - 1; ++k)
277 {
278 indices = {{(i + j * n) * n + k,
279 (i + 1 + j * n) * n + k,
280 (i + 1 + (j + 1) * n) * n + k,
281 (i + (j + 1) * n) * n + k,
282
283 (i + j * n) * n + k + 1,
284 (i + 1 + j * n) * n + k + 1,
285 (i + 1 + (j + 1) * n) * n + k + 1,
286 (i + (j + 1) * n) * n + k + 1}};
287
288 tmp = {{map[indices[1 - 1]], map[indices[2 - 1]], map[indices[4 - 1]], map[indices[5 - 1]]}};
289 add_tet(tmp, V, index, T);
290
291 tmp = {{map[indices[6 - 1]], map[indices[3 - 1]], map[indices[7 - 1]], map[indices[8 - 1]]}};
292 add_tet(tmp, V, index, T);
293
294 tmp = {{map[indices[5 - 1]], map[indices[2 - 1]], map[indices[6 - 1]], map[indices[4 - 1]]}};
295 add_tet(tmp, V, index, T);
296
297 tmp = {{map[indices[5 - 1]], map[indices[4 - 1]], map[indices[8 - 1]], map[indices[6 - 1]]}};
298 add_tet(tmp, V, index, T);
299
300 tmp = {{map[indices[4 - 1]], map[indices[2 - 1]], map[indices[6 - 1]], map[indices[3 - 1]]}};
301 add_tet(tmp, V, index, T);
302
303 tmp = {{map[indices[3 - 1]], map[indices[4 - 1]], map[indices[8 - 1]], map[indices[6 - 1]]}};
304 add_tet(tmp, V, index, T);
305 }
306 }
307 }
308
309 T.conservativeResize(index, 4);
310
311 F.resize(4 * index, 3);
312
313 F.block(0, 0, index, 1) = T.col(1);
314 F.block(0, 1, index, 1) = T.col(0);
315 F.block(0, 2, index, 1) = T.col(2);
316
317 F.block(index, 0, index, 1) = T.col(0);
318 F.block(index, 1, index, 1) = T.col(1);
319 F.block(index, 2, index, 1) = T.col(3);
320
321 F.block(2 * index, 0, index, 1) = T.col(1);
322 F.block(2 * index, 1, index, 1) = T.col(2);
323 F.block(2 * index, 2, index, 1) = T.col(3);
324
325 F.block(3 * index, 0, index, 1) = T.col(2);
326 F.block(3 * index, 1, index, 1) = T.col(0);
327 F.block(3 * index, 2, index, 1) = T.col(3);
328 }
329
330 void RefElementSampler::init(const bool is_volume, const int n_elements, double target_rel_area)
331 {
332 is_volume_ = is_volume;
333
334 area_param_ = target_rel_area * n_elements;
335#ifndef NDEBUG
336 area_param_ *= 10.0;
337#endif
338
339 build();
340 }
341
342 void RefElementSampler::build()
343 {
344 using namespace Eigen;
345
346 if (is_volume_)
347 {
348 // cube
349 {
350 MatrixXd pts(8, 3);
351 pts << 0, 0, 0,
352 0, 1, 0,
353 1, 1, 0,
354 1, 0, 0,
355
356 // 4
357 0, 0, 1,
358 0, 1, 1,
359 1, 1, 1,
360 1, 0, 1;
361
362 Eigen::MatrixXi faces(12, 3);
363 faces << 1, 2, 0,
364 0, 2, 3,
365
366 5, 4, 6,
367 4, 7, 6,
368
369 1, 0, 4,
370 1, 4, 5,
371
372 2, 1, 5,
373 2, 5, 6,
374
375 3, 2, 6,
376 3, 6, 7,
377
378 0, 3, 7,
379 0, 7, 4;
380
381 regular_3d_grid(std::max(2., round(1. / pow(area_param_, 1. / 3.) + 1) / 2.), false, false, false, cube_points_, cube_faces_, cube_tets_);
382
383 // Extract sampled edges matching the base element edges
384 Eigen::MatrixXi edges(12, 2);
385 edges << 0, 1,
386 1, 2,
387 2, 3,
388 3, 0,
389
390 4, 5,
391 5, 6,
392 6, 7,
393 7, 4,
394
395 0, 4,
396 1, 5,
397 2, 6,
398 3, 7;
399 igl::edges(cube_faces_, cube_edges_);
400 extract_parent_edges(cube_points_, cube_edges_, pts, edges, cube_edges_);
401
402 // Same local order as in FEMBasis3d
403 cube_corners_.resize(8, 3);
404 cube_corners_ << 0, 0, 0,
405 1, 0, 0,
406 1, 1, 0,
407 0, 1, 0,
408 0, 0, 1,
409 1, 0, 1,
410 1, 1, 1,
411 0, 1, 1;
412 }
413
414 // tet
415 {
416 MatrixXd pts(4, 3);
417 pts << 0, 0, 0,
418 1, 0, 0,
419 0, 1, 0,
420 0, 0, 1;
421
422 Eigen::MatrixXi faces(4, 3);
423 faces << 0, 1, 2,
424
425 3, 1, 0,
426 2, 1, 3,
427 0, 2, 3;
428
429 regular_3d_grid(std::max(2., round(1. / pow(area_param_, 1. / 3.) + 1)), true, false, false, simplex_points_, simplex_faces_, simplex_tets_);
430
431 // Extract sampled edges matching the base element edges
432 Eigen::MatrixXi edges;
433 igl::edges(faces, edges);
434 igl::edges(simplex_faces_, simplex_edges_);
435 extract_parent_edges(simplex_points_, simplex_edges_, pts, edges, simplex_edges_);
436
437 // Same local order as in FEMBasis3d
438 simplex_corners_.resize(4, 3);
439 simplex_corners_ << 0, 0, 0,
440 1, 0, 0,
441 0, 1, 0,
442 0, 0, 1;
443 }
444
445 // prism
446 {
447 MatrixXd pts(6, 3);
448 pts << 0, 0, 0,
449 0, 1, 0,
450 1, 0, 0,
451 0, 0, 1,
452 0, 1, 1,
453 1, 0, 1;
454
455 regular_3d_grid(std::max(2., round(1. / pow(area_param_, 1. / 3.) + 1) / 2.),
456 false, true, false, prism_points_, prism_faces_, prism_tets_);
457
458 // Extract sampled edges matching the base element edges
459 Eigen::MatrixXi edges(9, 2);
460 edges << 0, 1,
461 1, 2,
462 2, 0,
463
464 3, 4,
465 4, 5,
466 5, 3,
467
468 0, 3,
469 1, 4,
470 2, 5;
471
472 igl::edges(prism_faces_, prism_edges_);
473 extract_parent_edges(prism_points_, prism_edges_, pts, edges, prism_edges_);
474
475 // Same local order as in FEMBasis3d
476 prism_corners_.resize(6, 3);
477 prism_corners_ << 0, 0, 0,
478 0, 1, 0,
479 1, 0, 0,
480 0, 0, 1,
481 0, 1, 1,
482 1, 0, 1;
483 }
484
485 // pyramid
486 {
487 MatrixXd pts(5, 3);
488 pts << 0, 0, 0,
489 1, 0, 0,
490 1, 1, 0,
491 0, 1, 0,
492 0, 0, 1;
493
494 regular_3d_grid(std::max(2., round(1. / pow(area_param_, 1. / 3.) + 1) / 2.),
495 false, false, true, pyramid_points_, pyramid_faces_, pyramid_tets_);
496
497 // Extract sampled edges matching the base element edges
498 Eigen::MatrixXi edges(8, 2);
499 edges << 0, 1,
500 1, 2,
501 2, 3,
502 3, 0,
503 0, 4,
504 1, 4,
505 2, 4,
506 3, 4;
507
508 igl::edges(pyramid_faces_, pyramid_edges_);
509 Eigen::MatrixXd edge_pts = pts;
510 avoid_pyramid_apex(edge_pts);
511 extract_parent_edges(pyramid_points_, pyramid_edges_, edge_pts, edges, pyramid_edges_);
512
513 // Same local order as in FEMBasis3d
514 pyramid_corners_.resize(5, 3);
515 pyramid_corners_ << 0, 0, 0,
516 1, 0, 0,
517 1, 1, 0,
518 0, 1, 0,
519 0, 0, 1;
520 }
521 }
522 else
523 {
524 {
525 MatrixXd pts(4, 2);
526 pts << 0, 0,
527 0, 1,
528 1, 1,
529 1, 0;
530
531 MatrixXi E(4, 2);
532 E << 0, 1,
533 1, 2,
534 2, 3,
535 3, 0;
536
537 regular_2d_grid(std::max(2., round(1. / sqrt(area_param_) + 1)), false, cube_points_, cube_faces_);
538
539 // Extract sampled edges matching the base element edges
540 igl::edges(cube_faces_, cube_edges_);
541 extract_parent_edges(cube_points_, cube_edges_, pts, E, cube_edges_);
542
543 // Same local order as in FEMBasis2d
544 cube_corners_.resize(4, 2);
545 cube_corners_ << 0, 0,
546 1, 0,
547 1, 1,
548 0, 1;
549 }
550 {
551 MatrixXd pts(3, 2);
552 pts << 0, 0,
553 0, 1,
554 1, 0;
555
556 MatrixXi E(3, 2);
557 E << 0, 1,
558 1, 2,
559 2, 0;
560
561 regular_2d_grid(std::max(2., round(1. / sqrt(area_param_) + 1)), true, simplex_points_, simplex_faces_);
562
563 // Extract sampled edges matching the base element edges
564 igl::edges(simplex_faces_, simplex_edges_);
565 extract_parent_edges(simplex_points_, simplex_edges_, pts, E, simplex_edges_);
566
567 // Same local order as in FEMBasis2d
568 simplex_corners_.resize(3, 2);
569 simplex_corners_ << 0, 0,
570 1, 0,
571 0, 1;
572 }
573 }
574 }
575
576 void RefElementSampler::sample_polygon(const Eigen::MatrixXd &poly, Eigen::MatrixXd &pts, Eigen::MatrixXi &faces, Eigen::MatrixXi &edges) const
577 {
578
579#ifdef POLYFEM_WITH_TRIANGLE
580 Eigen::MatrixXi E(poly.rows(), 2);
581 const Eigen::MatrixXd H(0, 2);
582 const std::string flags = "Qzqa" + std::to_string(area_param_ / 10.0);
583
584 for (int i = 0; i < poly.rows(); ++i)
585 E.row(i) << i, (i + 1) % poly.rows();
586
587 igl::triangle::triangulate(poly, E, H, flags, pts, faces);
588#else
589
590 const Eigen::MatrixXi rt = Eigen::MatrixXi::Zero(poly.rows(), 1);
591 faces.resize(0, 0);
592 Eigen::VectorXi I;
593 Eigen::MatrixXd area;
594 igl::predicates::ear_clipping(poly, rt, faces, I);
595
596 igl::doublearea(poly, faces, area);
597
598 const double area_avg = area.array().sum() / poly.rows() / 2;
599
600 const int n_refs = area_avg / area_param_ * 40;
601
602 Eigen::MatrixXi new_faces;
603
604 igl::upsample(poly, faces, pts, new_faces, n_refs);
605
606 faces = new_faces;
607#endif
608 std::vector<int> loop;
609 igl::default_num_threads(1);
610 igl::boundary_loop(faces, loop);
611 igl::default_num_threads(0);
612 edges.resize(loop.size(), 2);
613 for (int i = 0; i < loop.size(); ++i)
614 edges.row(i) << loop[i], loop[(i + 1) % loop.size()];
615 }
616
617 void RefElementSampler::sample_polyhedron(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &f, Eigen::MatrixXd &pts, Eigen::MatrixXi &tets, Eigen::MatrixXi &faces) const
618 {
619 const Eigen::MatrixXd kernel = vertices.colwise().mean();
620
621 polyfem::tertrahedralize_star_shaped_surface(vertices, f, kernel, pts, faces, tets);
622 }
623 } // namespace utils
624} // namespace polyfem
V
int V
Definition AdjointNLProblem.cpp:47
faces
std::vector< Eigen::VectorXi > faces
Definition PressureAssembler.cpp:52
z
int z
Definition SplineBasis3d.cpp:55
polyfem::utils::RefElementSampler::init
void init(const bool is_volume, const int n_elements, const double target_rel_area)
Definition RefElementSampler.cpp:330
polyfem::utils::RefElementSampler::sample_polygon
void sample_polygon(const Eigen::MatrixXd &poly, Eigen::MatrixXd &pts, Eigen::MatrixXi &faces, Eigen::MatrixXi &edges) const
Definition RefElementSampler.cpp:576
polyfem::utils::RefElementSampler::sample_polyhedron
void sample_polyhedron(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &f, Eigen::MatrixXd &pts, Eigen::MatrixXi &faces, Eigen::MatrixXi &edges) const
Definition RefElementSampler.cpp:617
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.tmp
list tmp
Definition p_bases.py:366
polyfem::mesh::extract_parent_edges
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...
Definition MeshUtils.cpp:1147
polyfem::utils::regular_2d_grid
void regular_2d_grid(const int n, bool tri, Eigen::MatrixXd &V, Eigen::MatrixXi &F)
Generate a canonical triangle/quad subdivided from a regular grid.
Definition RefElementSampler.cpp:34
polyfem::utils::regular_3d_grid
void regular_3d_grid(const int nn, bool tet, bool prism, bool pyramid, Eigen::MatrixXd &V, Eigen::MatrixXi &F, Eigen::MatrixXi &T)
Generate a canonical tet/hex subdivided from a regular grid.
Definition RefElementSampler.cpp:134
polyfem::tertrahedralize_star_shaped_surface
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.
Definition PolyhedronQuadrature.cpp:17