PolyFEM
Loading...
Searching...
No Matches
MeshNodes.cpp
Go to the documentation of this file.
1
2#include <polyfem/mesh/MeshNodes.hpp>
3#include <polyfem/mesh/mesh2D/CMesh2D.hpp>
4#include <polyfem/mesh/mesh2D/NCMesh2D.hpp>
5#include <polyfem/mesh/mesh3D/CMesh3D.hpp>
6#include <polyfem/mesh/mesh3D/NCMesh3D.hpp>
7
8#include <polyfem/utils/Logger.hpp>
10
11namespace polyfem::mesh
12{
13 namespace
14 {
15
16 std::vector<bool> interface_edges(const Mesh2D &mesh)
17 {
18 std::vector<bool> is_interface(mesh.n_edges(), false);
19 for (int f = 0; f < mesh.n_faces(); ++f)
20 {
21 if (mesh.is_cube(f))
22 {
23 continue;
24 } // Skip quads
25
26 auto index = mesh.get_index_from_face(f, 0);
27 for (int lv = 0; lv < mesh.n_face_vertices(f); ++lv)
28 {
29 auto index2 = mesh.switch_face(index);
30 if (index2.face >= 0)
31 {
32 is_interface[index.edge] = true;
33 }
34 index = mesh.next_around_face(index);
35 }
36 }
37 return is_interface;
38 }
39
40 std::vector<bool> interface_faces(const Mesh3D &mesh)
41 {
42 std::vector<bool> is_interface(mesh.n_faces(), false);
43 for (int c = 0; c < mesh.n_cells(); ++c)
44 {
45 if (mesh.is_cube(c))
46 {
47 continue;
48 } // Skip hexes
49
50 for (int lf = 0; lf < mesh.n_cell_faces(c); ++lf)
51 {
52 auto index = mesh.get_index_from_element(c, lf, 0);
53 auto index2 = mesh.switch_element(index);
54 if (index2.element >= 0)
55 {
56 is_interface[index.face] = true;
57 }
58 }
59 }
60 return is_interface;
61 }
62
63 } // anonymous namespace
64
65 // -----------------------------------------------------------------------------
66
67 MeshNodes::MeshNodes(const Mesh &mesh, const bool has_poly, const bool connect_nodes, const int max_nodes_per_edge, const int max_nodes_per_face, const int max_nodes_per_cell)
68 : mesh_(mesh), connect_nodes_(connect_nodes), edge_offset_(mesh.n_vertices()), face_offset_(edge_offset_ + mesh.n_edges() * max_nodes_per_edge), cell_offset_(face_offset_ + mesh.n_faces() * max_nodes_per_face)
69
70 ,
71 max_nodes_per_edge_(max_nodes_per_edge), max_nodes_per_face_(max_nodes_per_face), max_nodes_per_cell_(max_nodes_per_cell)
72 {
73 // Initialization
74 int n_nodes = cell_offset_ + mesh.n_cells() * max_nodes_per_cell;
75 primitive_to_node_.assign(n_nodes, -1); // #v + #e + #f + #c
76 nodes_.resize(n_nodes, mesh.dimension());
77 is_boundary_.assign(n_nodes, false);
78
79 if (has_poly)
80 is_interface_.assign(n_nodes, false);
81
82 // Vertex nodes
83 for (int v = 0; v < mesh.n_vertices(); ++v)
84 {
85 // nodes_.row(v) = mesh.point(v);
86 is_boundary_[v] = mesh.is_boundary_vertex(v);
87 }
88 // if (!vertices_only) {
89 // Eigen::MatrixXd bary;
90 // Edge nodes
91 // mesh.edge_barycenters(bary);
92
93 for (int e = 0; e < mesh.n_edges(); ++e)
94 {
95 // nodes_.row(edge_offset_ + e) = bary.row(e);
96 const bool is_boundary = mesh.is_boundary_edge(e);
97 for (int tmp = 0; tmp < max_nodes_per_edge; ++tmp)
98 is_boundary_[edge_offset_ + max_nodes_per_edge * e + tmp] = is_boundary;
99 }
100 // Face nodes
101 // mesh.face_barycenters(bary);
102 for (int f = 0; f < mesh.n_faces(); ++f)
103 {
104 // nodes_.row(face_offset_ + f) = bary.row(f);
105 for (int tmp = 0; tmp < max_nodes_per_face; ++tmp)
106 {
107 if (mesh.is_volume())
108 {
109 is_boundary_[face_offset_ + max_nodes_per_face * f + tmp] = mesh.is_boundary_face(f);
110 }
111 else
112 {
113 is_boundary_[face_offset_ + max_nodes_per_face * f + tmp] = false;
114 }
115 }
116 }
117 // Cell nodes
118 // mesh.cell_barycenters(bary);
119 for (int c = 0; c < mesh.n_cells(); ++c)
120 {
121 for (int tmp = 0; tmp < max_nodes_per_cell; ++tmp)
122 is_boundary_[cell_offset_ + max_nodes_per_cell * c + tmp] = false;
123 }
124 // }
125
126 // Vertices only, no need to compute interface marker
127 // if (vertices_only) { return; }
128
129 // Compute edges/faces that are at interface with a polytope
130 if (has_poly)
131 {
132 if (mesh.is_volume())
133 {
134 const Mesh3D *mesh3d = dynamic_cast<const Mesh3D *>(&mesh);
135 assert(mesh3d);
136 auto is_interface_face = interface_faces(*mesh3d);
137 for (int f = 0; f < mesh.n_faces(); ++f)
138 {
139 for (int tmp = 0; tmp < max_nodes_per_face; ++tmp)
140 is_interface_[face_offset_ + max_nodes_per_face * f + tmp] = is_interface_face[f];
141 }
142 }
143 else
144 {
145 const Mesh2D *mesh2d = dynamic_cast<const Mesh2D *>(&mesh);
146 assert(mesh2d);
147 auto is_interface_edge = interface_edges(*mesh2d);
148 for (int e = 0; e < mesh.n_edges(); ++e)
149 {
150 for (int tmp = 0; tmp < max_nodes_per_edge; ++tmp)
151 is_interface_[edge_offset_ + max_nodes_per_edge * e + tmp] = is_interface_edge[e];
152 }
153 }
154 }
155 }
156
158
159 int MeshNodes::node_id_from_primitive(int primitive_id)
160 {
161 if (primitive_to_node_[primitive_id] < 0 || !connect_nodes_)
162 {
163 primitive_to_node_[primitive_id] = n_nodes();
164
165 in_ordered_vertices_.push_back(primitive_id);
166 node_to_primitive_.push_back(primitive_id);
167 node_to_primitive_gid_.push_back(primitive_id);
168 assert(in_ordered_vertices_.size() == n_nodes());
169
170 if (primitive_id < edge_offset_)
171 nodes_.row(primitive_id) = mesh_.point(primitive_id);
172 else if (primitive_id < face_offset_)
173 nodes_.row(primitive_id) = mesh_.edge_barycenter(primitive_id - edge_offset_);
174 else if (primitive_id < cell_offset_)
175 nodes_.row(primitive_id) = mesh_.face_barycenter(primitive_id - face_offset_);
176 else
177 nodes_.row(primitive_id) = mesh_.cell_barycenter(primitive_id - cell_offset_);
178 }
179 return primitive_to_node_[primitive_id];
180 }
181
182 std::vector<int> MeshNodes::node_ids_from_edge(const Navigation::Index &index, const int n_new_nodes)
183 {
184 std::vector<int> res;
185 if (n_new_nodes <= 0)
186 return res;
187
188 const Mesh2D *mesh2d = dynamic_cast<const Mesh2D *>(&mesh_);
189 int start;
190
191 if (connect_nodes_)
192 start = edge_offset_ + index.edge * max_nodes_per_edge_;
193 else
194 {
195 if (mesh2d->is_boundary_edge(index.edge) || mesh2d->switch_face(index).face > index.face)
196 start = edge_offset_ + index.edge * max_nodes_per_edge_;
197 else
198 start = edge_offset_ + index.edge * max_nodes_per_edge_ + max_nodes_per_edge_ / 2;
199 }
200
201 // const auto v1 = mesh2d->point(index.vertex);
202 // const auto v2 = mesh2d->point(mesh2d->switch_vertex(index).vertex);
203 assert(start < primitive_to_node_.size());
204 const int start_node_id = primitive_to_node_[start];
205#ifndef NDEBUG
206 if (!connect_nodes_)
207 {
208 assert(start_node_id < 0);
209 }
210#endif
211 if (start_node_id < 0 || !connect_nodes_)
212 {
213 for (int i = 1; i <= n_new_nodes; ++i)
214 {
215 // const double t = i/(n_new_nodes + 1.0);
216
217 const int primitive_id = start + i - 1;
218 assert(primitive_id < primitive_to_node_.size());
219
220 const auto [node, node_id] = mesh2d->edge_node(index, n_new_nodes, i);
221
222 primitive_to_node_[primitive_id] = n_nodes();
223
224 in_ordered_vertices_.push_back(node_id);
225 node_to_primitive_.push_back(primitive_id);
226 node_to_primitive_gid_.push_back(index.edge);
227
228 nodes_.row(primitive_id) = node;
229
230 res.push_back(primitive_to_node_[primitive_id]);
231 assert(in_ordered_vertices_.size() == n_nodes());
232 }
233 }
234 else
235 {
236 const auto [v, _] = mesh2d->edge_node(index, n_new_nodes, 1);
237
238 if ((node_position(start_node_id) - v).norm() < 1e-10)
239 {
240 for (int i = 0; i < n_new_nodes; ++i)
241 {
242 const int primitive_id = start + i;
243 res.push_back(primitive_to_node_[primitive_id]);
244 }
245 }
246 else
247 {
248 for (int i = n_new_nodes - 1; i >= 0; --i)
249 {
250 const int primitive_id = start + i;
251 res.push_back(primitive_to_node_[primitive_id]);
252 }
253 }
254 }
255
256 assert(res.size() == size_t(n_new_nodes));
257 return res;
258 }
259
260 std::vector<int> MeshNodes::node_ids_from_edge(const Navigation3D::Index &index, const int n_new_nodes)
261 {
262 std::vector<int> res;
263 if (n_new_nodes <= 0)
264 return res;
265
266 const int start = edge_offset_ + index.edge * max_nodes_per_edge_;
267
268 const Mesh3D *mesh3d = dynamic_cast<const Mesh3D *>(&mesh_);
269
270 const int start_node_id = primitive_to_node_[start];
271 if (start_node_id < 0)
272 {
273 for (int i = 1; i <= n_new_nodes; ++i)
274 {
275 const auto [node, node_id] = mesh3d->edge_node(index, n_new_nodes, i);
276
277 const int primitive_id = start + i - 1;
278 primitive_to_node_[primitive_id] = n_nodes();
279
280 in_ordered_vertices_.push_back(node_id);
281 node_to_primitive_.push_back(primitive_id);
282 node_to_primitive_gid_.push_back(index.edge);
283
284 nodes_.row(primitive_id) = node;
285
286 res.push_back(primitive_to_node_[primitive_id]);
287 assert(in_ordered_vertices_.size() == n_nodes());
288 }
289 }
290 else
291 {
292 const auto [v, _] = mesh3d->edge_node(index, n_new_nodes, 1);
293 if ((node_position(start_node_id) - v).norm() < 1e-10)
294 {
295 for (int i = 0; i < n_new_nodes; ++i)
296 {
297 const int primitive_id = start + i;
298 res.push_back(primitive_to_node_[primitive_id]);
299 }
300 }
301 else
302 {
303 for (int i = n_new_nodes - 1; i >= 0; --i)
304 {
305 const int primitive_id = start + i;
306 res.push_back(primitive_to_node_[primitive_id]);
307 }
308 }
309 }
310
311 assert(res.size() == size_t(n_new_nodes));
312 return res;
313 }
314
315 std::vector<int> MeshNodes::node_ids_from_face(const Navigation::Index &index, const int n_new_nodes)
316 {
317 std::vector<int> res;
318 if (n_new_nodes <= 0)
319 return res;
320
321 const int start = face_offset_ + index.face * max_nodes_per_face_;
322 const int start_node_id = primitive_to_node_[start];
323
324 const Mesh2D *mesh2d = dynamic_cast<const Mesh2D *>(&mesh_);
325
326 if (start_node_id < 0 || !connect_nodes_)
327 {
328 int loc_index = 0;
329 for (int i = 1; i <= n_new_nodes; ++i)
330 {
331 const int end = mesh2d->is_simplex(index.face) ? (n_new_nodes - i + 1) : n_new_nodes;
332 for (int j = 1; j <= end; ++j)
333 {
334 const auto [node, node_id] = mesh2d->face_node(index, n_new_nodes, i, j);
335
336 const int primitive_id = start + loc_index;
337 primitive_to_node_[primitive_id] = n_nodes();
338
339 in_ordered_vertices_.push_back(node_id);
340 node_to_primitive_.push_back(primitive_id);
341 node_to_primitive_gid_.push_back(index.face);
342
343 nodes_.row(primitive_id) = node;
344
345 res.push_back(primitive_to_node_[primitive_id]);
346
347 assert(in_ordered_vertices_.size() == n_nodes());
348
349 ++loc_index;
350 }
351 }
352 }
353 else
354 {
355 assert(false);
356 }
357
358#ifndef NDEBUG
359 if (mesh2d->is_simplex(index.face))
360 assert(res.size() == size_t(n_new_nodes * (n_new_nodes + 1) / 2));
361 else
362 assert(res.size() == size_t(n_new_nodes * n_new_nodes));
363#endif
364 return res;
365 }
366
367 std::vector<int> MeshNodes::node_ids_from_face(const Navigation3D::Index &index1, const int n_new_nodes, const int n_new_nodesq)
368 {
369 Navigation3D::Index index = index1;
370
371 std::vector<int> res;
372 if (n_new_nodes <= 0)
373 return res;
374
375 const Mesh3D *mesh3d = dynamic_cast<const Mesh3D *>(&mesh_);
376 assert(mesh3d->is_prism(index.element));
377
378 int start;
379 const bool is_tri = mesh3d->n_face_vertices(index.face) == 3;
380
381 if (connect_nodes_)
382 start = face_offset_ + index.face * max_nodes_per_face_;
383 else
384 {
385 if (mesh3d->is_boundary_face(index.face) || mesh3d->switch_element(index).element > index.element)
386 start = face_offset_ + index.face * max_nodes_per_face_;
387 else
388 start = face_offset_ + index.face * max_nodes_per_face_ + max_nodes_per_face_ / 2;
389 }
390
391 const int start_node_id = primitive_to_node_[start];
392
393#ifndef NDEBUG
394 if (!connect_nodes_)
395 {
396 assert(start_node_id < 0);
397 }
398#endif
399
400 if (start_node_id < 0 || !connect_nodes_)
401 {
402 int loc_index = 0;
403
404 if (is_tri)
405 {
406 for (int i = 1; i <= n_new_nodes; ++i)
407 {
408 const int end = n_new_nodes - i + 1;
409 for (int j = 1; j <= end; ++j)
410 {
411 const int primitive_id = start + loc_index;
412 const auto [node, node_id] = mesh3d->face_node(index, n_new_nodes, -1, i, j);
413
414 primitive_to_node_[primitive_id] = n_nodes();
415
416 in_ordered_vertices_.push_back(node_id);
417 node_to_primitive_.push_back(primitive_id);
418 node_to_primitive_gid_.push_back(index.face);
419
420 nodes_.row(primitive_id) = node;
421
422 res.push_back(primitive_to_node_[primitive_id]);
423 assert(in_ordered_vertices_.size() == n_nodes());
424
425 ++loc_index;
426 }
427 }
428 }
429 else
430 {
431 const bool is_on_tri_face = mesh3d->n_face_vertices(mesh3d->switch_face(index).face) == 3; // on quad face, next to tri face
432 if (is_on_tri_face)
433 index = mesh3d->next_around_face(index);
434 assert(mesh3d->n_face_vertices(mesh3d->switch_face(index).face) == 4); // on quad face, next to quad face
435
436 for (int i = 1; i <= n_new_nodes; ++i)
437 {
438 for (int j = 1; j <= n_new_nodesq; ++j)
439 {
440 const int primitive_id = start + loc_index;
441 const auto [node, node_id] = mesh3d->face_node(index, n_new_nodes, n_new_nodesq, j, i);
442
443 primitive_to_node_[primitive_id] = n_nodes();
444
445 in_ordered_vertices_.push_back(node_id);
446 node_to_primitive_.push_back(primitive_id);
447 node_to_primitive_gid_.push_back(index.face);
448
449 nodes_.row(primitive_id) = node;
450
451 res.push_back(primitive_to_node_[primitive_id]);
452 assert(in_ordered_vertices_.size() == n_nodes());
453
454 ++loc_index;
455 }
456 }
457 }
458 }
459 else // find existing nodes, todo
460 {
461 if (is_tri)
462 {
463 const int total_nodes = n_new_nodes * (n_new_nodes + 1) / 2;
464 for (int i = 1; i <= n_new_nodes; ++i)
465 {
466 const int end = n_new_nodes - i + 1;
467 for (int j = 1; j <= end; ++j)
468 {
469 const auto [p, _] = mesh3d->face_node(index, n_new_nodes, n_new_nodesq, i, j);
470
471 bool found = false;
472 for (int k = start; k < start + total_nodes; ++k)
473 {
474 const double dist = (nodes_.row(k) - p).norm();
475 if (dist < 1e-10)
476 {
477 res.push_back(primitive_to_node_[k]);
478 found = true;
479 break;
480 }
481 }
482
483 assert(found);
484 }
485 }
486 }
487 else
488 {
489 const bool is_on_tri_face = mesh3d->n_face_vertices(mesh3d->switch_face(index).face) == 3; // on quad face, next to tri face
490 if (is_on_tri_face)
491 index = mesh3d->next_around_face(index);
492 assert(mesh3d->n_face_vertices(mesh3d->switch_face(index).face) == 4); // on quad face, next to quad face
493
494 const int total_nodes = n_new_nodes * n_new_nodesq;
495 for (int i = 1; i <= n_new_nodes; ++i)
496 {
497 const int end = n_new_nodesq;
498 for (int j = 1; j <= end; ++j)
499 {
500 const auto [p, _] = mesh3d->face_node(index, n_new_nodes, n_new_nodesq, j, i);
501
502 bool found = false;
503 for (int k = start; k < start + total_nodes; ++k)
504 {
505 const double dist = (nodes_.row(k) - p).norm();
506 if (dist < 1e-10)
507 {
508 res.push_back(primitive_to_node_[k]);
509 found = true;
510 break;
511 }
512 }
513
514 assert(found);
515 }
516 }
517 }
518 }
519
520#ifndef NDEBUG
521 if (is_tri)
522 assert(res.size() == size_t(n_new_nodes * (n_new_nodes + 1) / 2));
523 else
524 assert(res.size() == size_t(n_new_nodes * n_new_nodesq));
525#endif
526 return res;
527 }
528
529 std::vector<int> MeshNodes::node_ids_from_face(const Navigation3D::Index &index, const int n_new_nodes)
530 {
531 std::vector<int> res;
532 if (n_new_nodes <= 0)
533 return res;
534
535 const Mesh3D *mesh3d = dynamic_cast<const Mesh3D *>(&mesh_);
536 assert(!mesh3d->is_prism(index.element));
537
538 int start;
539 const bool is_tri = mesh3d->is_simplex(index.element) || (mesh3d->is_prism(index.element) && mesh3d->n_face_vertices(index.face) == 3) || (mesh3d->is_pyramid(index.element) && mesh3d->n_face_vertices(index.face) == 3);
540
541 if (connect_nodes_)
542 start = face_offset_ + index.face * max_nodes_per_face_;
543 else
544 {
545 if (mesh3d->is_boundary_face(index.face) || mesh3d->switch_element(index).element > index.element)
546 start = face_offset_ + index.face * max_nodes_per_face_;
547 else
548 start = face_offset_ + index.face * max_nodes_per_face_ + max_nodes_per_face_ / 2;
549 }
550
551 const int start_node_id = primitive_to_node_[start];
552
553#ifndef NDEBUG
554 if (!connect_nodes_)
555 {
556 assert(start_node_id < 0);
557 }
558#endif
559
560 if (start_node_id < 0 || !connect_nodes_)
561 {
562 int loc_index = 0;
563 for (int i = 1; i <= n_new_nodes; ++i)
564 {
565 const int end = is_tri ? (n_new_nodes - i + 1) : n_new_nodes;
566 for (int j = 1; j <= end; ++j)
567 {
568 const int primitive_id = start + loc_index;
569 const auto [node, node_id] = mesh3d->face_node(index, n_new_nodes, i, j);
570
571 primitive_to_node_[primitive_id] = n_nodes();
572
573 in_ordered_vertices_.push_back(node_id);
574 node_to_primitive_.push_back(primitive_id);
575 node_to_primitive_gid_.push_back(index.face);
576
577 nodes_.row(primitive_id) = node;
578
579 res.push_back(primitive_to_node_[primitive_id]);
580 assert(in_ordered_vertices_.size() == n_nodes());
581
582 ++loc_index;
583 }
584 }
585 }
586 else
587 {
588 if (n_new_nodes == 1)
589 {
590 res.push_back(start_node_id);
591 }
592 else
593 {
594 const int total_nodes = is_tri ? (n_new_nodes * (n_new_nodes + 1) / 2) : (n_new_nodes * n_new_nodes);
595 for (int i = 1; i <= n_new_nodes; ++i)
596 {
597 const int end = is_tri ? (n_new_nodes - i + 1) : n_new_nodes;
598 for (int j = 1; j <= end; ++j)
599 {
600 const auto [p, _] = mesh3d->face_node(index, n_new_nodes, i, j);
601
602 bool found = false;
603 for (int k = start; k < start + total_nodes; ++k)
604 {
605 const double dist = (nodes_.row(k) - p).norm();
606 if (dist < 1e-10)
607 {
608 res.push_back(primitive_to_node_[k]);
609 found = true;
610 break;
611 }
612 }
613
614 assert(found);
615 }
616 }
617 }
618 }
619
620#ifndef NDEBUG
621 if (is_tri)
622 assert(res.size() == size_t(n_new_nodes * (n_new_nodes + 1) / 2));
623 else
624 assert(res.size() == size_t(n_new_nodes * n_new_nodes));
625#endif
626 return res;
627 }
628
629 std::vector<int> MeshNodes::node_ids_from_cell(const Navigation3D::Index &index, const int n_new_nodes)
630 {
631 std::vector<int> res;
632 const int start = cell_offset_ + index.element * max_nodes_per_cell_;
633 const Mesh3D *mesh3d = dynamic_cast<const Mesh3D *>(&mesh_);
634
635 if (mesh3d->is_simplex(index.element))
636 {
637 int loc_index = 0;
638 for (int i = 1; i <= n_new_nodes; ++i)
639 {
640 const int endj = (n_new_nodes - i + 1);
641 for (int j = 1; j <= endj; ++j)
642 {
643 const int endk = (n_new_nodes - i - j + 2);
644 for (int k = 1; k <= endk; ++k)
645 {
646 const int primitive_id = start + loc_index;
647
648 auto [node, node_id] = mesh3d->cell_node(index, n_new_nodes, i, j, k);
649
650 primitive_to_node_[primitive_id] = n_nodes();
651
652 in_ordered_vertices_.push_back(node_id);
653 node_to_primitive_.push_back(primitive_id);
654 node_to_primitive_gid_.push_back(index.element);
655
656 nodes_.row(primitive_id) = node;
657 res.push_back(primitive_to_node_[primitive_id]);
658 assert(in_ordered_vertices_.size() == n_nodes());
659
660 ++loc_index;
661 }
662 }
663 }
664 }
665 else if (mesh3d->is_cube(index.element))
666 {
667 int loc_index = 0;
668 for (int i = 1; i <= n_new_nodes; ++i)
669 {
670 for (int j = 1; j <= n_new_nodes; ++j)
671 {
672 for (int k = 1; k <= n_new_nodes; ++k)
673 {
674 const int primitive_id = start + loc_index;
675
676 const auto [node, node_id] = mesh3d->cell_node(index, n_new_nodes, i, j, k);
677 primitive_to_node_[primitive_id] = n_nodes();
678
679 in_ordered_vertices_.push_back(node_id);
680 node_to_primitive_.push_back(primitive_id);
681 node_to_primitive_gid_.push_back(index.element);
682
683 nodes_.row(primitive_id) = node;
684 res.push_back(primitive_to_node_[primitive_id]);
685 assert(in_ordered_vertices_.size() == n_nodes());
686
687 ++loc_index;
688 }
689 }
690 }
691 }
692 else if (mesh3d->is_prism(index.element))
693 {
694 assert(n_new_nodes == 1);
695 int loc_index = 0;
696
697 const int primitive_id = start + loc_index;
698
699 const auto [node, node_id] = mesh3d->cell_node(index, n_new_nodes, 0, 0, 0);
700 primitive_to_node_[primitive_id] = n_nodes();
701
702 in_ordered_vertices_.push_back(node_id);
703 node_to_primitive_.push_back(primitive_id);
704 node_to_primitive_gid_.push_back(index.element);
705
706 nodes_.row(primitive_id) = node;
707 res.push_back(primitive_to_node_[primitive_id]);
708 assert(in_ordered_vertices_.size() == n_nodes());
709 }
710 else if (mesh3d->is_pyramid(index.element))
711 {
712 // n_new_nodes = p-1, mirrors create_point_set interior loop in pyramid_bases.py
713 int loc_index = 0;
714 for (int k = 1; k < n_new_nodes; ++k)
715 {
716 const int grid = (k == 1) ? 1 : k; // 1 node for k=1, k×k for k>1
717 for (int i = 1; i <= grid; ++i)
718 {
719 for (int j = 1; j <= grid; ++j)
720 {
721 const int primitive_id = start + loc_index;
722 const auto [node, node_id] = mesh3d->cell_node(index, n_new_nodes, i, j, k);
723
724 primitive_to_node_[primitive_id] = n_nodes();
725 in_ordered_vertices_.push_back(node_id);
726 node_to_primitive_.push_back(primitive_id);
727 node_to_primitive_gid_.push_back(index.element);
728 nodes_.row(primitive_id) = node;
729
730 res.push_back(primitive_to_node_[primitive_id]);
731 assert(in_ordered_vertices_.size() == n_nodes());
732 ++loc_index;
733 }
734 }
735 }
736 }
737
738#ifndef NDEBUG
739 if (res.size() == 1 && connect_nodes_)
740 {
741 const int idx = node_id_from_cell(index.element);
742 assert(idx == res.front());
743 }
744#endif
745
746#ifndef NDEBUG
747 if (mesh3d->is_simplex(index.element))
748 {
749 int n_cell_nodes = 0;
750 for (int pp = 0; pp <= n_new_nodes; ++pp)
751 n_cell_nodes += (pp * (pp + 1) / 2);
752 assert(res.size() == size_t(n_cell_nodes));
753 }
754 else if (mesh3d->is_cube(index.element))
755 assert(res.size() == size_t(n_new_nodes * n_new_nodes * n_new_nodes));
756#endif
757
758 return res;
759 }
760
761 int MeshNodes::node_id_from_vertex(int v)
762 {
763 return node_id_from_primitive(v);
764 }
765
770
775
780
782
783 int MeshNodes::vertex_from_node_id(int node_id) const
784 {
785 const int res = node_to_primitive_[node_id];
786 assert(res >= 0 && res < edge_offset_);
787 return res;
788 }
789
790 int MeshNodes::edge_from_node_id(int node_id) const
791 {
792 const int res = node_to_primitive_[node_id];
793 assert(res >= edge_offset_ && res < face_offset_);
794 return res - edge_offset_;
795 }
796
797 int MeshNodes::face_from_node_id(int node_id) const
798 {
799 const int res = node_to_primitive_[node_id];
800 assert(res >= face_offset_ && res < cell_offset_);
801 return res - face_offset_;
802 }
803
804 int MeshNodes::cell_from_node_id(int node_id) const
805 {
806 const int res = node_to_primitive_[node_id];
807 assert(res >= cell_offset_ && res < n_nodes());
808 return res - cell_offset_;
809 }
810
812
813 std::vector<int> MeshNodes::boundary_nodes() const
814 {
815 std::vector<int> res;
816 res.reserve(n_nodes());
817 for (size_t prim_id = 0; prim_id < primitive_to_node_.size(); ++prim_id)
818 {
819 int node_id = primitive_to_node_[prim_id];
820 if (node_id >= 0 && is_boundary_[prim_id])
821 {
822 res.push_back(node_id);
823 }
824 }
825 return res;
826 }
827
828 int MeshNodes::count_nonnegative_nodes(int start_i, int end_i) const
829 {
830 int count = 0;
831 for (int i = start_i; i < end_i; i++)
832 count += primitive_to_node_[i] >= 0;
833 return count;
834 }
835
836} // namespace polyfem::mesh
polyfem::mesh::Mesh2D::edge_node
virtual std::pair< RowVectorNd, int > edge_node(const Navigation::Index &index, const int n_new_nodes, const int i) const =0
polyfem::mesh::Mesh2D::face_node
virtual std::pair< RowVectorNd, int > face_node(const Navigation::Index &index, const int n_new_nodes, const int i, const int j) const =0
polyfem::mesh::Mesh2D::switch_face
virtual Navigation::Index switch_face(Navigation::Index idx) const =0
polyfem::mesh::Mesh3D::switch_element
virtual Navigation3D::Index switch_element(Navigation3D::Index idx) const =0
polyfem::mesh::Mesh3D::cell_node
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:340
polyfem::mesh::Mesh3D::face_node
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:234
polyfem::mesh::Mesh3D::edge_node
std::pair< RowVectorNd, int > edge_node(const Navigation3D::Index &index, const int n_new_nodes, const int i) const
Definition Mesh3D.cpp:80
polyfem::mesh::Mesh3D::next_around_face
virtual Navigation3D::Index next_around_face(Navigation3D::Index idx) const =0
polyfem::mesh::Mesh3D::switch_face
virtual Navigation3D::Index switch_face(Navigation3D::Index idx) const =0
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:49
polyfem::mesh::Mesh::n_vertices
virtual int n_vertices() const =0
number of vertices
polyfem::mesh::Mesh::edge_barycenter
virtual RowVectorNd edge_barycenter(const int e) const =0
edge barycenter
polyfem::mesh::Mesh::face_barycenter
virtual RowVectorNd face_barycenter(const int f) const =0
face barycenter
polyfem::mesh::Mesh::is_cube
bool is_cube(const int el_id) const
checks if element is cube compatible
Definition Mesh.cpp:437
polyfem::mesh::Mesh::point
virtual RowVectorNd point(const int global_index) const =0
point coordinates
polyfem::mesh::Mesh::is_boundary_face
virtual bool is_boundary_face(const int face_global_id) const =0
is face boundary
polyfem::mesh::Mesh::is_boundary_vertex
virtual bool is_boundary_vertex(const int vertex_global_id) const =0
is vertex boundary
polyfem::mesh::Mesh::is_simplex
bool is_simplex(const int el_id) const
checks if element is simplex
Definition Mesh.cpp:507
polyfem::mesh::Mesh::is_prism
bool is_prism(const int el_id) const
checks if element is a prism
Definition Mesh.cpp:512
polyfem::mesh::Mesh::is_volume
virtual bool is_volume() const =0
checks if mesh is volume
polyfem::mesh::Mesh::is_boundary_edge
virtual bool is_boundary_edge(const int edge_global_id) const =0
is edge boundary
polyfem::mesh::Mesh::dimension
int dimension() const
utily for dimension
Definition Mesh.hpp:164
polyfem::mesh::Mesh::n_cells
virtual int n_cells() const =0
number of cells
polyfem::mesh::Mesh::n_faces
virtual int n_faces() const =0
number of faces
polyfem::mesh::Mesh::is_pyramid
bool is_pyramid(const int el_id) const
checks if element is a pyramid
Definition Mesh.cpp:517
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::n_face_vertices
virtual int n_face_vertices(const int f_id) const =0
number of vertices of a face
polyfem::mesh::MeshNodes::is_boundary
bool is_boundary(int node_id) const
Definition MeshNodes.hpp:71
polyfem::mesh::MeshNodes::node_id_from_primitive
int node_id_from_primitive(int primitive_id)
Definition MeshNodes.cpp:159
polyfem::mesh::MeshNodes::primitive_to_node_
std::vector< int > primitive_to_node_
Definition MeshNodes.hpp:99
polyfem::mesh::MeshNodes::node_ids_from_face
std::vector< int > node_ids_from_face(const Navigation::Index &index, const int n_new_nodes)
Definition MeshNodes.cpp:315
polyfem::mesh::MeshNodes::node_ids_from_cell
std::vector< int > node_ids_from_cell(const Navigation3D::Index &index, const int n_new_nodes)
Definition MeshNodes.cpp:629
polyfem::mesh::MeshNodes::cell_from_node_id
int cell_from_node_id(int node_id) const
Definition MeshNodes.cpp:804
polyfem::mesh::MeshNodes::node_position
RowVectorNd node_position(int node_id) const
Definition MeshNodes.hpp:68
polyfem::mesh::MeshNodes::MeshNodes
MeshNodes(const Mesh &mesh, const bool has_poly, const bool connect_nodes, const int max_nodes_per_edge, const int max_nodes_per_face, const int max_nodes_per_cell=0)
Definition MeshNodes.cpp:67
polyfem::mesh::MeshNodes::node_to_primitive_gid_
std::vector< int > node_to_primitive_gid_
Definition MeshNodes.hpp:101
polyfem::mesh::MeshNodes::nodes_
Eigen::MatrixXd nodes_
Definition MeshNodes.hpp:104
polyfem::mesh::MeshNodes::in_ordered_vertices_
std::vector< int > in_ordered_vertices_
Definition MeshNodes.hpp:109
polyfem::mesh::MeshNodes::face_from_node_id
int face_from_node_id(int node_id) const
Definition MeshNodes.cpp:797
polyfem::mesh::MeshNodes::vertex_from_node_id
int vertex_from_node_id(int node_id) const
Definition MeshNodes.cpp:783
polyfem::mesh::MeshNodes::count_nonnegative_nodes
int count_nonnegative_nodes(int start_i, int end_i) const
Definition MeshNodes.cpp:828
polyfem::mesh::MeshNodes::boundary_nodes
std::vector< int > boundary_nodes() const
Definition MeshNodes.cpp:813
polyfem::mesh::MeshNodes::edge_from_node_id
int edge_from_node_id(int node_id) const
Definition MeshNodes.cpp:790
polyfem::mesh::MeshNodes::node_to_primitive_
std::vector< int > node_to_primitive_
Definition MeshNodes.hpp:100
polyfem::mesh::MeshNodes::node_ids_from_edge
std::vector< int > node_ids_from_edge(const Navigation::Index &index, const int n_new_nodes)
Definition MeshNodes.cpp:182
polyfem::mesh::MeshNodes::is_boundary_
std::vector< bool > is_boundary_
Definition MeshNodes.hpp:105
polyfem::mesh::MeshNodes::is_interface_
std::vector< bool > is_interface_
Definition MeshNodes.hpp:106
elasticity_rhs.f
f
Definition elasticity_rhs.py:173