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 std::cout << "Pyramid sampler points: " << pyramid_points_.rows() << std::endl;
384
385 // Extract sampled edges matching the base element edges
386 Eigen::MatrixXi edges(8, 2);
387 edges << 0, 1,
388 1, 2,
389 2, 3,
390 3, 0,
391 0, 4,
392 1, 4,
393 2, 4,
394 3, 4;
395
396 igl::edges(pyramid_faces_, pyramid_edges_);
397 extract_parent_edges(pyramid_points_, pyramid_edges_, pts, edges, pyramid_edges_);
398
399 // Same local order as in FEMBasis3d
400 pyramid_corners_.resize(5, 3);
401 pyramid_corners_ << 0, 0, 0,
402 0, 1, 0,
403 1, 1, 0,
404 0, 1, 0,
405 1, 1, 1;
406 }
407 }
408 else
409 {
410 {
411 MatrixXd pts(4, 2);
412 pts << 0, 0,
413 0, 1,
414 1, 1,
415 1, 0;
416
417 MatrixXi E(4, 2);
418 E << 0, 1,
419 1, 2,
420 2, 3,
421 3, 0;
422
423 regular_2d_grid(std::max(2., round(1. / sqrt(area_param_) + 1)), false, cube_points_, cube_faces_);
424
425 // Extract sampled edges matching the base element edges
426 igl::edges(cube_faces_, cube_edges_);
427 extract_parent_edges(cube_points_, cube_edges_, pts, E, cube_edges_);
428
429 // Same local order as in FEMBasis2d
430 cube_corners_.resize(4, 2);
431 cube_corners_ << 0, 0,
432 1, 0,
433 1, 1,
434 0, 1;
435 }
436 {
437 MatrixXd pts(3, 2);
438 pts << 0, 0,
439 0, 1,
440 1, 0;
441
442 MatrixXi E(3, 2);
443 E << 0, 1,
444 1, 2,
445 2, 0;
446
447 regular_2d_grid(std::max(2., round(1. / sqrt(area_param_) + 1)), true, simplex_points_, simplex_faces_);
448
449 // Extract sampled edges matching the base element edges
450 igl::edges(simplex_faces_, simplex_edges_);
451 extract_parent_edges(simplex_points_, simplex_edges_, pts, E, simplex_edges_);
452
453 // Same local order as in FEMBasis2d
454 simplex_corners_.resize(3, 2);
455 simplex_corners_ << 0, 0,
456 1, 0,
457 0, 1;
458 }
459 }
460 }
461
462 void RefElementSampler::sample_polygon(const Eigen::MatrixXd &poly, Eigen::MatrixXd &pts, Eigen::MatrixXi &faces, Eigen::MatrixXi &edges) const
463 {
464
465#ifdef POLYFEM_WITH_TRIANGLE
466 Eigen::MatrixXi E(poly.rows(), 2);
467 const Eigen::MatrixXd H(0, 2);
468 const std::string flags = "Qzqa" + std::to_string(area_param_ / 10.0);
469
470 for (int i = 0; i < poly.rows(); ++i)
471 E.row(i) << i, (i + 1) % poly.rows();
472
473 igl::triangle::triangulate(poly, E, H, flags, pts, faces);
474#else
475
476 const Eigen::MatrixXi rt = Eigen::MatrixXi::Zero(poly.rows(), 1);
477 faces.resize(0, 0);
478 Eigen::VectorXi I;
479 Eigen::MatrixXd area;
480 igl::predicates::ear_clipping(poly, rt, faces, I);
481
482 igl::doublearea(poly, faces, area);
483
484 const double area_avg = area.array().sum() / poly.rows() / 2;
485
486 const int n_refs = area_avg / area_param_ * 40;
487
488 Eigen::MatrixXi new_faces;
489
490 igl::upsample(poly, faces, pts, new_faces, n_refs);
491
492 faces = new_faces;
493#endif
494 std::vector<int> loop;
495 igl::default_num_threads(1);
496 igl::boundary_loop(faces, loop);
497 igl::default_num_threads(0);
498 edges.resize(loop.size(), 2);
499 for (int i = 0; i < loop.size(); ++i)
500 edges.row(i) << loop[i], loop[(i + 1) % loop.size()];
501 }
502
503 void RefElementSampler::sample_polyhedron(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &f, Eigen::MatrixXd &pts, Eigen::MatrixXi &tets, Eigen::MatrixXi &faces) const
504 {
505 const Eigen::MatrixXd kernel = vertices.colwise().mean();
506
507 polyfem::tertrahedralize_star_shaped_surface(vertices, f, kernel, pts, faces, tets);
508 }
509 } // namespace utils
510} // 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:462
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:503
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