PolyFEM
Loading...
Searching...
No Matches
Mesh3D.cpp
Go to the documentation of this file.
1#include <polyfem/mesh/mesh3D/Mesh3D.hpp>
2#include <polyfem/mesh/MeshUtils.hpp>
3#include <polyfem/utils/StringUtils.hpp>
4
5#include <polyfem/utils/Logger.hpp>
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 RowVectorNd 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 (1 - t) * v1 + t * v2;
90 }
91
92 const auto &n = edge_nodes_[index.edge];
93 if (n.v1 == index.vertex)
94 return n.nodes.row(i - 1);
95 else
96 return n.nodes.row(n.nodes.rows() - i);
97 }
98
99 RowVectorNd 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 if (is_simplex(index.element))
103 {
104 if (orders_.size() <= 0 || orders_(index.element) == 1 || orders_(index.element) == 2 || face_nodes_.empty() || face_nodes_[index.face].nodes.rows() != tmp)
105 {
106 const auto v1 = point(index.vertex);
107 const auto v2 = point(switch_vertex(index).vertex);
108 const auto v3 = point(switch_vertex(switch_edge(index)).vertex);
109
110 const double b2 = i / (n_new_nodes + 2.0);
111 const double b3 = j / (n_new_nodes + 2.0);
112 const double b1 = 1 - b3 - b2;
113 assert(b3 < 1);
114 assert(b3 > 0);
115
116 return b1 * v1 + b2 * v2 + b3 * v3;
117 }
118
119 const int ii = i - 1;
120 const int jj = j - 1;
121
122 assert(orders_(index.element) == 3 || orders_(index.element) == 4);
123 const auto &n = face_nodes_[index.face];
124 int lindex = jj * n_new_nodes + ii;
125
126 if (orders_(index.element) == 4)
127 {
128 static const std::array<int, 3> remapping = {{0, 2, 1}};
129 if (n.v1 == index.vertex)
130 {
131 if (n.v2 != next_around_face(index).vertex)
132 {
133 lindex = remapping[lindex];
134 assert(n.v3 == next_around_face(index).vertex);
135 }
136 else
137 {
138 assert(n.v2 == next_around_face(index).vertex);
139 }
140 }
141 else if (n.v2 == index.vertex)
142 {
143
144 if (n.v3 != next_around_face(index).vertex)
145 {
146 lindex = remapping[lindex];
147 assert(n.v1 == next_around_face(index).vertex);
148 }
149 else
150 {
151 assert(n.v3 == switch_vertex(index).vertex);
152 }
153
154 lindex = (lindex + 1) % 3;
155 }
156 else if (n.v3 == index.vertex)
157 {
158
159 if (n.v1 != next_around_face(index).vertex)
160 {
161 lindex = remapping[lindex];
162 assert(n.v2 == next_around_face(index).vertex);
163 }
164 else
165 {
166 assert(n.v1 == switch_vertex(index).vertex);
167 }
168
169 lindex = (lindex + 2) % 3;
170 }
171 else
172 {
173 // assert(false);
174 }
175 }
176
177 return n.nodes.row(lindex);
178 }
179 else if (is_cube(index.element))
180 {
181 // supports only blilinear quads
182 assert(orders_.size() <= 0 || orders_(index.element) == 1);
183
184 const auto v1 = point(index.vertex);
185 const auto v2 = point(switch_vertex(index).vertex);
186 const auto v3 = point(switch_vertex(switch_edge(switch_vertex(index))).vertex);
187 const auto v4 = point(switch_vertex(switch_edge(index)).vertex);
188
189 const double b1 = i / (n_new_nodes + 1.0);
190 const double b2 = j / (n_new_nodes + 1.0);
191
192 return v1 * (1 - b1) * (1 - b2) + v2 * b1 * (1 - b2) + v3 * b1 * b2 + v4 * (1 - b1) * b2;
193 }
194
195 assert(false);
196 return RowVectorNd(3, 1);
197 }
198
199 RowVectorNd Mesh3D::cell_node(const Navigation3D::Index &index, const int n_new_nodes, const int i, const int j, const int k) const
200 {
201 if (is_simplex(index.element) && orders_.size() > 0 && orders_(index.element) == n_new_nodes + 3)
202 {
203 assert(n_new_nodes == 1); // test higher than 4 order meshes
204 const auto &n = cell_nodes_[index.element];
205 assert(n.nodes.rows() == 1);
206 return n.nodes;
207 }
208
209 if (n_new_nodes == 1)
210 return cell_barycenter(index.element);
211
212 if (is_simplex(index.element))
213 {
214 if (n_new_nodes == 1)
215 return cell_barycenter(index.element);
216 else
217 {
218 const auto v1 = point(index.vertex);
219 const auto v2 = point(switch_vertex(index).vertex);
220 const auto v3 = point(switch_vertex(switch_edge(switch_vertex(index))).vertex);
221 const auto v4 = point(switch_vertex(switch_edge(switch_face(index))).vertex);
222
223 const double w1 = double(i) / (n_new_nodes + 3);
224 const double w2 = double(j) / (n_new_nodes + 3);
225 const double w3 = double(k) / (n_new_nodes + 3);
226 const double w4 = 1 - w1 - w2 - w3;
227
228 // return v1 * w3
229 // + v2 * w4
230 // + v3 * w1
231 // + v4 * w2;
232
233 return w4 * v1 + w1 * v2 + w2 * v3 + w3 * v4;
234 }
235 }
236 else if (is_cube(index.element))
237 {
238 // supports only blilinear hexes
239 assert(orders_.size() <= 0 || orders_(index.element) == 1);
240
241 const auto v1 = point(index.vertex);
242 const auto v2 = point(switch_vertex(index).vertex);
243 const auto v3 = point(switch_vertex(switch_edge(switch_vertex(index))).vertex);
244 const auto v4 = point(switch_vertex(switch_edge(index)).vertex);
245
246 const Navigation3D::Index index1 = switch_face(switch_edge(switch_vertex(switch_edge(switch_face(index)))));
247 const auto v5 = point(index1.vertex);
248 const auto v6 = point(switch_vertex(index1).vertex);
249 const auto v7 = point(switch_vertex(switch_edge(switch_vertex(index1))).vertex);
250 const auto v8 = point(switch_vertex(switch_edge(index1)).vertex);
251
252 const double b1 = i / (n_new_nodes + 1.0);
253 const double b2 = j / (n_new_nodes + 1.0);
254
255 const double b3 = k / (n_new_nodes + 1.0);
256
257 RowVectorNd blin1 = v1 * (1 - b1) * (1 - b2) + v2 * b1 * (1 - b2) + v3 * b1 * b2 + v4 * (1 - b1) * b2;
258 RowVectorNd blin2 = v5 * (1 - b1) * (1 - b2) + v6 * b1 * (1 - b2) + v7 * b1 * b2 + v8 * (1 - b1) * b2;
259
260 return (1 - b3) * blin1 + b3 * blin2;
261 }
262
263 assert(false);
264 return RowVectorNd(3, 1);
265 }
266
267 void Mesh3D::to_face_functions(std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 6> &to_face) const
268 {
269 // top
270 to_face[0] = [&](Navigation3D::Index idx) { return switch_face(switch_edge(switch_vertex(switch_edge(switch_face(idx))))); };
271 // bottom
272 to_face[1] = [&](Navigation3D::Index idx) { return idx; };
273
274 // left
275 to_face[2] = [&](Navigation3D::Index idx) { return switch_face(switch_edge(switch_vertex(idx))); };
276 // right
277 to_face[3] = [&](Navigation3D::Index idx) { return switch_face(switch_edge(idx)); };
278
279 // back
280 to_face[4] = [&](Navigation3D::Index idx) { return switch_face(switch_edge(switch_vertex(switch_edge(switch_vertex(idx))))); };
281 // front
282 to_face[5] = [&](Navigation3D::Index idx) { return switch_face(idx); };
283 }
284
285 void Mesh3D::to_vertex_functions(std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 8> &to_vertex) const
286 {
287 to_vertex[0] = [&](Navigation3D::Index idx) { return idx; };
288 to_vertex[1] = [&](Navigation3D::Index idx) { return switch_vertex(idx); };
289 to_vertex[2] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(switch_vertex(idx))); };
290 to_vertex[3] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(idx)); };
291
292 to_vertex[4] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(switch_face(idx))); };
293 to_vertex[5] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(switch_vertex(switch_edge(switch_face(idx))))); };
294 to_vertex[6] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(switch_face(switch_vertex(switch_edge(switch_vertex(idx)))))); };
295 to_vertex[7] = [&](Navigation3D::Index idx) { return switch_vertex(switch_edge(switch_face(switch_vertex(switch_edge(idx))))); };
296 }
297
298 void Mesh3D::to_edge_functions(std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 12> &to_edge) const
299 {
300 to_edge[0] = [&](Navigation3D::Index idx) { return idx; };
301 to_edge[1] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(idx)); };
302 to_edge[2] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(switch_edge(switch_vertex(idx)))); };
303 to_edge[3] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(switch_edge(switch_vertex(switch_edge(switch_vertex(idx)))))); };
304
305 to_edge[4] = [&](Navigation3D::Index idx) { return switch_edge(switch_face(idx)); };
306 to_edge[5] = [&](Navigation3D::Index idx) { return switch_edge(switch_face(switch_edge(switch_vertex(idx)))); };
307 to_edge[6] = [&](Navigation3D::Index idx) { return switch_edge(switch_face(switch_edge(switch_vertex(switch_edge(switch_vertex(idx)))))); };
308 to_edge[7] = [&](Navigation3D::Index idx) { return switch_edge(switch_face(switch_edge(switch_vertex(switch_edge(switch_vertex(switch_edge(switch_vertex(idx)))))))); };
309
310 to_edge[8] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(switch_edge(switch_face(idx)))); };
311 to_edge[9] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(switch_edge(switch_face(switch_edge(switch_vertex(idx)))))); };
312 to_edge[10] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(switch_edge(switch_face(switch_edge(switch_vertex(switch_edge(switch_vertex(idx)))))))); };
313 to_edge[11] = [&](Navigation3D::Index idx) { return switch_edge(switch_vertex(switch_edge(switch_face(switch_edge(switch_vertex(switch_edge(switch_vertex(switch_edge(switch_vertex(idx)))))))))); };
314 }
315
316 // v7────v6
317 // ╱┆ ╱│
318 // v4─┼──v5 │
319 // │v3┄┄┄┼v2
320 // │╱ │╱
321 // v0────v1
322 std::array<int, 8> Mesh3D::get_ordered_vertices_from_hex(const int element_index) const
323 {
324 assert(is_cube(element_index));
325 auto idx = get_index_from_element(element_index);
326 std::array<int, 8> v;
327
328 std::array<std::function<Navigation3D::Index(Navigation3D::Index)>, 8> to_vertex;
329 to_vertex_functions(to_vertex);
330 for (int i = 0; i < 8; ++i)
331 v[i] = to_vertex[i](idx).vertex;
332
333 // for (int lv = 0; lv < 4; ++lv) {
334 // v[lv] = idx.vertex;
335 // idx = next_around_face_of_element(idx);
336 // }
337 // // assert(idx == get_index_from_element(element_index));
338 // idx = switch_face(switch_edge(switch_vertex(switch_edge(switch_face(idx)))));
339 // for (int lv = 0; lv < 4; ++lv) {
340 // v[4+lv] = idx.vertex;
341 // idx = next_around_face_of_element(idx);
342 // }
343 return v;
344 }
345
346 std::array<int, 4> Mesh3D::get_ordered_vertices_from_tet(const int element_index) const
347 {
348 auto idx = get_index_from_element(element_index);
349 std::array<int, 4> v;
350
351 for (int lv = 0; lv < 3; ++lv)
352 {
353 v[lv] = idx.vertex;
354 idx = next_around_face(idx);
355 }
356 // assert(idx == get_index_from_element(element_index));
357 idx = switch_vertex(switch_edge(switch_face(idx)));
358 v[3] = idx.vertex;
359
360 // std::array<GEO::vec3, 4> vertices;
361
362 // for(int lv = 0; lv < 4; ++lv)
363 // {
364 // auto pt = point(v[lv]);
365 // for(int d = 0; d < 3; ++d)
366 // {
367 // vertices[lv][d] = pt(d);
368 // }
369 // }
370
371 // const double vol = GEO::Geom::tetra_signed_volume(vertices[0], vertices[1], vertices[2], vertices[3]);
372 // if(vol < 0)
373 // {
374 // std::cout << "negative vol" << std::endl;
375 // // idx = switch_vertex(get_index_from_element(element_index));
376 // // for (int lv = 0; lv < 3; ++lv) {
377 // // v[lv] = idx.vertex;
378 // // idx = next_around_face(idx);
379 // // }
380 // //// assert(idx == get_index_from_element(element_index));
381 // // idx = switch_vertex(switch_edge(switch_face(idx)));
382 // // v[3] = idx.vertex;
383 // }
384
385 return v;
386 }
387
388 void Mesh3D::elements_boxes(std::vector<std::array<Eigen::Vector3d, 2>> &boxes) const
389 {
390 boxes.resize(n_elements());
391
392 for (int i = 0; i < n_elements(); ++i)
393 {
394 auto &box = boxes[i];
395 box[0].setConstant(std::numeric_limits<double>::max());
396 box[1].setConstant(std::numeric_limits<double>::min());
397
398 for (int j = 0; j < n_cell_vertices(i); ++j)
399 {
400 const int v_id = cell_vertex(i, j);
401 for (int d = 0; d < 3; ++d)
402 {
403 box[0][d] = std::min(box[0][d], point(v_id)[d]);
404 box[1][d] = std::max(box[1][d], point(v_id)[d]);
405 }
406 }
407 }
408 }
409
410 void Mesh3D::barycentric_coords(const RowVectorNd &p, const int el_id, Eigen::MatrixXd &coord) const
411 {
412 assert(is_simplex(el_id));
413
414 const auto indices = get_ordered_vertices_from_tet(el_id);
415
416 const auto A = point(indices[0]);
417 const auto B = point(indices[1]);
418 const auto C = point(indices[2]);
419 const auto D = point(indices[3]);
420
421 igl::barycentric_coordinates(p, A, B, C, D, coord);
422 }
423
424 void Mesh3D::compute_cell_jacobian(const int el_id, const Eigen::MatrixXd &reference_map, Eigen::MatrixXd &jacobian) const
425 {
426 assert(is_simplex(el_id));
427
428 const auto indices = get_ordered_vertices_from_tet(el_id);
429
430 const auto A = point(indices[0]);
431 const auto B = point(indices[1]);
432 const auto C = point(indices[2]);
433 const auto D = point(indices[3]);
434
435 Eigen::MatrixXd coords(4, 4);
436 coords << A, 1, B, 1, C, 1, D, 1;
437 coords.transposeInPlace();
438
439 jacobian = coords * reference_map;
440
441 assert(jacobian.determinant() > 0);
442 }
443
444 } // namespace mesh
445} // namespace polyfem
polyfem::mesh::Mesh3D::to_vertex_functions
void to_vertex_functions(std::array< std::function< Navigation3D::Index(Navigation3D::Index)>, 8 > &to_vertex) const
Definition Mesh3D.cpp:285
polyfem::mesh::Mesh3D::get_index_from_element
virtual Navigation3D::Index get_index_from_element(int hi, int lf, int lv) const =0
polyfem::mesh::Mesh3D::tri_area
double tri_area(const int gid) const override
area of a tri face of a tet mesh
Definition Mesh3D.cpp:18
polyfem::mesh::Mesh3D::edge_node
RowVectorNd edge_node(const Navigation3D::Index &index, const int n_new_nodes, const int i) const
Definition Mesh3D.cpp:80
polyfem::mesh::Mesh3D::compute_cell_jacobian
void compute_cell_jacobian(const int el_id, const Eigen::MatrixXd &reference_map, Eigen::MatrixXd &jacobian) const
Definition Mesh3D.cpp:424
polyfem::mesh::Mesh3D::n_cell_edges
virtual int n_cell_edges(const int c_id) const =0
polyfem::mesh::Mesh3D::cell_edge
virtual int cell_edge(const int c_id, const int le_id) const =0
polyfem::mesh::Mesh3D::elements_boxes
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:388
polyfem::mesh::Mesh3D::get_ordered_vertices_from_hex
std::array< int, 8 > get_ordered_vertices_from_hex(const int element_index) const
Definition Mesh3D.cpp:322
polyfem::mesh::Mesh3D::cell_node
RowVectorNd cell_node(const Navigation3D::Index &index, const int n_new_nodes, const int i, const int j, const int k) const
Definition Mesh3D.cpp:199
polyfem::mesh::Mesh3D::get_edges
void get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1) const override
Get all the edges.
Definition Mesh3D.cpp:33
polyfem::mesh::Mesh3D::switch_edge
virtual Navigation3D::Index switch_edge(Navigation3D::Index idx) const =0
polyfem::mesh::Mesh3D::barycentric_coords
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:410
polyfem::mesh::Mesh3D::get_ordered_vertices_from_tet
virtual std::array< int, 4 > get_ordered_vertices_from_tet(const int element_index) const
Definition Mesh3D.cpp:346
polyfem::mesh::Mesh3D::to_edge_functions
void to_edge_functions(std::array< std::function< Navigation3D::Index(Navigation3D::Index)>, 12 > &to_edge) const
Definition Mesh3D.cpp:298
polyfem::mesh::Mesh3D::next_around_face
virtual Navigation3D::Index next_around_face(Navigation3D::Index idx) const =0
polyfem::mesh::Mesh3D::face_node
RowVectorNd face_node(const Navigation3D::Index &index, const int n_new_nodes, const int i, const int j) const
Definition Mesh3D.cpp:99
polyfem::mesh::Mesh3D::switch_face
virtual Navigation3D::Index switch_face(Navigation3D::Index idx) const =0
polyfem::mesh::Mesh3D::to_face_functions
void to_face_functions(std::array< std::function< Navigation3D::Index(Navigation3D::Index)>, 6 > &to_face) const
Definition Mesh3D.cpp:267
polyfem::mesh::Mesh3D::switch_vertex
virtual Navigation3D::Index switch_vertex(Navigation3D::Index idx) const =0
polyfem::mesh::Mesh::n_elements
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
Definition Mesh.hpp:158
polyfem::mesh::Mesh::n_vertices
virtual int n_vertices() const =0
number of vertices
polyfem::mesh::Mesh::orders_
Eigen::MatrixXi orders_
list of geometry orders, one per cell
Definition Mesh.hpp:666
polyfem::mesh::Mesh::is_cube
bool is_cube(const int el_id) const
checks if element is cube compatible
Definition Mesh.cpp:348
polyfem::mesh::Mesh::point
virtual RowVectorNd point(const int global_index) const =0
point coordinates
polyfem::mesh::Mesh::is_simplex
bool is_simplex(const int el_id) const
checks if element is simples compatible
Definition Mesh.cpp:418
polyfem::mesh::Mesh::cell_nodes_
std::vector< CellNodes > cell_nodes_
high-order nodes associates to cells
Definition Mesh.hpp:675
polyfem::mesh::Mesh::edge_vertex
virtual int edge_vertex(const int e_id, const int lv_id) const =0
id of the edge vertex
polyfem::mesh::Mesh::face_nodes_
std::vector< FaceNodes > face_nodes_
high-order nodes associates to faces
Definition Mesh.hpp:673
polyfem::mesh::Mesh::edge_nodes_
std::vector< EdgeNodes > edge_nodes_
high-order nodes associates to edges
Definition Mesh.hpp:671
polyfem::mesh::Mesh::cell_barycenter
virtual RowVectorNd cell_barycenter(const int c) const =0
cell barycenter
polyfem::mesh::Mesh::n_edges
virtual int n_edges() const =0
number of edges
polyfem::mesh::Mesh::cell_vertex
virtual int cell_vertex(const int f_id, const int lv_id) const =0
id of the vertex of a cell
polyfem::mesh::Mesh::n_face_vertices
virtual int n_face_vertices(const int f_id) const =0
number of vertices of a face
polyfem::mesh::Mesh::n_cell_vertices
virtual int n_cell_vertices(const int c_id) const =0
number of vertices of a cell
polyfem::mesh::Mesh::face_vertex
virtual int face_vertex(const int f_id, const int lv_id) const =0
id of the face vertex
polyfem::RowVectorNd
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:11