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
19namespace polyfem
20{
21 using namespace mesh;
22 namespace utils
23 {
24
33 void regular_2d_grid(const int n, bool tri, Eigen::MatrixXd &V, Eigen::MatrixXi &F)
34 {
35
36 V.resize(n * n, 2);
37 F.resize((n - 1) * (n - 1) * 2, 3);
38 const double delta = 1. / (n - 1.);
39 std::vector<int> map(n * n, -1);
40
41 int index = 0;
42 for (int i = 0; i < n; ++i)
43 {
44 for (int j = 0; j < n; ++j)
45 {
46 if (tri && i + j >= n)
47 continue;
48 map[i + j * n] = index;
49 V.row(index) << i * delta, j * delta;
50 ++index;
51 }
52 }
53
54 V.conservativeResize(index, 2);
55
56 std::array<int, 3> tmp;
57
58 index = 0;
59 for (int i = 0; i < n - 1; ++i)
60 {
61 for (int j = 0; j < n - 1; ++j)
62 {
63 tmp = {{map[i + j * n], map[i + 1 + j * n], map[i + (j + 1) * n]}};
64 if (tmp[0] >= 0 && tmp[1] >= 0 && tmp[2] >= 0)
65 {
66 F.row(index) << tmp[0], tmp[1], tmp[2];
67 ++index;
68 }
69
70 tmp = {{map[i + 1 + j * n], map[i + 1 + (j + 1) * n], map[i + (j + 1) * n]}};
71 if (tmp[0] >= 0 && tmp[1] >= 0 && tmp[2] >= 0)
72 {
73 F.row(index) << tmp[0], tmp[1], tmp[2];
74 ++index;
75 }
76 }
77 }
78
79 F.conservativeResize(index, 3);
80 }
81
82 namespace
83 {
84
85 void add_tet(const std::array<int, 4> &tmp, const Eigen::MatrixXd &V, int &index, Eigen::MatrixXi &T)
86 {
87 if (tmp[0] >= 0 && tmp[1] >= 0 && tmp[2] >= 0 && tmp[3] >= 0)
88 {
89 const Eigen::Vector3d e0 = V.row(tmp[1]) - V.row(tmp[0]);
90 const Eigen::Vector3d e1 = V.row(tmp[2]) - V.row(tmp[0]);
91 const Eigen::Vector3d e2 = V.row(tmp[3]) - V.row(tmp[0]);
92 double vol = (e0.cross(e1)).dot(e2);
93 if (vol < 0)
94 T.row(index) << tmp[0], tmp[1], tmp[2], tmp[3];
95 else
96 T.row(index) << tmp[0], tmp[1], tmp[3], tmp[2];
97#ifndef NDEBUG
98 const Eigen::Vector3d ed0 = V.row(T(index, 1)) - V.row(T(index, 0));
99 const Eigen::Vector3d ed1 = V.row(T(index, 2)) - V.row(T(index, 0));
100 const Eigen::Vector3d ed2 = V.row(T(index, 3)) - V.row(T(index, 0));
101 assert((ed0.cross(ed1)).dot(ed2) < 0);
102#endif
103 ++index;
104 }
105 }
106
107 } // anonymous namespace
108
120 void regular_3d_grid(const int nn,
121 bool tet,
122 bool prism,
123 bool pyramid,
124 Eigen::MatrixXd &V,
125 Eigen::MatrixXi &F,
126 Eigen::MatrixXi &T)
127 {
128 const int n = nn;
129 const double delta = 1. / (n - 1.);
130 T.resize((n - 1) * (n - 1) * (n - 1) * 6, 4);
131 V.resize(n * n * n, 3);
132 std::vector<int> map(n * n * n, -1);
133
134 int index = 0;
135 for (int i = 0; i < n; ++i)
136 {
137 for (int j = 0; j < n; ++j)
138 {
139 for (int k = 0; k < n; ++k)
140 {
141 if (tet && i + j + k >= n)
142 continue;
143 if (prism && i + j >= n)
144 continue;
145 if (pyramid && (i + k >= n || j + k >= n))
146 continue;
147 // TODO
148 map[(i + j * n) * n + k] = index;
149 V.row(index) << i * delta, j * delta, k * delta;
150 ++index;
151 }
152 }
153 }
154 V.conservativeResize(index, 3);
155
156 std::array<int, 8> indices;
157 std::array<int, 4> tmp;
158 index = 0;
159 for (int i = 0; i < n - 1; ++i)
160 {
161 for (int j = 0; j < n - 1; ++j)
162 {
163 for (int k = 0; k < n - 1; ++k)
164 {
165 indices = {{(i + j * n) * n + k,
166 (i + 1 + j * n) * n + k,
167 (i + 1 + (j + 1) * n) * n + k,
168 (i + (j + 1) * n) * n + k,
169
170 (i + j * n) * n + k + 1,
171 (i + 1 + j * n) * n + k + 1,
172 (i + 1 + (j + 1) * n) * n + k + 1,
173 (i + (j + 1) * n) * n + k + 1}};
174
175 tmp = {{map[indices[1 - 1]], map[indices[2 - 1]], map[indices[4 - 1]], map[indices[5 - 1]]}};
176 add_tet(tmp, V, index, T);
177
178 tmp = {{map[indices[6 - 1]], map[indices[3 - 1]], map[indices[7 - 1]], map[indices[8 - 1]]}};
179 add_tet(tmp, V, index, T);
180
181 tmp = {{map[indices[5 - 1]], map[indices[2 - 1]], map[indices[6 - 1]], map[indices[4 - 1]]}};
182 add_tet(tmp, V, index, T);
183
184 tmp = {{map[indices[5 - 1]], map[indices[4 - 1]], map[indices[8 - 1]], map[indices[6 - 1]]}};
185 add_tet(tmp, V, index, T);
186
187 tmp = {{map[indices[4 - 1]], map[indices[2 - 1]], map[indices[6 - 1]], map[indices[3 - 1]]}};
188 add_tet(tmp, V, index, T);
189
190 tmp = {{map[indices[3 - 1]], map[indices[4 - 1]], map[indices[8 - 1]], map[indices[6 - 1]]}};
191 add_tet(tmp, V, index, T);
192 }
193 }
194 }
195
196 T.conservativeResize(index, 4);
197
198 F.resize(4 * index, 3);
199
200 F.block(0, 0, index, 1) = T.col(1);
201 F.block(0, 1, index, 1) = T.col(0);
202 F.block(0, 2, index, 1) = T.col(2);
203
204 F.block(index, 0, index, 1) = T.col(0);
205 F.block(index, 1, index, 1) = T.col(1);
206 F.block(index, 2, index, 1) = T.col(3);
207
208 F.block(2 * index, 0, index, 1) = T.col(1);
209 F.block(2 * index, 1, index, 1) = T.col(2);
210 F.block(2 * index, 2, index, 1) = T.col(3);
211
212 F.block(3 * index, 0, index, 1) = T.col(2);
213 F.block(3 * index, 1, index, 1) = T.col(0);
214 F.block(3 * index, 2, index, 1) = T.col(3);
215 }
216
217 void RefElementSampler::init(const bool is_volume, const int n_elements, double target_rel_area)
218 {
219 is_volume_ = is_volume;
220
221 area_param_ = target_rel_area * n_elements;
222#ifndef NDEBUG
223 area_param_ *= 10.0;
224#endif
225
226 build();
227 }
228
229 void RefElementSampler::build()
230 {
231 using namespace Eigen;
232
233 if (is_volume_)
234 {
235 // cube
236 {
237 MatrixXd pts(8, 3);
238 pts << 0, 0, 0,
239 0, 1, 0,
240 1, 1, 0,
241 1, 0, 0,
242
243 // 4
244 0, 0, 1,
245 0, 1, 1,
246 1, 1, 1,
247 1, 0, 1;
248
249 Eigen::MatrixXi faces(12, 3);
250 faces << 1, 2, 0,
251 0, 2, 3,
252
253 5, 4, 6,
254 4, 7, 6,
255
256 1, 0, 4,
257 1, 4, 5,
258
259 2, 1, 5,
260 2, 5, 6,
261
262 3, 2, 6,
263 3, 6, 7,
264
265 0, 3, 7,
266 0, 7, 4;
267
268 regular_3d_grid(std::max(2., round(1. / pow(area_param_, 1. / 3.) + 1) / 2.), false, false, false, cube_points_, cube_faces_, cube_tets_);
269
270 // Extract sampled edges matching the base element edges
271 Eigen::MatrixXi edges(12, 2);
272 edges << 0, 1,
273 1, 2,
274 2, 3,
275 3, 0,
276
277 4, 5,
278 5, 6,
279 6, 7,
280 7, 4,
281
282 0, 4,
283 1, 5,
284 2, 6,
285 3, 7;
286 igl::edges(cube_faces_, cube_edges_);
287 extract_parent_edges(cube_points_, cube_edges_, pts, edges, cube_edges_);
288
289 // Same local order as in FEMBasis3d
290 cube_corners_.resize(8, 3);
291 cube_corners_ << 0, 0, 0,
292 1, 0, 0,
293 1, 1, 0,
294 0, 1, 0,
295 0, 0, 1,
296 1, 0, 1,
297 1, 1, 1,
298 0, 1, 1;
299 }
300
301 // tet
302 {
303 MatrixXd pts(4, 3);
304 pts << 0, 0, 0,
305 1, 0, 0,
306 0, 1, 0,
307 0, 0, 1;
308
309 Eigen::MatrixXi faces(4, 3);
310 faces << 0, 1, 2,
311
312 3, 1, 0,
313 2, 1, 3,
314 0, 2, 3;
315
316 regular_3d_grid(std::max(2., round(1. / pow(area_param_, 1. / 3.) + 1)), true, false, false, simplex_points_, simplex_faces_, simplex_tets_);
317
318 // Extract sampled edges matching the base element edges
319 Eigen::MatrixXi edges;
320 igl::edges(faces, edges);
321 igl::edges(simplex_faces_, simplex_edges_);
322 extract_parent_edges(simplex_points_, simplex_edges_, pts, edges, simplex_edges_);
323
324 // Same local order as in FEMBasis3d
325 simplex_corners_.resize(4, 3);
326 simplex_corners_ << 0, 0, 0,
327 1, 0, 0,
328 0, 1, 0,
329 0, 0, 1;
330 }
331
332 // prism
333 {
334 MatrixXd pts(6, 3);
335 pts << 0, 0, 0,
336 0, 1, 0,
337 1, 0, 0,
338 0, 0, 1,
339 0, 1, 1,
340 1, 0, 1;
341
342 regular_3d_grid(std::max(2., round(1. / pow(area_param_, 1. / 3.) + 1) / 2.),
343 false, true, false, prism_points_, prism_faces_, prism_tets_);
344
345 // Extract sampled edges matching the base element edges
346 Eigen::MatrixXi edges(9, 2);
347 edges << 0, 1,
348 1, 2,
349 2, 0,
350
351 3, 4,
352 4, 5,
353 5, 3,
354
355 0, 3,
356 1, 4,
357 2, 5;
358
359 igl::edges(prism_faces_, prism_edges_);
360 extract_parent_edges(prism_points_, prism_edges_, pts, edges, prism_edges_);
361
362 // Same local order as in FEMBasis3d
363 prism_corners_.resize(6, 3);
364 prism_corners_ << 0, 0, 0,
365 0, 1, 0,
366 1, 0, 0,
367 0, 0, 1,
368 0, 1, 1,
369 1, 0, 1;
370 }
371
372 // pyramid
373 {
374 MatrixXd pts(5, 3);
375 pts << 0, 0, 0,
376 0, 1, 0,
377 1, 1, 0,
378 0, 1, 0,
379 1, 1, 1;
380
381 regular_3d_grid(std::max(2., round(1. / pow(area_param_, 1. / 3.) + 1) / 2.),
382 false, false, true, pyramid_points_, pyramid_faces_, pyramid_tets_);
383
384 // Extract sampled edges matching the base element edges
385 Eigen::MatrixXi edges(8, 2);
386 edges << 0, 1,
387 1, 2,
388 2, 3,
389 3, 0,
390 0, 4,
391 1, 4,
392 2, 4,
393 3, 4;
394
395 igl::edges(pyramid_faces_, pyramid_edges_);
396 extract_parent_edges(pyramid_points_, pyramid_edges_, pts, edges, pyramid_edges_);
397
398 // Same local order as in FEMBasis3d
399 pyramid_corners_.resize(5, 3);
400 pyramid_corners_ << 0, 0, 0,
401 0, 1, 0,
402 1, 1, 0,
403 0, 1, 0,
404 1, 1, 1;
405 }
406 }
407 else
408 {
409 {
410 MatrixXd pts(4, 2);
411 pts << 0, 0,
412 0, 1,
413 1, 1,
414 1, 0;
415
416 MatrixXi E(4, 2);
417 E << 0, 1,
418 1, 2,
419 2, 3,
420 3, 0;
421
422 regular_2d_grid(std::max(2., round(1. / sqrt(area_param_) + 1)), false, cube_points_, cube_faces_);
423
424 // Extract sampled edges matching the base element edges
425 igl::edges(cube_faces_, cube_edges_);
426 extract_parent_edges(cube_points_, cube_edges_, pts, E, cube_edges_);
427
428 // Same local order as in FEMBasis2d
429 cube_corners_.resize(4, 2);
430 cube_corners_ << 0, 0,
431 1, 0,
432 1, 1,
433 0, 1;
434 }
435 {
436 MatrixXd pts(3, 2);
437 pts << 0, 0,
438 0, 1,
439 1, 0;
440
441 MatrixXi E(3, 2);
442 E << 0, 1,
443 1, 2,
444 2, 0;
445
446 regular_2d_grid(std::max(2., round(1. / sqrt(area_param_) + 1)), true, simplex_points_, simplex_faces_);
447
448 // Extract sampled edges matching the base element edges
449 igl::edges(simplex_faces_, simplex_edges_);
450 extract_parent_edges(simplex_points_, simplex_edges_, pts, E, simplex_edges_);
451
452 // Same local order as in FEMBasis2d
453 simplex_corners_.resize(3, 2);
454 simplex_corners_ << 0, 0,
455 1, 0,
456 0, 1;
457 }
458 }
459 }
460
461 void RefElementSampler::sample_polygon(const Eigen::MatrixXd &poly, Eigen::MatrixXd &pts, Eigen::MatrixXi &faces, Eigen::MatrixXi &edges) const
462 {
463
464#ifdef POLYFEM_WITH_TRIANGLE
465 Eigen::MatrixXi E(poly.rows(), 2);
466 const Eigen::MatrixXd H(0, 2);
467 const std::string flags = "Qzqa" + std::to_string(area_param_ / 10.0);
468
469 for (int i = 0; i < poly.rows(); ++i)
470 E.row(i) << i, (i + 1) % poly.rows();
471
472 igl::triangle::triangulate(poly, E, H, flags, pts, faces);
473#else
474
475 const Eigen::MatrixXi rt = Eigen::MatrixXi::Zero(poly.rows(), 1);
476 faces.resize(0, 0);
477 Eigen::VectorXi I;
478 Eigen::MatrixXd area;
479 igl::predicates::ear_clipping(poly, rt, faces, I);
480
481 igl::doublearea(poly, faces, area);
482
483 const double area_avg = area.array().sum() / poly.rows() / 2;
484
485 const int n_refs = area_avg / area_param_ * 40;
486
487 Eigen::MatrixXi new_faces;
488
489 igl::upsample(poly, faces, pts, new_faces, n_refs);
490
491 faces = new_faces;
492#endif
493 std::vector<int> loop;
494 igl::default_num_threads(1);
495 igl::boundary_loop(faces, loop);
496 igl::default_num_threads(0);
497 edges.resize(loop.size(), 2);
498 for (int i = 0; i < loop.size(); ++i)
499 edges.row(i) << loop[i], loop[(i + 1) % loop.size()];
500 }
501
502 void RefElementSampler::sample_polyhedron(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &f, Eigen::MatrixXd &pts, Eigen::MatrixXi &tets, Eigen::MatrixXi &faces) const
503 {
504 const Eigen::MatrixXd kernel = vertices.colwise().mean();
505
506 polyfem::tertrahedralize_star_shaped_surface(vertices, f, kernel, pts, faces, tets);
507 }
508 } // namespace utils
509} // namespace polyfem
V
int V
Definition AdjointNLProblem.cpp:58
faces
std::vector< Eigen::VectorXi > faces
Definition PressureAssembler.cpp:52
polyfem::utils::RefElementSampler::init
void init(const bool is_volume, const int n_elements, const double target_rel_area)
Definition RefElementSampler.cpp:217
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:461
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:502
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:930
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:33
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:120
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