PolyFEM
Loading...
Searching...
No Matches
Mesh3D.cpp
Go to the documentation of this file.
4
6
7#include <igl/barycentric_coordinates.h>
8
9#include <geogram/mesh/mesh_io.h>
10#include <fstream>
11
12namespace polyfem
13{
14 using namespace utils;
15
16 namespace mesh
17 {
18 double Mesh3D::tri_area(const int gid) const
19 {
20 const int n_vertices = n_face_vertices(gid);
21 assert(n_vertices == 3);
22
23 const auto v1 = point(face_vertex(gid, 0));
24 const auto v2 = point(face_vertex(gid, 1));
25 const auto v3 = point(face_vertex(gid, 2));
26
27 const Eigen::Vector3d e0 = (v2 - v1).transpose();
28 const Eigen::Vector3d e1 = (v3 - v1).transpose();
29
30 return e0.cross(e1).norm() / 2;
31 }
32
33 void Mesh3D::get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1) const
34 {
35 p0.resize(n_edges(), 3);
36 p1.resize(p0.rows(), p0.cols());
37
38 for (std::size_t e = 0; e < n_edges(); ++e)
39 {
40 const int v0 = edge_vertex(e, 0);
41 const int v1 = edge_vertex(e, 1);
42
43 p0.row(e) = point(v0);
44 p1.row(e) = point(v1);
45 }
46 }
47
48 void Mesh3D::get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1, const std::vector<bool> &valid_elements) const
49 {
50 int count = 0;
51 for (size_t i = 0; i < valid_elements.size(); ++i)
52 {
53 if (valid_elements[i])
54 {
55 count += n_cell_edges(i);
56 }
57 }
58
59 p0.resize(count, 3);
60 p1.resize(count, 3);
61
62 count = 0;
63
64 for (size_t i = 0; i < valid_elements.size(); ++i)
65 {
66 if (!valid_elements[i])
67 continue;
68
69 for (size_t ei = 0; ei < n_cell_edges(i); ++ei)
70 {
71 const int e = cell_edge(i, ei);
72 p0.row(count) = point(edge_vertex(e, 0));
73 p1.row(count) = point(edge_vertex(e, 1));
74
75 ++count;
76 }
77 }
78 }
79
80 std::pair<RowVectorNd, int> Mesh3D::edge_node(const Navigation3D::Index &index, const int n_new_nodes, const int i) const
81 {
82 if (orders_.size() <= 0 || orders_(index.element) == 1 || edge_nodes_.empty() || edge_nodes_[index.edge].nodes.rows() != n_new_nodes)
83 {
84 const auto v1 = point(index.vertex);
85 const auto v2 = point(switch_vertex(index).vertex);
86
87 const double t = i / (n_new_nodes + 1.0);
88
89 return std::make_pair((1 - t) * v1 + t * v2, -1);
90 }
91
92 const auto &n = edge_nodes_[index.edge];
93 if (n.v1 == index.vertex)
94 return std::make_pair(n.nodes.row(i - 1), n.nodes_ids[i - 1]);
95 else
96 return std::make_pair(n.nodes.row(n.nodes.rows() - i), n.nodes_ids[n.nodes_ids.size() - i]);
97 }
98
99 std::pair<RowVectorNd, int> Mesh3D::face_node(const Navigation3D::Index &index, const int n_new_nodes, const int i, const int j) const
100 {
101 const int tmp = n_new_nodes == 2 ? 3 : n_new_nodes;
102
103 const bool is_prism_tri = is_prism(index.element) && n_face_vertices(index.face) == 3;
104 const bool is_prism_quad = is_prism(index.element) && n_face_vertices(index.face) == 4;
105
106 const bool is_pyramid_tri = is_pyramid(index.element) && n_face_vertices(index.face) == 3;
107 const bool is_pyramid_quad = is_pyramid(index.element) && n_face_vertices(index.face) == 4;
108
109 if (is_simplex(index.element) || is_prism_tri || is_pyramid_tri)
110 {
111 if (orders_.size() <= 0 || orders_(index.element) == 1 || orders_(index.element) == 2 || face_nodes_.empty() || face_nodes_[index.face].nodes.rows() != tmp)
112 {
113 const auto v1 = point(index.vertex);
114 const auto v2 = point(switch_vertex(index).vertex);
115 const auto v3 = point(switch_vertex(switch_edge(index)).vertex);
116
117 const double b2 = i / (n_new_nodes + 2.0);
118 const double b3 = j / (n_new_nodes + 2.0);
119 const double b1 = 1 - b3 - b2;
120 assert(b3 < 1);
121 assert(b3 > 0);
122
123 return std::make_pair(b1 * v1 + b2 * v2 + b3 * v3, -1);
124 }
125
126 const int ii = i - 1;
127 const int jj = j - 1;
128
129 assert(orders_(index.element) == 3 || orders_(index.element) == 4);
130 const auto &n = face_nodes_[index.face];
131 int lindex = jj * n_new_nodes + ii;
132
133 if (orders_(index.element) == 4)
134 {
135 static const std::array<int, 3> remapping = {{0, 2, 1}};
136 if (n.v1 == index.vertex)
137 {
138 if (n.v2 != next_around_face(index).vertex)
139 {
140 lindex = remapping[lindex];
141 assert(n.v3 == next_around_face(index).vertex);
142 }
143 else
144 {
145 assert(n.v2 == next_around_face(index).vertex);
146 }
147 }
148 else if (n.v2 == index.vertex)
149 {
150
151 if (n.v3 != next_around_face(index).vertex)
152 {
153 lindex = remapping[lindex];
154 assert(n.v1 == next_around_face(index).vertex);
155 }
156 else
157 {
158 assert(n.v3 == switch_vertex(index).vertex);
159 }
160
161 lindex = (lindex + 1) % 3;
162 }
163 else if (n.v3 == index.vertex)
164 {
165
166 if (n.v1 != next_around_face(index).vertex)
167 {
168 lindex = remapping[lindex];
169 assert(n.v2 == next_around_face(index).vertex);
170 }
171 else
172 {
173 assert(n.v1 == switch_vertex(index).vertex);
174 }
175
176 lindex = (lindex + 2) % 3;
177 }
178 else
179 {
180 // assert(false);
181 }
182 }
183
184 return std::make_pair(n.nodes.row(lindex), n.nodes_ids[lindex]);
185 }
186 else if (is_cube(index.element) || is_prism_quad || is_pyramid_quad)
187 {
188 // supports only blilinear quads
189 assert(orders_.size() <= 0 || orders_(index.element) == 1);
190
191 const auto v1 = point(index.vertex);
192 const auto v2 = point(switch_vertex(index).vertex);
193 const auto v3 = point(switch_vertex(switch_edge(switch_vertex(index))).vertex);
194 const auto v4 = point(switch_vertex(switch_edge(index)).vertex);
195
196 const double b1 = i / (n_new_nodes + 1.0);
197 const double b2 = j / (n_new_nodes + 1.0);
198
199 return std::make_pair(v1 * (1 - b1) * (1 - b2) + v2 * b1 * (1 - b2) + v3 * b1 * b2 + v4 * (1 - b1) * b2, -1);
200 }
201
202 assert(false);
203 return std::make_pair(RowVectorNd(3, 1), -1);
204 }
205
206 std::pair<RowVectorNd, int> Mesh3D::cell_node(const Navigation3D::Index &index, const int n_new_nodes, const int i, const int j, const int k) const
207 {
208 if (is_simplex(index.element) && orders_.size() > 0 && orders_(index.element) == n_new_nodes + 3)
209 {
210 assert(n_new_nodes == 1); // test higher than 4 order meshes
211 const auto &n = cell_nodes_[index.element];
212 assert(n.nodes.rows() == 1);
213 return std::make_pair(n.nodes, n.nodes_ids[0]);
214 }
215
216 if (n_new_nodes == 1)
217 return std::make_pair(cell_barycenter(index.element), -1);
218
219 if (is_simplex(index.element))
220 {
221 if (n_new_nodes == 1)
222 return std::make_pair(cell_barycenter(index.element), -1);
223 else
224 {
225 const auto v1 = point(index.vertex);
226 const auto v2 = point(switch_vertex(index).vertex);
227 const auto v3 = point(switch_vertex(switch_edge(switch_vertex(index))).vertex);
228 const auto v4 = point(switch_vertex(switch_edge(switch_face(index))).vertex);
229
230 const double w1 = double(i) / (n_new_nodes + 3);
231 const double w2 = double(j) / (n_new_nodes + 3);
232 const double w3 = double(k) / (n_new_nodes + 3);
233 const double w4 = 1 - w1 - w2 - w3;
234
235 // return v1 * w3
236 // + v2 * w4
237 // + v3 * w1
238 // + v4 * w2;
239
240 return std::make_pair(w4 * v1 + w1 * v2 + w2 * v3 + w3 * v4, -1);
241 }
242 }
243 else if (is_cube(index.element))
244 {
245 // supports only blilinear hexes
246 assert(orders_.size() <= 0 || orders_(index.element) == 1);
247
248 const auto v1 = point(index.vertex);
249 const auto v2 = point(switch_vertex(index).vertex);
250 const auto v3 = point(switch_vertex(switch_edge(switch_vertex(index))).vertex);
251 const auto v4 = point(switch_vertex(switch_edge(index)).vertex);
252
254 const auto v5 = point(index1.vertex);
255 const auto v6 = point(switch_vertex(index1).vertex);
256 const auto v7 = point(switch_vertex(switch_edge(switch_vertex(index1))).vertex);
257 const auto v8 = point(switch_vertex(switch_edge(index1)).vertex);
258
259 const double b1 = i / (n_new_nodes + 1.0);
260 const double b2 = j / (n_new_nodes + 1.0);
261
262 const double b3 = k / (n_new_nodes + 1.0);
263
264 RowVectorNd blin1 = v1 * (1 - b1) * (1 - b2) + v2 * b1 * (1 - b2) + v3 * b1 * b2 + v4 * (1 - b1) * b2;
265 RowVectorNd blin2 = v5 * (1 - b1) * (1 - b2) + v6 * b1 * (1 - b2) + v7 * b1 * b2 + v8 * (1 - b1) * b2;
266
267 return std::make_pair((1 - b3) * blin1 + b3 * blin2, -1);
268 }
269 else if (is_pyramid(index.element))
270 {
271 // n_new_nodes = p-1, k = z-level (1..n-1), i,j = grid indices
272 const double n = n_new_nodes;
273 const double zv = 1.0 - k / n; // z coordinate
274
275 double xv, yv;
276 if (k == 1) // single node at centroid of slice
277 {
278 xv = 0.5 * (1.0 - zv);
279 yv = 0.5 * (1.0 - zv);
280 }
281 else // k×k grid
282 {
283 const double hk = 1.0 / (k + 1);
284 xv = i * hk * (1.0 - zv);
285 yv = j * hk * (1.0 - zv);
286 }
287 return std::make_pair(RowVectorNd{{xv, yv, zv}}, -1);
288 }
289
290 assert(false);
291 return std::make_pair(RowVectorNd(3, 1), -1);
292 }
293
294 void Mesh3D::to_face_functions(std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 6> &to_face) const
295 {
296 // top
297 to_face[0] = [&](Navigation3D::Index idx) { return switch_face(switch_edge(switch_vertex(switch_edge(switch_face(idx))))); };
298 // bottom
299 to_face[1] = [&](Navigation3D::Index idx) { return idx; };
300
301 // left
302 to_face[2] = [&](Navigation3D::Index idx) { return switch_face(switch_edge(switch_vertex(idx))); };
303 // right
304 to_face[3] = [&](Navigation3D::Index idx) { return switch_face(switch_edge(idx)); };
305
306 // back
307 to_face[4] = [&](Navigation3D::Index idx) { return switch_face(switch_edge(switch_vertex(switch_edge(switch_vertex(idx))))); };
308 // front
309 to_face[5] = [&](Navigation3D::Index idx) { return switch_face(idx); };
310 }
311
312 void Mesh3D::to_vertex_functions(std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 8> &to_vertex) const
313 {
314 to_vertex[0] = [&](Navigation3D::Index idx) { return idx; };
315 to_vertex[1] = [&](Navigation3D::Index idx) { return switch_vertex(idx); };
316 to_vertex[2] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(switch_vertex(idx))); };
317 to_vertex[3] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(idx)); };
318
319 to_vertex[4] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(switch_face(idx))); };
320 to_vertex[5] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(switch_vertex(switch_edge(switch_face(idx))))); };
322 to_vertex[7] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(switch_face(switch_vertex(switch_edge(idx))))); };
323 }
324
325 void Mesh3D::to_edge_functions(std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 12> &to_edge) const
326 {
327 to_edge[0] = [&](Navigation3D::Index idx) { return idx; };
328 to_edge[1] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(idx)); };
329 to_edge[2] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(switch_edge(switch_vertex(idx)))); };
331
332 to_edge[4] = [&](Navigation3D::Index idx) { return switch_edge(switch_face(idx)); };
333 to_edge[5] = [&](Navigation3D::Index idx) { return switch_edge(switch_face(switch_edge(switch_vertex(idx)))); };
336
337 to_edge[8] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(switch_edge(switch_face(idx)))); };
341 }
342
343 // v7────v6
344 // ╱┆ ╱│
345 // v4─┼──v5 │
346 // │v3┄┄┄┼v2
347 // │╱ │╱
348 // v0────v1
349 std::array<int, 8> Mesh3D::get_ordered_vertices_from_hex(const int element_index) const
350 {
351 assert(is_cube(element_index));
352 auto idx = get_index_from_element(element_index);
353 std::array<int, 8> v;
354
355 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 8> to_vertex;
356 to_vertex_functions(to_vertex);
357 for (int i = 0; i < 8; ++i)
358 v[i] = to_vertex[i](idx).vertex;
359
360 // for (int lv = 0; lv < 4; ++lv) {
361 // v[lv] = idx.vertex;
362 // idx = next_around_face_of_element(idx);
363 // }
364 // // assert(idx == get_index_from_element(element_index));
365 // idx = switch_face(switch_edge(switch_vertex(switch_edge(switch_face(idx)))));
366 // for (int lv = 0; lv < 4; ++lv) {
367 // v[4+lv] = idx.vertex;
368 // idx = next_around_face_of_element(idx);
369 // }
370 return v;
371 }
372
373 std::array<int, 4> Mesh3D::get_ordered_vertices_from_tet(const int element_index) const
374 {
375 auto idx = get_index_from_element(element_index);
376 std::array<int, 4> v;
377
378 for (int lv = 0; lv < 3; ++lv)
379 {
380 v[lv] = idx.vertex;
381 idx = next_around_face(idx);
382 }
383 // assert(idx == get_index_from_element(element_index));
385 v[3] = idx.vertex;
386
387 return v;
388 }
389
390 // v5
391 // ╱┆ \
392 // v3─┼──v4
393 // │v2 |
394 // │╱ \ │
395 // v0────v1
396 std::array<int, 6> Mesh3D::get_ordered_vertices_from_prism(const int element_index) const
397 {
398 assert(is_prism(element_index));
399 auto idx = get_index_from_element(element_index);
400 std::array<int, 6> v;
401
402 Navigation3D::Index start = idx;
403
404 for (int i = 0; i < 3; ++i)
405 idx = next_around_face(idx);
406
407 if (idx.vertex != start.vertex)
408 start = switch_face(idx);
409
410 for (int i = 0; i < 3; ++i)
411 {
412 v[i] = start.vertex;
413 start = next_around_face(start);
414 }
415 assert(start.vertex == v[0]);
416
418 for (int i = 0; i < 3; ++i)
419 {
420 v[i + 3] = start.vertex;
421 start = next_around_face(start);
422 }
423 assert(start.vertex == v[3]);
424 return v;
425 }
426
427 std::array<int, 5> Mesh3D::get_ordered_vertices_from_pyramid(const int element_index) const
428 {
429 auto idx = get_index_from_element(element_index);
430 std::array<int, 5> v;
431
432 Navigation3D::Index base; // we want to start navigation from base
433 if (n_face_vertices(idx.face) == 4)
434 {
435 base = idx;
436 }
437 else
438 {
439 Navigation3D::Index tmp = idx;
440 bool found = false;
441 for (int k = 0; k < 3; ++k)
442 {
443 auto nb = switch_face(tmp);
444 if (n_face_vertices(nb.face) == 4)
445 {
446 base = nb;
447 found = true;
448 break;
449 }
450 tmp = next_around_face(tmp);
451 }
452 assert(found);
453 }
454
455 idx = base;
456 for (int lv = 0; lv < 4; ++lv)
457 {
458 v[lv] = idx.vertex;
459 idx = next_around_face(idx);
460 }
461 // assert(idx == get_index_from_element(element_index));
463 v[4] = idx.vertex;
464
465 return v;
466 }
467
468 void Mesh3D::elements_boxes(std::vector<std::array<Eigen::Vector3d, 2>> &boxes) const
469 {
470 boxes.resize(n_elements());
471
472 for (int i = 0; i < n_elements(); ++i)
473 {
474 auto &box = boxes[i];
475 box[0].setConstant(std::numeric_limits<double>::max());
476 box[1].setConstant(std::numeric_limits<double>::min());
477
478 for (int j = 0; j < n_cell_vertices(i); ++j)
479 {
480 const int v_id = cell_vertex(i, j);
481 for (int d = 0; d < 3; ++d)
482 {
483 box[0][d] = std::min(box[0][d], point(v_id)[d]);
484 box[1][d] = std::max(box[1][d], point(v_id)[d]);
485 }
486 }
487 }
488 }
489
490 void Mesh3D::barycentric_coords(const RowVectorNd &p, const int el_id, Eigen::MatrixXd &coord) const
491 {
492 assert(is_simplex(el_id));
493
494 const auto indices = get_ordered_vertices_from_tet(el_id);
495
496 const auto A = point(indices[0]);
497 const auto B = point(indices[1]);
498 const auto C = point(indices[2]);
499 const auto D = point(indices[3]);
500
501 igl::barycentric_coordinates(p, A, B, C, D, coord);
502 }
503
504 void Mesh3D::compute_cell_jacobian(const int el_id, const Eigen::MatrixXd &reference_map, Eigen::MatrixXd &jacobian) const
505 {
506 assert(is_simplex(el_id));
507
508 const auto indices = get_ordered_vertices_from_tet(el_id);
509
510 const auto A = point(indices[0]);
511 const auto B = point(indices[1]);
512 const auto C = point(indices[2]);
513 const auto D = point(indices[3]);
514
515 Eigen::MatrixXd coords(4, 4);
516 coords << A, 1, B, 1, C, 1, D, 1;
517 coords.transposeInPlace();
518
519 jacobian = coords * reference_map;
520
521 assert(jacobian.determinant() > 0);
522 }
523
524 } // namespace mesh
525} // namespace polyfem
void to_vertex_functions(std::array< std::function< Navigation3D::Index(Navigation3D::Index)>, 8 > &to_vertex) const
Definition Mesh3D.cpp:312
virtual Navigation3D::Index get_index_from_element(int hi, int lf, int lv) const =0
double tri_area(const int gid) const override
area of a tri face of a tet mesh
Definition Mesh3D.cpp:18
void compute_cell_jacobian(const int el_id, const Eigen::MatrixXd &reference_map, Eigen::MatrixXd &jacobian) const
Definition Mesh3D.cpp:504
virtual int n_cell_edges(const int c_id) const =0
virtual int cell_edge(const int c_id, const int le_id) const =0
std::array< int, 5 > get_ordered_vertices_from_pyramid(const int element_index) const
Definition Mesh3D.cpp:427
void elements_boxes(std::vector< std::array< Eigen::Vector3d, 2 > > &boxes) const override
constructs a box around every element (3d cell, 2d face)
Definition Mesh3D.cpp:468
std::array< int, 8 > get_ordered_vertices_from_hex(const int element_index) const
Definition Mesh3D.cpp:349
std::pair< RowVectorNd, int > cell_node(const Navigation3D::Index &index, const int n_new_nodes, const int i, const int j, const int k) const
Definition Mesh3D.cpp:206
void get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1) const override
Get all the edges.
Definition Mesh3D.cpp:33
virtual Navigation3D::Index switch_edge(Navigation3D::Index idx) const =0
void barycentric_coords(const RowVectorNd &p, const int el_id, Eigen::MatrixXd &coord) const override
constructs barycentric coodiantes for a point p.
Definition Mesh3D.cpp:490
std::pair< RowVectorNd, int > face_node(const Navigation3D::Index &index, const int n_new_nodes, const int i, const int j) const
Definition Mesh3D.cpp:99
std::pair< RowVectorNd, int > edge_node(const Navigation3D::Index &index, const int n_new_nodes, const int i) const
Definition Mesh3D.cpp:80
virtual std::array< int, 4 > get_ordered_vertices_from_tet(const int element_index) const
Definition Mesh3D.cpp:373
void to_edge_functions(std::array< std::function< Navigation3D::Index(Navigation3D::Index)>, 12 > &to_edge) const
Definition Mesh3D.cpp:325
std::array< int, 6 > get_ordered_vertices_from_prism(const int element_index) const
Definition Mesh3D.cpp:396
virtual Navigation3D::Index next_around_face(Navigation3D::Index idx) const =0
virtual Navigation3D::Index switch_face(Navigation3D::Index idx) const =0
void to_face_functions(std::array< std::function< Navigation3D::Index(Navigation3D::Index)>, 6 > &to_face) const
Definition Mesh3D.cpp:294
virtual Navigation3D::Index switch_vertex(Navigation3D::Index idx) const =0
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
Definition Mesh.hpp:163
virtual int n_vertices() const =0
number of vertices
Eigen::MatrixXi orders_
list of geometry orders, one per cell
Definition Mesh.hpp:705
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
bool is_simplex(const int el_id) const
checks if element is simplex
Definition Mesh.cpp:422
bool is_prism(const int el_id) const
checks if element is a prism
Definition Mesh.cpp:427
std::vector< CellNodes > cell_nodes_
high-order nodes associates to cells
Definition Mesh.hpp:714
virtual int edge_vertex(const int e_id, const int lv_id) const =0
id of the edge vertex
std::vector< FaceNodes > face_nodes_
high-order nodes associates to faces
Definition Mesh.hpp:712
std::vector< EdgeNodes > edge_nodes_
high-order nodes associates to edges
Definition Mesh.hpp:710
bool is_pyramid(const int el_id) const
checks if element is a pyramid
Definition Mesh.cpp:432
virtual RowVectorNd cell_barycenter(const int c) const =0
cell barycenter
virtual int n_edges() const =0
number of edges
virtual int cell_vertex(const int f_id, const int lv_id) const =0
id of the vertex of a cell
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
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13