PolyFEM
Loading...
Searching...
No Matches
OutData.cpp
Go to the documentation of this file.
1#include "OutData.hpp"
2
3#include <array>
4#include <map>
5
6#include "Evaluator.hpp"
7#include "MatrixIO.hpp"
8
9#include <polyfem/assembler/ElementAssemblyValues.hpp>
10#include <polyfem/assembler/AssemblyValues.hpp>
11#include <polyfem/assembler/MatParams.hpp>
12#include <polyfem/assembler/Mass.hpp>
13
14#include <polyfem/basis/ElementBases.hpp>
15
16#include <polyfem/mesh/Obstacle.hpp>
17#include <polyfem/mesh/MeshUtils.hpp>
18#include <polyfem/mesh/mesh2D/Mesh2D.hpp>
19#include <polyfem/mesh/mesh3D/Mesh3D.hpp>
20
21#include <polyfem/utils/getRSS.h>
22#include <polyfem/utils/EdgeSampler.hpp>
23#include <polyfem/utils/Logger.hpp>
24#include <polyfem/utils/par_for.hpp>
25#include <polyfem/utils/BoundarySampler.hpp>
26#include <polyfem/utils/Timer.hpp>
27#include <polyfem/utils/MaybeParallelFor.hpp>
28
29#include <polyfem/autogen/auto_p_bases.hpp>
30#include <polyfem/autogen/auto_q_bases.hpp>
31#include <polyfem/autogen/prism_bases.hpp>
32#include <polyfem/autogen/auto_pyramid_bases.hpp>
33
34#include <paraviewo/VTMWriter.hpp>
35#include <paraviewo/PVDWriter.hpp>
36
37#include <ipc/potentials/normal_adhesion_potential.hpp>
38#include <ipc/potentials/tangential_adhesion_potential.hpp>
39
40#include <SimpleBVH/BVH.hpp>
41
42#include <igl/write_triangle_mesh.h>
43#include <igl/edges.h>
44#include <igl/facet_adjacency_matrix.h>
45#include <igl/connected_components.h>
46
47#include <ipc/ipc.hpp>
48
49#include <algorithm>
50#include <cmath>
51
52#include <filesystem>
53
54namespace polyfem::io
55{
56 bool OutputFieldOptions::export_field(const std::string &field) const
57 {
58 return fields.empty() || std::find(fields.begin(), fields.end(), field) != fields.end();
59 }
60
61 using CellType = paraviewo::CellType;
62 using CellElement = paraviewo::CellElement;
63
64 namespace
65 {
66 void add_output_fields(
67 paraviewo::ParaviewWriter &writer,
68 const OutputSample &sample,
69 const OutputFieldFunction &output_fields)
70 {
71 if (!output_fields)
72 return;
73
74 for (const OutputField &field : output_fields(sample))
75 {
76 if (field.values.rows() <= 0)
77 continue;
78
79 const int expected_rows =
80 field.association == OutputField::Association::Cell
81 ? sample.cell_count
82 : sample.points.rows();
83 if (field.values.rows() != expected_rows)
84 {
85 logger().warn(
86 "Skipping output field '{}' with {} rows; expected {} {} rows",
87 field.name, field.values.rows(), expected_rows,
88 field.association == OutputField::Association::Cell ? "cell" : "point");
89 continue;
90 }
91
92 if (field.association == OutputField::Association::Cell)
93 writer.add_cell_field(field.name, field.values);
94 else
95 writer.add_field(field.name, field.values);
96 }
97 }
98
99 void avoid_pyramid_apex(Eigen::MatrixXd &points)
100 {
101 assert(points.cols() == 3);
102 constexpr double eps = 1e-8;
103 for (int i = 0; i < points.rows(); ++i)
104 {
105 if (std::abs(points(i, 2) - 1.0) < eps)
106 points(i, 2) = 1.0 - eps;
107 }
108 }
109
110 void pyramid_nodes_for_output(const int order, Eigen::MatrixXd &points)
111 {
112 autogen::pyramid_nodes_3d(order, points);
113 avoid_pyramid_apex(points);
114 }
115
116 // ------------------------------------------------------------------
117 // helpers for the hybrid (prism/pyramid) collision-proxy paths
118 // ------------------------------------------------------------------
119
120 // face/edge parameters are fractions k/n with n <= 8, so coordinates
121 // in units of 1/PROXY_PARAM_SCALE are exact integers (lcm of 1..8)
122 constexpr long PROXY_PARAM_SCALE = 840;
123
124 int prism_q_order(const basis::ElementBases &b)
125 {
126 const int p = b.bases.empty() ? 1 : b.bases.front().order();
127 const int n_tri = (p + 1) * (p + 2) / 2;
128 const int q = (n_tri > 0 && b.bases.size() % n_tri == 0) ? int(b.bases.size()) / n_tri - 1 : p;
129 return std::max(1, q);
130 }
131
132 // reference-node layout of an element (vertices first); false if unsupported
133 bool element_ref_nodes(const mesh::Mesh &mesh, const int el, const basis::ElementBases &b, Eigen::MatrixXd &ref_nodes)
134 {
135 const int p = b.bases.empty() ? 1 : b.bases.front().order();
136 if (mesh.is_simplex(el))
137 autogen::p_nodes_3d(p, ref_nodes);
138 else if (mesh.is_cube(el))
139 autogen::q_nodes_3d(p, ref_nodes);
140 else if (mesh.is_prism(el))
141 autogen::prism_nodes_3d(p, prism_q_order(b), ref_nodes);
142 else if (mesh.is_pyramid(el))
143 autogen::pyramid_nodes_3d(p, ref_nodes);
144 else
145 return false;
146 return ref_nodes.rows() == long(b.bases.size());
147 }
148
149 int element_n_vertices(const mesh::Mesh &mesh, const int el)
150 {
151 if (mesh.is_simplex(el))
152 return 4;
153 if (mesh.is_pyramid(el))
154 return 5;
155 if (mesh.is_prism(el))
156 return 6;
157 assert(mesh.is_cube(el));
158 return 8;
159 }
160
161 // local vertex pairs forming the element's edges, derived from the
162 // reference vertex coordinates (no vertex-ordering assumptions)
163 std::vector<std::pair<int, int>> element_ref_edges(const mesh::Mesh &mesh, const int el, const Eigen::MatrixXd &ref_nodes)
164 {
165 const int nv = element_n_vertices(mesh, el);
166 const auto n_coord_diffs = [&](int a, int b) {
167 int d = 0;
168 for (int c = 0; c < 3; ++c)
169 if (std::abs(ref_nodes(a, c) - ref_nodes(b, c)) > 1e-12)
170 ++d;
171 return d;
172 };
173 int apex = -1;
174 if (mesh.is_pyramid(el))
175 for (int a = 0; a < nv; ++a)
176 if (std::abs(ref_nodes(a, 2) - 1.0) < 1e-12)
177 apex = a;
178
179 std::vector<std::pair<int, int>> edges;
180 for (int a = 0; a < nv; ++a)
181 {
182 for (int b = a + 1; b < nv; ++b)
183 {
184 bool is_edge;
185 if (mesh.is_simplex(el)) // every vertex pair
186 is_edge = true;
187 else if (mesh.is_prism(el)) // triangle x [0,1]
188 is_edge = std::abs(ref_nodes(a, 2) - ref_nodes(b, 2)) < 1e-12
189 || (std::abs(ref_nodes(a, 0) - ref_nodes(b, 0)) < 1e-12
190 && std::abs(ref_nodes(a, 1) - ref_nodes(b, 1)) < 1e-12);
191 else if (mesh.is_pyramid(el)) // base square + base-to-apex
192 is_edge = (a == apex || b == apex) || n_coord_diffs(a, b) == 1;
193 else // hex
194 is_edge = n_coord_diffs(a, b) == 1;
195 if (is_edge)
196 edges.emplace_back(a, b);
197 }
198 }
199 return edges;
200 }
201
202 // For every mesh edge, keyed by its endpoint DOFs (smaller first): the
203 // owned edge-interior DOFs located on it as (parameter from the
204 // smaller-DOF endpoint in units of 1/PROXY_PARAM_SCALE, dof, node
205 // position), sorted by parameter. Built from ALL elements: an edge on
206 // the domain boundary can get its DOFs from an interior element (e.g. a
207 // prism wedged between two promoted tets). Stitched nodes
208 // (glob.size() != 1, no DOF of their own) are not registered.
209 using EdgeDofs = std::vector<std::tuple<long, int, Eigen::Vector3d>>;
210 std::map<std::pair<int, int>, EdgeDofs> build_edge_dofs(const mesh::Mesh &mesh, const std::vector<basis::ElementBases> &bases)
211 {
212 std::map<std::pair<int, int>, EdgeDofs> edge_dofs;
213 Eigen::MatrixXd ref_nodes;
214 for (int el = 0; el < int(bases.size()); ++el)
215 {
216 const basis::ElementBases &b = bases[el];
217 if (b.bases.empty() || !element_ref_nodes(mesh, el, b, ref_nodes))
218 continue;
219
220 const int nv = element_n_vertices(mesh, el);
221 std::vector<int> vd(nv, -1);
222 bool ok = true;
223 for (int i = 0; i < nv; ++i)
224 {
225 const auto &glob = b.bases[i].global();
226 assert(glob.size() == 1); // vertices always own their DOF
227 if (glob.size() != 1)
228 {
229 ok = false;
230 break;
231 }
232 vd[i] = glob.front().index;
233 }
234 if (!ok)
235 continue;
236
237 const auto edges = element_ref_edges(mesh, el, ref_nodes);
238 for (int j = nv; j < int(b.bases.size()); ++j)
239 {
240 const auto &glob = b.bases[j].global();
241 if (glob.size() != 1)
242 continue;
243 const Eigen::RowVector3d r = ref_nodes.row(j);
244 for (const auto &e : edges)
245 {
246 const Eigen::RowVector3d pa = ref_nodes.row(e.first);
247 const Eigen::RowVector3d d = ref_nodes.row(e.second) - pa;
248 const double t = (r - pa).dot(d) / d.squaredNorm();
249 if (t < 1e-9 || t > 1 - 1e-9 || ((r - pa) - t * d).norm() > 1e-9)
250 continue; // not on this edge
251 long tl = std::lround(t * PROXY_PARAM_SCALE);
252 assert(std::abs(t * PROXY_PARAM_SCALE - tl) < 1e-6);
253 int va = vd[e.first], vb = vd[e.second];
254 if (va > vb)
255 {
256 std::swap(va, vb);
257 tl = PROXY_PARAM_SCALE - tl;
258 }
259 EdgeDofs &ed = edge_dofs[{va, vb}];
260 const int dof = glob.front().index;
261 bool present = false;
262 for (const auto &existing : ed)
263 if (std::get<1>(existing) == dof)
264 {
265 assert(std::get<0>(existing) == tl);
266 present = true;
267 break;
268 }
269 if (!present)
270 ed.emplace_back(tl, dof, glob.front().node.transpose());
271 break;
272 }
273 }
274 }
275 for (auto &kv : edge_dofs)
276 std::sort(kv.second.begin(), kv.second.end(),
277 [](const EdgeDofs::value_type &x, const EdgeDofs::value_type &y) { return std::get<0>(x) < std::get<0>(y); });
278 return edge_dofs;
279 }
280
281 // Structured triangulation for faces whose points form a complete
282 // regular lattice (all faces away from mixed-order interfaces). The
283 // row-wise patterns keep vertex valence balanced -- the lexicographic
284 // sweep below concentrates fan triangles on lex-extreme corners, which
285 // at order 3+ can push one-ring sizes past the smooth-contact
286 // N_VERT_NEIGHBORS_3D cap. Returns false if the points are not a
287 // complete lattice (caller falls back to the sweep).
288 bool triangulate_lattice(const std::vector<std::array<long, 2>> &pts, const int nfv, std::vector<std::array<int, 3>> &tris)
289 {
290 constexpr long S = PROXY_PARAM_SCALE;
291 const int n = int(pts.size());
292 std::map<std::array<long, 2>, int> id;
293 for (int i = 0; i < n; ++i)
294 if (!id.emplace(pts[i], i).second)
295 return false;
296 tris.clear();
297
298 if (nfv == 3)
299 {
300 // n == (k+1)(k+2)/2 for a complete triangular lattice of order k
301 const int k = int(std::lround((std::sqrt(8.0 * n + 1.0) - 3.0) / 2.0));
302 if ((k + 1) * (k + 2) / 2 != n || k < 1 || S % k != 0)
303 return false;
304 const long h = S / k;
305 const auto at = [&](int i, int j) -> int {
306 const auto it = id.find({{i * h, j * h}});
307 return it == id.end() ? -1 : it->second;
308 };
309 for (int j = 0; j <= k; ++j)
310 for (int i = 0; i <= k - j; ++i)
311 if (at(i, j) < 0)
312 return false;
313 for (int j = 0; j < k; ++j)
314 {
315 for (int i = 0; i < k - j; ++i)
316 {
317 tris.push_back({{at(i, j), at(i + 1, j), at(i, j + 1)}});
318 if (i + j < k - 1)
319 tris.push_back({{at(i + 1, j), at(i + 1, j + 1), at(i, j + 1)}});
320 }
321 }
322 return true;
323 }
324
325 // quad: complete (k1+1) x (k2+1) tensor grid, possibly anisotropic
326 std::set<long> us, vs;
327 for (const auto &p : pts)
328 {
329 us.insert(p[0]);
330 vs.insert(p[1]);
331 }
332 const int k1 = int(us.size()) - 1, k2 = int(vs.size()) - 1;
333 if (k1 < 1 || k2 < 1 || (k1 + 1) * (k2 + 1) != n || S % k1 != 0 || S % k2 != 0)
334 return false;
335 const long h1 = S / k1, h2 = S / k2;
336 const auto at = [&](int i, int j) -> int {
337 const auto it = id.find({{i * h1, j * h2}});
338 return it == id.end() ? -1 : it->second;
339 };
340 for (int j = 0; j <= k2; ++j)
341 for (int i = 0; i <= k1; ++i)
342 if (at(i, j) < 0)
343 return false;
344 for (int j = 0; j < k2; ++j)
345 {
346 for (int i = 0; i < k1; ++i)
347 {
348 tris.push_back({{at(i, j), at(i + 1, j), at(i, j + 1)}});
349 tris.push_back({{at(i + 1, j + 1), at(i, j + 1), at(i + 1, j)}});
350 }
351 }
352 return true;
353 }
354
355 // Triangulate a 2D point set in convex position (boundary points lie on
356 // the convex hull, possibly collinear; interior points allowed) into
357 // strictly positive-area CCW triangles via an incremental lexicographic
358 // sweep. Collinear boundary chains are preserved -- no triangle edge
359 // ever spans a boundary point -- which is what keeps neighboring faces
360 // of the proxy conforming. Coordinates must be exact integers.
361 void triangulate_convex_pointset(const std::vector<std::array<long, 2>> &pts, std::vector<std::array<int, 3>> &tris)
362 {
363 const int n = int(pts.size());
364 tris.clear();
365 std::vector<int> order(n);
366 for (int i = 0; i < n; ++i)
367 order[i] = i;
368 std::sort(order.begin(), order.end(), [&](int a, int b) { return pts[a] < pts[b]; });
369 const auto orient = [&](int a, int b, int c) -> long long {
370 return (long long)(pts[b][0] - pts[a][0]) * (pts[c][1] - pts[a][1])
371 - (long long)(pts[b][1] - pts[a][1]) * (pts[c][0] - pts[a][0]);
372 };
373
374 // initial (possibly collinear) chain
375 std::vector<int> hull;
376 int k = 0;
377 for (; k < n; ++k)
378 {
379 if (hull.size() >= 2 && orient(hull[0], hull[1], order[k]) != 0)
380 break;
381 hull.push_back(order[k]);
382 }
383 if (k == n)
384 return; // all points collinear: nothing to emit
385
386 // attach the first off-line point: fan over the chain, make hull CCW
387 {
388 const int p = order[k];
389 const bool left = orient(hull[0], hull[1], p) > 0;
390 for (int i = 0; i + 1 < int(hull.size()); ++i)
391 {
392 if (left)
393 tris.push_back({{hull[i], hull[i + 1], p}});
394 else
395 tris.push_back({{hull[i + 1], hull[i], p}});
396 }
397 if (!left)
398 std::reverse(hull.begin(), hull.end());
399 hull.push_back(p);
400 ++k;
401 }
402
403 for (; k < n; ++k)
404 {
405 const int p = order[k];
406 const int m = int(hull.size());
407 // visible hull edges: p strictly right of a->b on the CCW cycle
408 std::vector<bool> vis(m);
409 for (int i = 0; i < m; ++i)
410 vis[i] = orient(hull[i], hull[(i + 1) % m], p) < 0;
411 int s = 0;
412 while (s < m && !(vis[s] && !vis[(s + m - 1) % m]))
413 ++s;
414 assert(s < m); // p is outside the hull (lexicographic order)
415 int e = s;
416 while (vis[e % m])
417 {
418 // orient(a, b, p) < 0 => (a, p, b) is CCW
419 tris.push_back({{hull[e % m], p, hull[(e + 1) % m]}});
420 ++e;
421 }
422 // replace the strictly-visible part of the hull with p
423 std::vector<int> new_hull;
424 new_hull.push_back(p);
425 for (int i = e % m; i != s; i = (i + 1) % m)
426 new_hull.push_back(hull[i]);
427 new_hull.push_back(hull[s]);
428 hull.swap(new_hull);
429 }
430 }
431
432 } // namespace
433
434 void OutGeometryData::extract_boundary_mesh_sampled(
435 const mesh::Mesh &mesh,
436 const int n_bases,
437 const std::vector<basis::ElementBases> &bases,
438 const std::vector<mesh::LocalBoundary> &total_local_boundary,
439 Eigen::MatrixXd &node_positions,
440 Eigen::MatrixXi &boundary_edges,
441 Eigen::MatrixXi &boundary_triangles,
442 std::vector<Eigen::Triplet<double>> &displacement_map_entries,
443 const int sampling_order)
444 {
445 using namespace polyfem::mesh;
446
447 // Sample every boundary face on a uniform lattice of the globally
448 // maximal order M through the element bases: the proxy is conforming
449 // and watertight regardless of order mismatches across interfaces
450 // (every shared edge is split into M chords on both sides), stitched
451 // interface nodes are handled intrinsically by the basis evaluation,
452 // and the displacement map carries each proxy vertex's (possibly
453 // weighted) dependence on the real DOFs.
454 if (!mesh.is_volume() || mesh.has_poly()
455 || !dynamic_cast<const Mesh3D &>(mesh).is_conforming())
456 {
457 logger().warn("max_order collision-proxy sampling requires a conforming volume mesh without polytopes; falling back to the standard boundary extraction");
458 extract_boundary_mesh(mesh, n_bases, bases, total_local_boundary,
459 node_positions, boundary_edges, boundary_triangles, displacement_map_entries);
460 return;
461 }
462
463 displacement_map_entries.clear();
464 const Mesh3D &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
465
466 int M = 1;
467 if (sampling_order > 0)
468 M = sampling_order; // explicit lattice order (finer OR coarser than the elements)
469 else
470 for (const LocalBoundary &lb : total_local_boundary)
471 {
472 const basis::ElementBases &b = bases[lb.element_id()];
473 if (b.bases.empty())
474 continue;
475 int o = b.bases.front().order();
476 if (mesh.is_prism(lb.element_id()))
477 o = std::max(o, prism_q_order(b));
478 M = std::max(M, o);
479 }
480
481 // samples shared between neighboring faces merge combinatorially, by
482 // the mesh primitive they lie on -- {0, vertex_id} for face corners,
483 // {1, vmin, vmax, step-from-vmin} for edge samples (every face uses
484 // the same global lattice order M, so steps coincide exactly), and
485 // {2, face_id, i, r} for face-interior samples (a boundary face
486 // belongs to exactly one element) -- no floating-point tolerance
487 std::map<std::array<int, 4>, int> vertex_id;
488 std::vector<Eigen::Vector3d> vertices;
489 std::vector<std::tuple<int, int, int>> proxy_tris;
490
491 for (const LocalBoundary &lb : total_local_boundary)
492 {
493 const int el = lb.element_id();
494 const basis::ElementBases &b = bases[el];
495 Eigen::MatrixXd ref_nodes;
496 if (b.bases.empty() || !element_ref_nodes(mesh, el, b, ref_nodes))
497 continue;
498
499 // the rational pyramid bases are 0/0 at the apex (z=1), which is a
500 // corner of every triangular pyramid face; apex samples are instead
501 // evaluated below via the nodal limit (apex basis 1, all others 0)
502 int apex_node = -1;
503 if (mesh.is_pyramid(el))
504 for (int i = 0; i < ref_nodes.rows(); ++i)
505 if (std::abs(ref_nodes(i, 2) - 1.0) < 1e-12)
506 {
507 apex_node = i;
508 break;
509 }
510
511 for (int j = 0; j < lb.size(); ++j)
512 {
513 const int eid = lb.global_primitive_id(j);
514 const Eigen::VectorXi nodes = b.local_nodes_for_primitive(eid, mesh3d);
515 const int nfv = mesh3d.n_face_vertices(eid);
516 assert(nfv == 3 || nfv == 4);
517 assert(nodes.size() >= nfv);
518
519 // face corners come first in the local ordering (cyclic for quads);
520 // the same order as a next_around_face walk from the lf-scan index
521 // (tet/prism/pyramid/hex *_face_local_nodes all use it), so walking
522 // it again recovers each corner's global mesh vertex id
523 std::vector<Eigen::RowVector3d> c(nfv);
524 for (int k = 0; k < nfv; ++k)
525 c[k] = ref_nodes.row(nodes(k));
526
527 Navigation3D::Index nav;
528 for (int lf = 0; lf < mesh3d.n_cell_faces(el); ++lf)
529 {
530 nav = mesh3d.get_index_from_element(el, lf, 0);
531 if (nav.face == eid)
532 break;
533 }
534 assert(nav.face == eid);
535 std::vector<int> gv(nfv);
536 {
537 Navigation3D::Index cur = nav;
538 for (int k = 0; k < nfv; ++k)
539 {
540 gv[k] = cur.vertex;
541 cur = mesh3d.next_around_face(cur);
542 }
543 }
544
545 Eigen::MatrixXd pts;
546 std::vector<std::array<int, 2>> coords;
547 std::vector<std::tuple<int, int, int>> local_tris;
548 if (nfv == 3)
549 {
550 pts.resize((M + 1) * (M + 2) / 2, 3);
551 coords.resize(pts.rows());
552 std::vector<int> off(M + 2, 0);
553 for (int r = 0; r <= M; ++r)
554 off[r + 1] = off[r] + (M + 1 - r);
555 for (int r = 0; r <= M; ++r)
556 for (int i = 0; i <= M - r; ++i)
557 {
558 pts.row(off[r] + i) = c[0] + (double(i) / M) * (c[1] - c[0]) + (double(r) / M) * (c[2] - c[0]);
559 coords[off[r] + i] = {{i, r}};
560 }
561 for (int r = 0; r < M; ++r)
562 {
563 for (int i = 0; i < M - r; ++i)
564 {
565 local_tris.emplace_back(off[r] + i, off[r] + i + 1, off[r + 1] + i);
566 if (i + r < M - 1)
567 local_tris.emplace_back(off[r] + i + 1, off[r + 1] + i + 1, off[r + 1] + i);
568 }
569 }
570 }
571 else
572 {
573 pts.resize((M + 1) * (M + 1), 3);
574 coords.resize(pts.rows());
575 const auto gid = [M](const int i, const int r) { return r * (M + 1) + i; };
576 for (int r = 0; r <= M; ++r)
577 {
578 for (int i = 0; i <= M; ++i)
579 {
580 const double u = double(i) / M, v = double(r) / M;
581 pts.row(gid(i, r)) = (1 - u) * (1 - v) * c[0] + u * (1 - v) * c[1] + u * v * c[2] + (1 - u) * v * c[3];
582 coords[gid(i, r)] = {{i, r}};
583 }
584 }
585 for (int r = 0; r < M; ++r)
586 {
587 for (int i = 0; i < M; ++i)
588 {
589 local_tris.emplace_back(gid(i, r), gid(i + 1, r), gid(i, r + 1));
590 local_tris.emplace_back(gid(i + 1, r + 1), gid(i, r + 1), gid(i + 1, r));
591 }
592 }
593 }
594
595 std::vector<polyfem::assembler::AssemblyValues> vals;
596 b.evaluate_bases(pts, vals);
597
598 const auto edge_key = [M](const int va, const int vb, const int step) {
599 return va < vb ? std::array<int, 4>{{1, va, vb, step}}
600 : std::array<int, 4>{{1, vb, va, M - step}};
601 };
602
603 std::vector<int> ids(pts.rows());
604 for (int s = 0; s < pts.rows(); ++s)
605 {
606 const int i = coords[s][0], r = coords[s][1];
607 std::array<int, 4> key;
608 if (nfv == 3)
609 {
610 if (r == 0 && i == 0)
611 key = {{0, gv[0], 0, 0}};
612 else if (r == 0 && i == M)
613 key = {{0, gv[1], 0, 0}};
614 else if (r == M)
615 key = {{0, gv[2], 0, 0}};
616 else if (r == 0)
617 key = edge_key(gv[0], gv[1], i);
618 else if (i == 0)
619 key = edge_key(gv[0], gv[2], r);
620 else if (i + r == M)
621 key = edge_key(gv[1], gv[2], r);
622 else
623 key = {{2, eid, i, r}};
624 }
625 else
626 {
627 if (i == 0 && r == 0)
628 key = {{0, gv[0], 0, 0}};
629 else if (i == M && r == 0)
630 key = {{0, gv[1], 0, 0}};
631 else if (i == M && r == M)
632 key = {{0, gv[2], 0, 0}};
633 else if (i == 0 && r == M)
634 key = {{0, gv[3], 0, 0}};
635 else if (r == 0)
636 key = edge_key(gv[0], gv[1], i);
637 else if (i == M)
638 key = edge_key(gv[1], gv[2], r);
639 else if (r == M)
640 key = edge_key(gv[3], gv[2], i);
641 else if (i == 0)
642 key = edge_key(gv[0], gv[3], r);
643 else
644 key = {{2, eid, i, r}};
645 }
646
647 const auto it = vertex_id.find(key);
648 if (it != vertex_id.end())
649 {
650 ids[s] = it->second;
651 continue;
652 }
653
654 Eigen::Vector3d pos = Eigen::Vector3d::Zero();
655 std::map<int, double> weights;
656 if (apex_node >= 0 && std::abs(pts(s, 2) - 1.0) < 1e-12)
657 {
658 for (const auto &g : b.bases[apex_node].global())
659 {
660 pos += g.val * g.node.transpose();
661 weights[g.index] += g.val;
662 }
663 }
664 else
665 for (size_t i2 = 0; i2 < vals.size(); ++i2)
666 {
667 const double Ni = vals[i2].val(s);
668 if (std::abs(Ni) < 1e-12)
669 continue;
670 for (const auto &g : b.bases[i2].global())
671 {
672 pos += Ni * g.val * g.node.transpose();
673 weights[g.index] += Ni * g.val;
674 }
675 }
676 assert(pos.allFinite());
677
678 const int vid = int(vertices.size());
679 vertex_id[key] = vid;
680 vertices.push_back(pos);
681 for (const auto &kv : weights)
682 if (std::abs(kv.second) > 1e-10)
683 displacement_map_entries.emplace_back(vid, kv.first, kv.second);
684 ids[s] = vid;
685 }
686
687 for (const auto &t : local_tris)
688 proxy_tris.emplace_back(ids[std::get<0>(t)], ids[std::get<1>(t)], ids[std::get<2>(t)]);
689 }
690 }
691
692 node_positions.resize(vertices.size(), 3);
693 for (int i = 0; i < int(vertices.size()); ++i)
694 node_positions.row(i) = vertices[i];
695
696 boundary_triangles.resize(proxy_tris.size(), 3);
697 for (int i = 0; i < int(proxy_tris.size()); ++i)
698 boundary_triangles.row(i) << std::get<0>(proxy_tris[i]), std::get<2>(proxy_tris[i]), std::get<1>(proxy_tris[i]);
699
700 if (boundary_triangles.rows() > 0)
701 igl::edges(boundary_triangles, boundary_edges);
702
703 if (const char *dump = getenv("POLYFEM_DUMP_COLLISION_PROXY"))
704 igl::write_triangle_mesh(dump, node_positions, boundary_triangles);
705 }
706
707 void OutGeometryData::extract_boundary_mesh(
708 const mesh::Mesh &mesh,
709 const int n_bases,
710 const std::vector<basis::ElementBases> &bases,
711 const std::vector<mesh::LocalBoundary> &total_local_boundary,
712 Eigen::MatrixXd &node_positions,
713 Eigen::MatrixXi &boundary_edges,
714 Eigen::MatrixXi &boundary_triangles,
715 std::vector<Eigen::Triplet<double>> &displacement_map_entries)
716 {
717 using namespace polyfem::mesh;
718
719 displacement_map_entries.clear();
720
721 if (mesh.is_volume())
722 {
723 if (mesh.has_poly())
724 {
725 logger().warn("Skipping as the mesh has polygons");
726 return;
727 }
728
729 const bool is_simplicial = mesh.is_simplicial();
730
731 std::vector<Eigen::Vector3d> node_positions_vec;
732 node_positions_vec.reserve(n_bases + (is_simplicial ? 0 : mesh.n_faces()));
733
734 // node_positions.resize(n_bases + (is_simplicial ? 0 : mesh.n_faces()), 3);
735 // node_positions.setZero();
736 const Mesh3D &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
737
738 std::vector<std::tuple<int, int, int>> tris;
739
740 std::vector<bool> visited_node(n_bases, false);
741
742 std::stringstream print_warning;
743
744 // Hybrid (prism/pyramid) meshes: element orders can differ across an
745 // interface (anisotropic prisms promote their tet/pyramid neighbors),
746 // so tessellating each boundary face at its own order produces
747 // T-junctions along shared edges, and stitched interface nodes have
748 // no DOF of their own. Build the proxy from the DOFs on the
749 // boundary instead: every proxy vertex IS a global DOF (weight-1
750 // displacement map), each face edge is subdivided by the DOFs
751 // located on it -- a property of the edge, not of the element
752 // looking at it, so neighboring faces conform by construction --
753 // and face interiors by the element's own owned face nodes.
754 // Stitched nodes are skipped; the DOFs they depend on subdivide the
755 // edge instead.
756 bool has_prism_or_pyramid = false;
757 for (const LocalBoundary &lb : total_local_boundary)
758 {
759 if (mesh.is_prism(lb.element_id()) || mesh.is_pyramid(lb.element_id()))
760 {
761 has_prism_or_pyramid = true;
762 break;
763 }
764 }
765
766 if (has_prism_or_pyramid && mesh3d.is_conforming())
767 {
768 const auto edge_dofs = build_edge_dofs(mesh, bases);
769 constexpr long S = PROXY_PARAM_SCALE;
770
771 const auto emit_dof = [&](const int gindex, const Eigen::Vector3d &pos) {
772 if (gindex >= int(node_positions_vec.size()))
773 node_positions_vec.resize(gindex + 1, Eigen::Vector3d::Zero());
774 node_positions_vec[gindex] = pos;
775 if (!visited_node[gindex])
776 displacement_map_entries.emplace_back(gindex, gindex, 1);
777 visited_node[gindex] = true;
778 };
779
780 Eigen::MatrixXd ref_nodes;
781 for (const LocalBoundary &lb : total_local_boundary)
782 {
783 const int el = lb.element_id();
784 const basis::ElementBases &b = bases[el];
785 if (b.bases.empty() || !element_ref_nodes(mesh, el, b, ref_nodes))
786 continue;
787
788 for (int j = 0; j < lb.size(); ++j)
789 {
790 const int eid = lb.global_primitive_id(j);
791 const Eigen::VectorXi nodes = b.local_nodes_for_primitive(eid, mesh3d);
792 const int nfv = mesh3d.n_face_vertices(eid);
793 assert(nfv == 3 || nfv == 4);
794 assert(nodes.size() >= nfv);
795
796 // exact face parameters (units of 1/S) and the DOF of each point
797 std::vector<std::array<long, 2>> pts;
798 std::vector<int> dof;
799
800 // face corners come first in the local ordering (cyclic for quads)
801 std::array<std::array<long, 2>, 4> cp;
802 if (nfv == 3)
803 cp = {{{{0, 0}}, {{S, 0}}, {{0, S}}, {{0, 0}}}};
804 else
805 cp = {{{{0, 0}}, {{S, 0}}, {{S, S}}, {{0, S}}}};
806 std::array<int, 4> cd{{-1, -1, -1, -1}};
807 bool ok = true;
808 for (int k = 0; k < nfv; ++k)
809 {
810 const auto &glob = b.bases[nodes(k)].global();
811 assert(glob.size() == 1); // face corners always own their DOF
812 if (glob.size() != 1)
813 {
814 ok = false;
815 break;
816 }
817 cd[k] = glob.front().index;
818 pts.push_back(cp[k]);
819 dof.push_back(cd[k]);
820 emit_dof(cd[k], glob.front().node.transpose());
821 }
822 if (!ok)
823 continue;
824
825 // edge subdivision: the DOFs located on each face edge
826 for (int k = 0; k < nfv; ++k)
827 {
828 const int va = cd[k], vb = cd[(k + 1) % nfv];
829 const auto it = edge_dofs.find({std::min(va, vb), std::max(va, vb)});
830 if (it == edge_dofs.end())
831 continue;
832 for (const auto &ed : it->second)
833 {
834 const long s = va < vb ? std::get<0>(ed) : S - std::get<0>(ed);
835 const auto &A = cp[k];
836 const auto &B = cp[(k + 1) % nfv];
837 pts.push_back({{(A[0] * (S - s) + B[0] * s) / S, (A[1] * (S - s) + B[1] * s) / S}});
838 dof.push_back(std::get<1>(ed));
839 emit_dof(std::get<1>(ed), std::get<2>(ed));
840 }
841 }
842
843 // owned face-interior nodes at the element's own resolution
844 // (reference faces are planar parallelograms -> affine map)
845 const Eigen::RowVector3d c0 = ref_nodes.row(nodes(0));
846 const Eigen::RowVector3d A3 = ref_nodes.row(nodes(1)) - c0;
847 const Eigen::RowVector3d B3 = ref_nodes.row(nodes(nfv - 1)) - c0;
848 const double aa = A3.squaredNorm(), bb = B3.squaredNorm(), ab = A3.dot(B3);
849 const double det = aa * bb - ab * ab;
850 for (long n = nfv; n < nodes.size(); ++n)
851 {
852 const auto &glob = b.bases[nodes(n)].global();
853 if (glob.size() != 1)
854 continue; // stitched interface node
855 const Eigen::RowVector3d d3 = ref_nodes.row(nodes(n)) - c0;
856 const double du = d3.dot(A3), dv = d3.dot(B3);
857 const double u = (du * bb - dv * ab) / det;
858 const double v = (dv * aa - du * ab) / det;
859 const long lu = std::lround(u * S), lv = std::lround(v * S);
860 assert(std::abs(u * S - lu) < 1e-6 && std::abs(v * S - lv) < 1e-6);
861 // nodes on a face edge already subdivide the edge above
862 const bool on_edge = nfv == 3
863 ? (lu == 0 || lv == 0 || lu + lv == S)
864 : (lu == 0 || lu == S || lv == 0 || lv == S);
865 if (on_edge)
866 continue;
867 pts.push_back({{lu, lv}});
868 dof.push_back(glob.front().index);
869 emit_dof(glob.front().index, glob.front().node.transpose());
870 }
871
872 std::vector<std::array<int, 3>> local_tris;
873 if (!triangulate_lattice(pts, nfv, local_tris))
874 triangulate_convex_pointset(pts, local_tris);
875 for (const auto &t : local_tris)
876 tris.emplace_back(dof[t[0]], dof[t[1]], dof[t[2]]);
877 }
878 }
879
880 // downstream consumers (e.g. the shape-derivative code) expect every
881 // FE node to have a row, as the pre-split extraction guaranteed
882 node_positions_vec.resize(
883 std::max(node_positions_vec.size(), size_t(n_bases)), Eigen::Vector3d::Zero());
884
885 node_positions.resize(node_positions_vec.size(), 3);
886 for (int i = 0; i < int(node_positions_vec.size()); ++i)
887 node_positions.row(i) = node_positions_vec[i];
888
889 boundary_triangles.resize(tris.size(), 3);
890 for (int i = 0; i < int(tris.size()); ++i)
891 boundary_triangles.row(i) << std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]);
892
893 if (boundary_triangles.rows() > 0)
894 igl::edges(boundary_triangles, boundary_edges);
895
896 if (const char *dump = getenv("POLYFEM_DUMP_COLLISION_PROXY"))
897 igl::write_triangle_mesh(dump, node_positions, boundary_triangles);
898
899 return;
900 }
901
902 for (const LocalBoundary &lb : total_local_boundary)
903 {
904 const basis::ElementBases &b = bases[lb.element_id()];
905
906 for (int j = 0; j < lb.size(); ++j)
907 {
908 const int eid = lb.global_primitive_id(j);
909 const int lid = lb[j];
910 const Eigen::VectorXi nodes = b.local_nodes_for_primitive(eid, mesh3d);
911
912 if (mesh.is_cube(lb.element_id()))
913 {
914 assert(!is_simplicial);
915 assert(!mesh.has_poly());
916 std::vector<int> loc_nodes;
917 RowVectorNd bary = RowVectorNd::Zero(3);
918
919 for (long n = 0; n < nodes.size(); ++n)
920 {
921 auto &bs = b.bases[nodes(n)];
922 const auto &glob = bs.global();
923 if (glob.size() != 1)
924 continue;
925
926 int gindex = glob.front().index;
927 node_positions_vec.resize(std::max(int(node_positions_vec.size()), gindex + 1));
928 node_positions_vec[gindex] = glob.front().node;
929 bary += glob.front().node;
930 loc_nodes.push_back(gindex);
931 }
932
933 if (loc_nodes.size() != 4)
934 {
935 logger().trace("skipping element {} since it is not Q1", eid);
936 continue;
937 }
938
939 bary /= 4;
940
941 const int new_node = n_bases + eid;
942 node_positions_vec.resize(std::max(int(node_positions_vec.size()), new_node + 1));
943 node_positions_vec[new_node] = bary;
944 tris.emplace_back(loc_nodes[1], loc_nodes[0], new_node);
945 tris.emplace_back(loc_nodes[2], loc_nodes[1], new_node);
946 tris.emplace_back(loc_nodes[3], loc_nodes[2], new_node);
947 tris.emplace_back(loc_nodes[0], loc_nodes[3], new_node);
948
949 for (int q = 0; q < 4; ++q)
950 {
951 if (!visited_node[loc_nodes[q]])
952 displacement_map_entries.emplace_back(loc_nodes[q], loc_nodes[q], 1);
953
954 visited_node[loc_nodes[q]] = true;
955 displacement_map_entries.emplace_back(new_node, loc_nodes[q], 0.25);
956 }
957
958 continue;
959 }
960 else if (mesh.is_prism(lb.element_id()))
961 {
962 assert(!is_simplicial);
963 assert(!mesh.has_poly());
964 std::vector<int> loc_nodes;
965 std::vector<int> loc_local_nodes;
966
967 for (long n = 0; n < nodes.size(); ++n)
968 {
969 auto &bs = b.bases[nodes(n)];
970 const auto &glob = bs.global();
971 if (glob.size() != 1)
972 continue;
973
974 int gindex = glob.front().index;
975 node_positions_vec.resize(std::max(int(node_positions_vec.size()), gindex + 1));
976 node_positions_vec[gindex] = glob.front().node;
977 loc_nodes.push_back(gindex);
978 loc_local_nodes.push_back(nodes(n));
979 }
980
981 auto update_mapping = [&displacement_map_entries, &visited_node](const std::vector<int> &loc_nodes) {
982 for (int k = 0; k < loc_nodes.size(); ++k)
983 {
984 if (!visited_node[loc_nodes[k]])
985 displacement_map_entries.emplace_back(loc_nodes[k], loc_nodes[k], 1);
986
987 visited_node[loc_nodes[k]] = true;
988 }
989 };
990
991 // tri face
992 if (lid < 2)
993 {
994 if (loc_nodes.size() == 3)
995 {
996 tris.emplace_back(loc_nodes[0], loc_nodes[1], loc_nodes[2]);
997
998 update_mapping(loc_nodes);
999 }
1000 else if (loc_nodes.size() == 6)
1001 {
1002 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[5]);
1003 tris.emplace_back(loc_nodes[3], loc_nodes[1], loc_nodes[4]);
1004 tris.emplace_back(loc_nodes[4], loc_nodes[2], loc_nodes[5]);
1005 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[5]);
1006
1007 update_mapping(loc_nodes);
1008 }
1009 else if (loc_nodes.size() == 10)
1010 {
1011 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[8]);
1012 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[9]);
1013 tris.emplace_back(loc_nodes[4], loc_nodes[1], loc_nodes[5]);
1014 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[9]);
1015 tris.emplace_back(loc_nodes[6], loc_nodes[2], loc_nodes[7]);
1016 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[9]);
1017 tris.emplace_back(loc_nodes[8], loc_nodes[3], loc_nodes[9]);
1018 tris.emplace_back(loc_nodes[9], loc_nodes[4], loc_nodes[5]);
1019 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[9]);
1020 update_mapping(loc_nodes);
1021 }
1022 else
1023 {
1024 logger().trace("skipping element {} since it is not linear, it has {} nodes", eid, loc_nodes.size());
1025 }
1026 }
1027 else
1028 {
1029 if (loc_nodes.size() < 4 || loc_local_nodes.size() < 4)
1030 {
1031 logger().trace("skipping prism quad face {} since it has only {} complete nodes", eid, loc_nodes.size());
1032 continue;
1033 }
1034
1035 const int p = b.bases.empty() ? -1 : b.bases.front().order();
1036 const int n_tri_nodes = (p + 1) * (p + 2) / 2;
1037 const int q = n_tri_nodes > 0 && b.bases.size() % n_tri_nodes == 0 ? int(b.bases.size()) / n_tri_nodes - 1 : -1;
1038
1039 if (p < 1 || p > 3 || q < 1 || q > 3 || (p == 3 && q == 3))
1040 {
1041 logger().trace("skipping prism quad face {} with unsupported p={}, q={}", eid, p, q);
1042 continue;
1043 }
1044
1045 auto is_vertical_prism_edge = [](const int a, const int b) {
1046 return (a >= 0 && a < 3 && b == a + 3) || (b >= 0 && b < 3 && a == b + 3);
1047 };
1048
1049 std::vector<int> edge_orders(4);
1050 for (int k = 0; k < 4; ++k)
1051 edge_orders[k] = is_vertical_prism_edge(loc_local_nodes[k], loc_local_nodes[(k + 1) % 4]) ? q : p;
1052
1053 const int u_order = edge_orders[0];
1054 const int v_order = edge_orders[1];
1055 const int expected_nodes = (u_order + 1) * (v_order + 1);
1056 if (loc_nodes.size() != expected_nodes || edge_orders[0] != edge_orders[2] || edge_orders[1] != edge_orders[3])
1057 {
1058 logger().trace("skipping prism quad face {} with p={}, q={} and {} nodes", eid, p, q, loc_nodes.size());
1059 continue;
1060 }
1061
1062 std::vector<int> grid(expected_nodes, -1);
1063 auto grid_index = [u_order](const int i, const int j) {
1064 return j * (u_order + 1) + i;
1065 };
1066
1067 grid[grid_index(0, 0)] = loc_nodes[0];
1068 grid[grid_index(u_order, 0)] = loc_nodes[1];
1069 grid[grid_index(u_order, v_order)] = loc_nodes[2];
1070 grid[grid_index(0, v_order)] = loc_nodes[3];
1071
1072 int node_index = 4;
1073 for (int i = 1; i < u_order; ++i)
1074 grid[grid_index(i, 0)] = loc_nodes[node_index++];
1075 for (int j = 1; j < v_order; ++j)
1076 grid[grid_index(u_order, j)] = loc_nodes[node_index++];
1077 for (int i = u_order - 1; i > 0; --i)
1078 grid[grid_index(i, v_order)] = loc_nodes[node_index++];
1079 for (int j = v_order - 1; j > 0; --j)
1080 grid[grid_index(0, j)] = loc_nodes[node_index++];
1081
1082 for (int j = 1; j < v_order; ++j)
1083 for (int i = 1; i < u_order; ++i)
1084 grid[grid_index(i, j)] = loc_nodes[node_index++];
1085
1086 assert(node_index == loc_nodes.size());
1087 assert(std::all_of(grid.begin(), grid.end(), [](const int n) { return n >= 0; }));
1088
1089 for (int j = 0; j < v_order; ++j)
1090 {
1091 for (int i = 0; i < u_order; ++i)
1092 {
1093 tris.emplace_back(grid[grid_index(i, j)], grid[grid_index(i + 1, j)], grid[grid_index(i, j + 1)]);
1094 tris.emplace_back(grid[grid_index(i + 1, j + 1)], grid[grid_index(i, j + 1)], grid[grid_index(i + 1, j)]);
1095 }
1096 }
1097
1098 update_mapping(loc_nodes);
1099 }
1100
1101 continue;
1102 }
1103 else if (mesh.is_pyramid(lb.element_id()))
1104 {
1105 assert(!is_simplicial);
1106 assert(!mesh.has_poly());
1107 std::vector<int> loc_nodes;
1108 std::vector<int> loc_local_nodes;
1109
1110 for (long n = 0; n < nodes.size(); ++n)
1111 {
1112 auto &bs = b.bases[nodes(n)];
1113 const auto &glob = bs.global();
1114 if (glob.size() != 1)
1115 continue;
1116
1117 int gindex = glob.front().index;
1118 node_positions_vec.resize(std::max(int(node_positions_vec.size()), gindex + 1));
1119 node_positions_vec[gindex] = glob.front().node;
1120 loc_nodes.push_back(gindex);
1121 loc_local_nodes.push_back(nodes(n));
1122 }
1123
1124 auto update_mapping = [&displacement_map_entries, &visited_node](const std::vector<int> &loc_nodes) {
1125 for (int k = 0; k < loc_nodes.size(); ++k)
1126 {
1127 if (!visited_node[loc_nodes[k]])
1128 displacement_map_entries.emplace_back(loc_nodes[k], loc_nodes[k], 1);
1129
1130 visited_node[loc_nodes[k]] = true;
1131 }
1132 };
1133
1134 const int p = b.bases.empty() ? -1 : b.bases.front().order();
1135 if (p < 1 || p > 3)
1136 {
1137 logger().trace("skipping pyramid face {} with unsupported p={}", eid, p);
1138 continue;
1139 }
1140
1141 if (lid == 0)
1142 {
1143 const int expected_nodes = (p + 1) * (p + 1);
1144 if (loc_nodes.size() != expected_nodes || loc_local_nodes.size() != expected_nodes)
1145 {
1146 logger().trace("skipping pyramid quad face {} with p={} and {} nodes", eid, p, loc_nodes.size());
1147 continue;
1148 }
1149
1150 Eigen::MatrixXd pyramid_nodes;
1151 autogen::pyramid_nodes_3d(p, pyramid_nodes);
1152
1153 const Eigen::RowVector3d origin = pyramid_nodes.row(loc_local_nodes[0]);
1154 const Eigen::RowVector3d u_axis = pyramid_nodes.row(loc_local_nodes[1]) - origin;
1155 const Eigen::RowVector3d v_axis = pyramid_nodes.row(loc_local_nodes[3]) - origin;
1156
1157 std::vector<int> grid(expected_nodes, -1);
1158 auto grid_index = [p](const int i, const int j) {
1159 return j * (p + 1) + i;
1160 };
1161
1162 bool valid_grid = true;
1163 for (int n = 0; n < loc_nodes.size(); ++n)
1164 {
1165 const Eigen::RowVector3d rel = pyramid_nodes.row(loc_local_nodes[n]) - origin;
1166 const int i = int(std::lround(p * rel.dot(u_axis) / u_axis.squaredNorm()));
1167 const int j = int(std::lround(p * rel.dot(v_axis) / v_axis.squaredNorm()));
1168 if (i < 0 || i > p || j < 0 || j > p)
1169 {
1170 logger().trace("skipping pyramid quad face {} with invalid local grid coordinate ({}, {})", eid, i, j);
1171 valid_grid = false;
1172 break;
1173 }
1174 if (grid[grid_index(i, j)] >= 0)
1175 {
1176 logger().trace("skipping pyramid quad face {} with duplicate local grid coordinate ({}, {})", eid, i, j);
1177 valid_grid = false;
1178 break;
1179 }
1180 grid[grid_index(i, j)] = loc_nodes[n];
1181 }
1182
1183 if (!valid_grid || !std::all_of(grid.begin(), grid.end(), [](const int n) { return n >= 0; }))
1184 continue;
1185
1186 for (int j = 0; j < p; ++j)
1187 {
1188 for (int i = 0; i < p; ++i)
1189 {
1190 tris.emplace_back(grid[grid_index(i, j)], grid[grid_index(i + 1, j)], grid[grid_index(i, j + 1)]);
1191 tris.emplace_back(grid[grid_index(i + 1, j + 1)], grid[grid_index(i, j + 1)], grid[grid_index(i + 1, j)]);
1192 }
1193 }
1194
1195 update_mapping(loc_nodes);
1196 }
1197 else if (loc_nodes.size() == 3)
1198 {
1199 tris.emplace_back(loc_nodes[0], loc_nodes[1], loc_nodes[2]);
1200 update_mapping(loc_nodes);
1201 }
1202 else if (loc_nodes.size() == 6)
1203 {
1204 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[5]);
1205 tris.emplace_back(loc_nodes[3], loc_nodes[1], loc_nodes[4]);
1206 tris.emplace_back(loc_nodes[4], loc_nodes[2], loc_nodes[5]);
1207 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[5]);
1208 update_mapping(loc_nodes);
1209 }
1210 else if (loc_nodes.size() == 10)
1211 {
1212 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[8]);
1213 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[9]);
1214 tris.emplace_back(loc_nodes[4], loc_nodes[1], loc_nodes[5]);
1215 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[9]);
1216 tris.emplace_back(loc_nodes[6], loc_nodes[2], loc_nodes[7]);
1217 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[9]);
1218 tris.emplace_back(loc_nodes[8], loc_nodes[3], loc_nodes[9]);
1219 tris.emplace_back(loc_nodes[9], loc_nodes[4], loc_nodes[5]);
1220 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[9]);
1221 update_mapping(loc_nodes);
1222 }
1223 else
1224 {
1225 logger().trace("skipping pyramid tri face {} with p={} and {} nodes", eid, p, loc_nodes.size());
1226 continue;
1227 }
1228
1229 continue;
1230 }
1231
1232 if (!mesh.is_simplex(lb.element_id()))
1233 {
1234 logger().trace("skipping element {} since it is not a simplex or hex", eid);
1235 continue;
1236 }
1237
1238 assert(mesh.is_simplex(lb.element_id()));
1239
1240 std::vector<int> loc_nodes;
1241
1242 bool is_follower = false;
1243 if (!mesh3d.is_conforming())
1244 {
1245 for (long n = 0; n < nodes.size(); ++n)
1246 {
1247 auto &bs = b.bases[nodes(n)];
1248 const auto &glob = bs.global();
1249 if (glob.size() != 1)
1250 {
1251 is_follower = true;
1252 break;
1253 }
1254 }
1255 }
1256
1257 if (is_follower)
1258 continue;
1259
1260 for (long n = 0; n < nodes.size(); ++n)
1261 {
1262 const basis::Basis &bs = b.bases[nodes(n)];
1263 const std::vector<basis::Local2Global> &glob = bs.global();
1264 if (glob.size() != 1)
1265 continue;
1266
1267 int gindex = glob.front().index;
1268 node_positions_vec.resize(std::max(int(node_positions_vec.size()), gindex + 1));
1269 node_positions_vec[gindex] = glob.front().node;
1270 loc_nodes.push_back(gindex);
1271 }
1272
1273 if (loc_nodes.size() == 3)
1274 {
1275 tris.emplace_back(loc_nodes[0], loc_nodes[1], loc_nodes[2]);
1276 }
1277 else if (loc_nodes.size() == 6)
1278 {
1279 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[5]);
1280 tris.emplace_back(loc_nodes[3], loc_nodes[1], loc_nodes[4]);
1281 tris.emplace_back(loc_nodes[4], loc_nodes[2], loc_nodes[5]);
1282 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[5]);
1283 }
1284 else if (loc_nodes.size() == 10)
1285 {
1286 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[8]);
1287 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[9]);
1288 tris.emplace_back(loc_nodes[4], loc_nodes[1], loc_nodes[5]);
1289 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[9]);
1290 tris.emplace_back(loc_nodes[6], loc_nodes[2], loc_nodes[7]);
1291 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[9]);
1292 tris.emplace_back(loc_nodes[8], loc_nodes[3], loc_nodes[9]);
1293 tris.emplace_back(loc_nodes[9], loc_nodes[4], loc_nodes[5]);
1294 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[9]);
1295 }
1296 else if (loc_nodes.size() == 15)
1297 {
1298 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[11]);
1299 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[12]);
1300 tris.emplace_back(loc_nodes[3], loc_nodes[12], loc_nodes[11]);
1301 tris.emplace_back(loc_nodes[12], loc_nodes[10], loc_nodes[11]);
1302 tris.emplace_back(loc_nodes[4], loc_nodes[5], loc_nodes[13]);
1303 tris.emplace_back(loc_nodes[4], loc_nodes[13], loc_nodes[12]);
1304 tris.emplace_back(loc_nodes[12], loc_nodes[13], loc_nodes[14]);
1305 tris.emplace_back(loc_nodes[12], loc_nodes[14], loc_nodes[10]);
1306 tris.emplace_back(loc_nodes[14], loc_nodes[9], loc_nodes[10]);
1307 tris.emplace_back(loc_nodes[5], loc_nodes[1], loc_nodes[6]);
1308 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[13]);
1309 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[13]);
1310 tris.emplace_back(loc_nodes[13], loc_nodes[7], loc_nodes[14]);
1311 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[14]);
1312 tris.emplace_back(loc_nodes[14], loc_nodes[8], loc_nodes[9]);
1313 tris.emplace_back(loc_nodes[8], loc_nodes[2], loc_nodes[9]);
1314 }
1315 else
1316 {
1317 print_warning << loc_nodes.size() << " ";
1318 // assert(false);
1319 }
1320
1321 if (!is_simplicial)
1322 {
1323 for (int k = 0; k < loc_nodes.size(); ++k)
1324 {
1325 if (!visited_node[loc_nodes[k]])
1326 displacement_map_entries.emplace_back(loc_nodes[k], loc_nodes[k], 1);
1327
1328 visited_node[loc_nodes[k]] = true;
1329 }
1330 }
1331 }
1332 }
1333
1334 if (print_warning.str().size() > 0)
1335 logger().warn("Skipping faces as theys have {} nodes, boundary export supported up to p4", print_warning.str());
1336
1337 // downstream consumers (e.g. the shape-derivative code) expect every
1338 // FE node to have a row, as the pre-split extraction guaranteed
1339 node_positions_vec.resize(
1340 std::max(node_positions_vec.size(), size_t(n_bases + (is_simplicial ? 0 : mesh.n_faces()))),
1341 Eigen::Vector3d::Zero());
1342
1343 node_positions.resize(node_positions_vec.size(), 3);
1344 for (int i = 0; i < node_positions_vec.size(); ++i)
1345 node_positions.row(i) = node_positions_vec[i];
1346
1347 boundary_triangles.resize(tris.size(), 3);
1348 for (int i = 0; i < tris.size(); ++i)
1349 {
1350 boundary_triangles.row(i) << std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]);
1351 }
1352
1353 if (boundary_triangles.rows() > 0)
1354 {
1355 igl::edges(boundary_triangles, boundary_edges);
1356 }
1357
1358 if (const char *dump = getenv("POLYFEM_DUMP_COLLISION_PROXY"))
1359 igl::write_triangle_mesh(dump, node_positions, boundary_triangles);
1360 }
1361 else
1362 {
1363 node_positions.resize(n_bases, 2);
1364 node_positions.setZero();
1365 const Mesh2D &mesh2d = dynamic_cast<const Mesh2D &>(mesh);
1366
1367 std::vector<std::pair<int, int>> edges;
1368
1369 for (const LocalBoundary &lb : total_local_boundary)
1370 {
1371 const basis::ElementBases &b = bases[lb.element_id()];
1372
1373 for (int j = 0; j < lb.size(); ++j)
1374 {
1375 const int eid = lb.global_primitive_id(j);
1376 const int lid = lb[j];
1377 const Eigen::VectorXi nodes = b.local_nodes_for_primitive(eid, mesh2d);
1378
1379 int prev_node = -1;
1380
1381 for (long n = 0; n < nodes.size(); ++n)
1382 {
1383 const basis::Basis &bs = b.bases[nodes(n)];
1384 const std::vector<basis::Local2Global> &glob = bs.global();
1385 if (glob.size() != 1)
1386 continue;
1387
1388 int gindex = glob.front().index;
1389 node_positions.row(gindex) = glob.front().node.head<2>();
1390
1391 if (prev_node >= 0)
1392 edges.emplace_back(prev_node, gindex);
1393
1394 prev_node = gindex;
1395 }
1396 }
1397 }
1398
1399 boundary_triangles.resize(0, 0);
1400 boundary_edges.resize(edges.size(), 2);
1401 for (int i = 0; i < edges.size(); ++i)
1402 {
1403 boundary_edges.row(i) << edges[i].first, edges[i].second;
1404 }
1405 }
1406 }
1407
1408 void OutGeometryData::build_vis_boundary_mesh(
1409 const mesh::Mesh &mesh,
1410 const std::vector<basis::ElementBases> &gbases,
1411 const std::vector<mesh::LocalBoundary> &total_local_boundary,
1412 Eigen::MatrixXd &boundary_vis_vertices,
1413 Eigen::MatrixXd &boundary_vis_local_vertices,
1414 Eigen::MatrixXi &boundary_vis_elements,
1415 Eigen::MatrixXi &boundary_vis_elements_ids,
1416 Eigen::MatrixXi &boundary_vis_primitive_ids,
1417 Eigen::MatrixXd &boundary_vis_normals) const
1418 {
1419 using namespace polyfem::mesh;
1420
1421 std::vector<Eigen::MatrixXd> lv, vertices, allnormals;
1422 std::vector<int> el_ids, global_primitive_ids;
1423 Eigen::MatrixXd uv, local_pts, tmp_n, normals;
1424 assembler::ElementAssemblyValues vals;
1425 const auto &sampler = ref_element_sampler;
1426 const int n_samples = sampler.num_samples();
1427 int size = 0;
1428
1429 std::vector<std::pair<int, int>> edges;
1430 std::vector<std::tuple<int, int, int>> tris;
1431
1432 for (auto it = total_local_boundary.begin(); it != total_local_boundary.end(); ++it)
1433 {
1434 const auto &lb = *it;
1435 const auto &gbs = gbases[lb.element_id()];
1436
1437 for (int k = 0; k < lb.size(); ++k)
1438 {
1439 switch (lb.type())
1440 {
1441 case BoundaryType::TRI_LINE:
1442 utils::BoundarySampler::normal_for_tri_edge(lb[k], tmp_n);
1443 utils::BoundarySampler::sample_parametric_tri_edge(lb[k], n_samples, uv, local_pts);
1444 break;
1445 case BoundaryType::QUAD_LINE:
1446 utils::BoundarySampler::normal_for_quad_edge(lb[k], tmp_n);
1447 utils::BoundarySampler::sample_parametric_quad_edge(lb[k], n_samples, uv, local_pts);
1448 break;
1449 case BoundaryType::QUAD:
1450 utils::BoundarySampler::normal_for_quad_face(lb[k], tmp_n);
1451 utils::BoundarySampler::sample_parametric_quad_face(lb[k], n_samples, uv, local_pts);
1452 break;
1453 case BoundaryType::TRI:
1454 utils::BoundarySampler::normal_for_tri_face(lb[k], tmp_n);
1455 utils::BoundarySampler::sample_parametric_tri_face(lb[k], n_samples, uv, local_pts);
1456 break;
1457 case BoundaryType::PRISM:
1458 utils::BoundarySampler::normal_for_prism_face(lb[k], tmp_n);
1459 utils::BoundarySampler::sample_parametric_prism_face(lb[k], n_samples, uv, local_pts);
1460 break;
1461 case BoundaryType::PYRAMID:
1462 utils::BoundarySampler::normal_for_pyramid_face(lb[k], tmp_n);
1463 utils::BoundarySampler::sample_parametric_pyramid_face(lb[k], n_samples, uv, local_pts);
1464 break;
1465 case BoundaryType::POLYGON:
1466 utils::BoundarySampler::normal_for_polygon_edge(lb.element_id(), lb.global_primitive_id(k), mesh, tmp_n);
1467 utils::BoundarySampler::sample_polygon_edge(lb.element_id(), lb.global_primitive_id(k), n_samples, mesh, uv, local_pts);
1468 break;
1469 case BoundaryType::POLYHEDRON:
1470 assert(false);
1471 break;
1472 case BoundaryType::INVALID:
1473 assert(false);
1474 break;
1475 default:
1476 assert(false);
1477 }
1478
1479 vertices.emplace_back();
1480 lv.emplace_back(local_pts);
1481 el_ids.push_back(lb.element_id());
1482 global_primitive_ids.push_back(lb.global_primitive_id(k));
1483 gbs.eval_geom_mapping(local_pts, vertices.back());
1484 vals.compute(lb.element_id(), mesh.is_volume(), local_pts, gbs, gbs);
1485 const int tris_start = tris.size();
1486
1487 if (mesh.is_volume())
1488 {
1489 const bool prism_quad = lb.type() == BoundaryType::PRISM && lb[k] >= 2;
1490 const bool prism_tri = lb.type() == BoundaryType::PRISM && lb[k] < 2;
1491
1492 const bool pyramid_quad = lb.type() == BoundaryType::PYRAMID && lb[k] == 0;
1493 const bool pyramid_tri = lb.type() == BoundaryType::PYRAMID && lb[k] > 0;
1494
1495 if (lb.type() == BoundaryType::QUAD || prism_quad || pyramid_quad)
1496 {
1497 const auto map = [n_samples, size](int i, int j) { return j * n_samples + i + size; };
1498
1499 for (int j = 0; j < n_samples - 1; ++j)
1500 {
1501 for (int i = 0; i < n_samples - 1; ++i)
1502 {
1503 tris.emplace_back(map(i, j), map(i + 1, j), map(i, j + 1));
1504 tris.emplace_back(map(i + 1, j + 1), map(i, j + 1), map(i + 1, j));
1505 }
1506 }
1507 }
1508 else if (lb.type() == BoundaryType::TRI || prism_tri || pyramid_tri)
1509 {
1510 int index = 0;
1511 std::vector<int> mapp(n_samples * n_samples, -1);
1512 for (int j = 0; j < n_samples; ++j)
1513 {
1514 for (int i = 0; i < n_samples - j; ++i)
1515 {
1516 mapp[j * n_samples + i] = index;
1517 ++index;
1518 }
1519 }
1520 const auto map = [mapp, n_samples](int i, int j) {
1521 if (j * n_samples + i >= mapp.size())
1522 return -1;
1523 return mapp[j * n_samples + i];
1524 };
1525
1526 for (int j = 0; j < n_samples - 1; ++j)
1527 {
1528 for (int i = 0; i < n_samples - j; ++i)
1529 {
1530 if (map(i, j) >= 0 && map(i + 1, j) >= 0 && map(i, j + 1) >= 0)
1531 tris.emplace_back(map(i, j) + size, map(i + 1, j) + size, map(i, j + 1) + size);
1532
1533 if (map(i + 1, j + 1) >= 0 && map(i, j + 1) >= 0 && map(i + 1, j) >= 0)
1534 tris.emplace_back(map(i + 1, j + 1) + size, map(i, j + 1) + size, map(i + 1, j) + size);
1535 }
1536 }
1537 }
1538 else
1539 {
1540 assert(false);
1541 }
1542 }
1543 else
1544 {
1545 for (int i = 0; i < vertices.back().rows() - 1; ++i)
1546 edges.emplace_back(i + size, i + size + 1);
1547 }
1548
1549 normals.resize(vals.jac_it.size(), tmp_n.cols());
1550
1551 for (int n = 0; n < vals.jac_it.size(); ++n)
1552 {
1553 normals.row(n) = tmp_n * vals.jac_it[n];
1554 normals.row(n).normalize();
1555 }
1556
1557 allnormals.push_back(normals);
1558
1559 tmp_n.setZero();
1560 for (int n = 0; n < vals.jac_it.size(); ++n)
1561 {
1562 tmp_n += normals.row(n);
1563 }
1564
1565 if (mesh.is_volume())
1566 {
1567 Eigen::Vector3d e1 = vertices.back().row(std::get<1>(tris.back()) - size) - vertices.back().row(std::get<0>(tris.back()) - size);
1568 Eigen::Vector3d e2 = vertices.back().row(std::get<2>(tris.back()) - size) - vertices.back().row(std::get<0>(tris.back()) - size);
1569
1570 Eigen::Vector3d n = e1.cross(e2);
1571 Eigen::Vector3d nn = tmp_n.transpose();
1572
1573 if (n.dot(nn) < 0)
1574 {
1575 for (int i = tris_start; i < tris.size(); ++i)
1576 {
1577 tris[i] = std::tuple<int, int, int>(std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]));
1578 }
1579 }
1580 }
1581
1582 size += vertices.back().rows();
1583 }
1584 }
1585
1586 boundary_vis_vertices.resize(size, vertices.front().cols());
1587 boundary_vis_local_vertices.resize(size, vertices.front().cols());
1588 boundary_vis_elements_ids.resize(size, 1);
1589 boundary_vis_primitive_ids.resize(size, 1);
1590 boundary_vis_normals.resize(size, vertices.front().cols());
1591
1592 if (mesh.is_volume())
1593 boundary_vis_elements.resize(tris.size(), 3);
1594 else
1595 boundary_vis_elements.resize(edges.size(), 2);
1596
1597 int index = 0;
1598 int ii = 0;
1599 for (const auto &v : vertices)
1600 {
1601 boundary_vis_vertices.block(index, 0, v.rows(), v.cols()) = v;
1602 boundary_vis_local_vertices.block(index, 0, v.rows(), v.cols()) = lv[ii];
1603 boundary_vis_elements_ids.block(index, 0, v.rows(), 1).setConstant(el_ids[ii]);
1604 boundary_vis_primitive_ids.block(index, 0, v.rows(), 1).setConstant(global_primitive_ids[ii++]);
1605 index += v.rows();
1606 }
1607
1608 index = 0;
1609 for (const auto &n : allnormals)
1610 {
1611 boundary_vis_normals.block(index, 0, n.rows(), n.cols()) = n;
1612 index += n.rows();
1613 }
1614
1615 index = 0;
1616 if (mesh.is_volume())
1617 {
1618 for (const auto &t : tris)
1619 {
1620 boundary_vis_elements.row(index) << std::get<0>(t), std::get<1>(t), std::get<2>(t);
1621 ++index;
1622 }
1623 }
1624 else
1625 {
1626 for (const auto &e : edges)
1627 {
1628 boundary_vis_elements.row(index) << e.first, e.second;
1629 ++index;
1630 }
1631 }
1632 }
1633
1634 void OutGeometryData::build_vis_mesh(
1635 const mesh::Mesh &mesh,
1636 const Eigen::VectorXi &disc_orders,
1637 const std::vector<basis::ElementBases> &gbases,
1638 const std::map<int, Eigen::MatrixXd> &polys,
1639 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
1640 const bool boundary_only,
1641 Eigen::MatrixXd &points,
1642 Eigen::MatrixXi &tets,
1643 Eigen::MatrixXi &el_id,
1644 Eigen::MatrixXd &discr,
1645 Eigen::MatrixXd &local_points) const
1646 {
1647 const auto &sampler = ref_element_sampler;
1648
1649 const auto &current_bases = gbases;
1650 int tet_total_size = 0;
1651 int pts_total_size = 0;
1652
1653 Eigen::MatrixXd vis_pts_poly;
1654 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
1655
1656 for (size_t i = 0; i < current_bases.size(); ++i)
1657 {
1658 const auto &bs = current_bases[i];
1659
1660 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
1661 continue;
1662
1663 if (mesh.is_simplex(i))
1664 {
1665 tet_total_size += sampler.simplex_volume().rows();
1666 pts_total_size += sampler.simplex_points().rows();
1667 }
1668 else if (mesh.is_cube(i))
1669 {
1670 tet_total_size += sampler.cube_volume().rows();
1671 pts_total_size += sampler.cube_points().rows();
1672 }
1673 else if (mesh.is_prism(i))
1674 {
1675 tet_total_size += sampler.prism_volume().rows();
1676 pts_total_size += sampler.prism_points().rows();
1677 }
1678 else if (mesh.is_pyramid(i))
1679 {
1680 tet_total_size += sampler.pyramid_volume().rows();
1681 pts_total_size += sampler.pyramid_points().rows();
1682 }
1683 else
1684 {
1685 if (mesh.is_volume())
1686 {
1687 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
1688
1689 tet_total_size += vis_faces_poly.rows();
1690 pts_total_size += vis_pts_poly.rows();
1691 }
1692 else
1693 {
1694 sampler.sample_polygon(polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
1695
1696 tet_total_size += vis_faces_poly.rows();
1697 pts_total_size += vis_pts_poly.rows();
1698 }
1699 }
1700 }
1701
1702 points.resize(pts_total_size, mesh.dimension());
1703 local_points.resize(pts_total_size, mesh.dimension());
1704 local_points.setZero();
1705 tets.resize(tet_total_size, mesh.is_volume() ? 4 : 3);
1706
1707 el_id.resize(pts_total_size, 1);
1708 discr.resize(pts_total_size, 1);
1709
1710 Eigen::MatrixXd mapped, tmp;
1711 int tet_index = 0, pts_index = 0;
1712
1713 for (size_t i = 0; i < current_bases.size(); ++i)
1714 {
1715 const auto &bs = current_bases[i];
1716
1717 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
1718 continue;
1719
1720 if (mesh.is_simplex(i))
1721 {
1722 bs.eval_geom_mapping(sampler.simplex_points(), mapped);
1723
1724 tets.block(tet_index, 0, sampler.simplex_volume().rows(), tets.cols()) = sampler.simplex_volume().array() + pts_index;
1725 tet_index += sampler.simplex_volume().rows();
1726
1727 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
1728 local_points.block(pts_index, 0, sampler.simplex_points().rows(), sampler.simplex_points().cols()) = sampler.simplex_points();
1729 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
1730 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
1731 pts_index += mapped.rows();
1732 }
1733 else if (mesh.is_cube(i))
1734 {
1735 bs.eval_geom_mapping(sampler.cube_points(), mapped);
1736
1737 tets.block(tet_index, 0, sampler.cube_volume().rows(), tets.cols()) = sampler.cube_volume().array() + pts_index;
1738 tet_index += sampler.cube_volume().rows();
1739
1740 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
1741 local_points.block(pts_index, 0, sampler.cube_points().rows(), sampler.cube_points().cols()) = sampler.cube_points();
1742 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
1743 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
1744 pts_index += mapped.rows();
1745 }
1746 else if (mesh.is_prism(i))
1747 {
1748 bs.eval_geom_mapping(sampler.prism_points(), mapped);
1749
1750 tets.block(tet_index, 0, sampler.prism_volume().rows(), tets.cols()) = sampler.prism_volume().array() + pts_index;
1751 tet_index += sampler.prism_volume().rows();
1752
1753 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
1754 local_points.block(pts_index, 0, sampler.prism_points().rows(), sampler.prism_points().cols()) = sampler.prism_points();
1755 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
1756 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
1757 pts_index += mapped.rows();
1758 }
1759 else if (mesh.is_pyramid(i))
1760 {
1761 bs.eval_geom_mapping(sampler.pyramid_points(), mapped);
1762
1763 tets.block(tet_index, 0, sampler.pyramid_volume().rows(), tets.cols()) = sampler.pyramid_volume().array() + pts_index;
1764 tet_index += sampler.pyramid_volume().rows();
1765
1766 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
1767 local_points.block(pts_index, 0, sampler.pyramid_points().rows(), sampler.pyramid_points().cols()) = sampler.pyramid_points();
1768 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
1769 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
1770 pts_index += mapped.rows();
1771 }
1772 else
1773 {
1774 if (mesh.is_volume())
1775 {
1776 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
1777 bs.eval_geom_mapping(vis_pts_poly, mapped);
1778
1779 tets.block(tet_index, 0, vis_faces_poly.rows(), tets.cols()) = vis_faces_poly.array() + pts_index;
1780 tet_index += vis_faces_poly.rows();
1781
1782 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
1783 local_points.block(pts_index, 0, vis_pts_poly.rows(), vis_pts_poly.cols()) = vis_pts_poly;
1784 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(-1);
1785 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
1786 pts_index += mapped.rows();
1787 }
1788 else
1789 {
1790 sampler.sample_polygon(polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
1791 bs.eval_geom_mapping(vis_pts_poly, mapped);
1792
1793 tets.block(tet_index, 0, vis_faces_poly.rows(), tets.cols()) = vis_faces_poly.array() + pts_index;
1794 tet_index += vis_faces_poly.rows();
1795
1796 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
1797 local_points.block(pts_index, 0, vis_pts_poly.rows(), vis_pts_poly.cols()) = vis_pts_poly;
1798 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(-1);
1799 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
1800 pts_index += mapped.rows();
1801 }
1802 }
1803 }
1804
1805 assert(pts_index == points.rows());
1806 assert(tet_index == tets.rows());
1807 }
1808
1809 void OutGeometryData::build_high_order_vis_mesh(
1810 const mesh::Mesh &mesh,
1811 const Eigen::VectorXi &output_orders,
1812 const std::vector<basis::ElementBases> &bases,
1813 Eigen::MatrixXd &points,
1814 std::vector<CellElement> &elements,
1815 Eigen::MatrixXi &el_id,
1816 Eigen::MatrixXd &discr,
1817 Eigen::MatrixXd &local_points) const
1818 {
1819 // if (!mesh)
1820 // {
1821 // logger().error("Load the mesh first!");
1822 // return;
1823 // }
1824 // if (n_bases <= 0)
1825 // {
1826 // logger().error("Build the bases first!");
1827 // return;
1828 // }
1829 // assert(mesh.is_linear());
1830
1831 std::vector<RowVectorNd> nodes;
1832 int pts_total_size = 0;
1833 elements.resize(bases.size());
1834 Eigen::MatrixXd ref_pts;
1835
1836 for (size_t i = 0; i < bases.size(); ++i)
1837 {
1838 const auto &bs = bases[i];
1839 if (mesh.is_volume())
1840 {
1841 if (mesh.is_simplex(i))
1842 autogen::p_nodes_3d(output_orders(i), ref_pts);
1843 else if (mesh.is_cube(i))
1844 autogen::q_nodes_3d(output_orders(i), ref_pts);
1845 else if (mesh.is_prism(i))
1846 {
1847 autogen::prism_nodes_3d(output_orders(i), output_orders(i), ref_pts);
1848 }
1849 else if (mesh.is_pyramid(i))
1850 {
1851 if (output_orders(i) == 1)
1852 pyramid_nodes_for_output(1, ref_pts);
1853 else
1854 ref_pts = ref_element_sampler.pyramid_points();
1855 }
1856 else
1857 continue;
1858 }
1859 else
1860 {
1861 if (mesh.is_simplex(i))
1862 autogen::p_nodes_2d(output_orders(i), ref_pts);
1863 else if (mesh.is_cube(i))
1864 autogen::q_nodes_2d(output_orders(i), ref_pts);
1865 else
1866 {
1867 const int n_v = static_cast<const mesh::Mesh2D &>(mesh).n_face_vertices(i);
1868 ref_pts.resize(n_v, 2);
1869 }
1870 }
1871
1872 pts_total_size += ref_pts.rows();
1873 }
1874
1875 points.resize(pts_total_size, mesh.dimension());
1876 local_points.resize(pts_total_size, mesh.dimension());
1877 local_points.setZero();
1878
1879 el_id.resize(pts_total_size, 1);
1880 discr.resize(pts_total_size, 1);
1881
1882 Eigen::MatrixXd mapped;
1883 int pts_index = 0;
1884
1885 std::string error_msg = "";
1886
1887 for (size_t i = 0; i < bases.size(); ++i)
1888 {
1889 const auto &bs = bases[i];
1890 if (mesh.is_volume())
1891 {
1892 if (mesh.is_simplex(i))
1893 autogen::p_nodes_3d(output_orders(i), ref_pts);
1894 else if (mesh.is_cube(i))
1895 autogen::q_nodes_3d(output_orders(i), ref_pts);
1896 else if (mesh.is_prism(i))
1897 {
1898 autogen::prism_nodes_3d(output_orders(i), output_orders(i), ref_pts);
1899 }
1900 else if (mesh.is_pyramid(i))
1901 {
1902 if (output_orders(i) == 1)
1903 pyramid_nodes_for_output(1, ref_pts);
1904 else
1905 ref_pts = ref_element_sampler.pyramid_points();
1906 }
1907 else
1908 continue;
1909 }
1910 else
1911 {
1912 if (mesh.is_simplex(i))
1913 autogen::p_nodes_2d(output_orders(i), ref_pts);
1914 else if (mesh.is_cube(i))
1915 autogen::q_nodes_2d(output_orders(i), ref_pts);
1916 else
1917 continue;
1918 }
1919
1920 bs.eval_geom_mapping(ref_pts, mapped);
1921
1922 for (int j = 0; j < mapped.rows(); ++j)
1923 {
1924 points.row(pts_index) = mapped.row(j);
1925 local_points.row(pts_index).leftCols(ref_pts.cols()) = ref_pts.row(j);
1926 el_id(pts_index) = i;
1927 discr(pts_index) = output_orders(i);
1928 elements[i].vertices.push_back(pts_index);
1929
1930 pts_index++;
1931 }
1932
1933 if (mesh.is_simplex(i))
1934 {
1935 if (mesh.is_volume())
1936 {
1937 const int n_nodes = elements[i].vertices.size();
1938 if (output_orders(i) >= 3)
1939 {
1940 std::swap(elements[i].vertices[16], elements[i].vertices[17]);
1941 std::swap(elements[i].vertices[17], elements[i].vertices[18]);
1942 std::swap(elements[i].vertices[18], elements[i].vertices[19]);
1943 }
1944 if (output_orders(i) > 4)
1945 error_msg = "Saving high-order meshes not implemented for P5+ elements!";
1946 }
1947 else
1948 {
1949 if (output_orders(i) == 4)
1950 {
1951 const int n_nodes = elements[i].vertices.size();
1952 std::swap(elements[i].vertices[n_nodes - 1], elements[i].vertices[n_nodes - 2]);
1953 }
1954 if (output_orders(i) > 4)
1955 error_msg = "Saving high-order meshes not implemented for P5+ elements!";
1956 }
1957 }
1958 else if (mesh.is_cube(i) && mesh.is_volume())
1959 {
1960 const int n_nodes = elements[i].vertices.size();
1961 if (output_orders(i) == 2) // Lagrange hex, order=2
1962 {
1963 std::swap(elements[i].vertices[12], elements[i].vertices[16]);
1964 std::swap(elements[i].vertices[13], elements[i].vertices[17]);
1965 std::swap(elements[i].vertices[14], elements[i].vertices[18]);
1966 std::swap(elements[i].vertices[15], elements[i].vertices[19]);
1967 std::swap(elements[i].vertices[18], elements[i].vertices[19]); // a hack fix
1968 }
1969 // if (disc_orders(i) == 3) // Incomplete fix, need to fix order on the edge
1970 // {
1971 // std::swap(elements[i].vertices[24], elements[i].vertices[16]);
1972 // std::swap(elements[i].vertices[25], elements[i].vertices[17]);
1973 // std::swap(elements[i].vertices[26], elements[i].vertices[18]);
1974 // std::swap(elements[i].vertices[27], elements[i].vertices[19]);
1975 // std::swap(elements[i].vertices[28], elements[i].vertices[20]);
1976 // std::swap(elements[i].vertices[29], elements[i].vertices[21]);
1977 // std::swap(elements[i].vertices[30], elements[i].vertices[22]);
1978 // std::swap(elements[i].vertices[31], elements[i].vertices[23]);
1979 // std::swap(elements[i].vertices[28], elements[i].vertices[30]); // hack
1980 // std::swap(elements[i].vertices[29], elements[i].vertices[31]); // hack
1981 // }
1982 if (output_orders(i) > 2)
1983 error_msg = "Saving high-order meshes not implemented for Q2+ elements!";
1984 }
1985 else if (output_orders(i) > 1)
1986 {
1987 if (mesh.is_cube(i))
1988 error_msg = "Saving high-order meshes not implemented for Q2+ elements!";
1989 }
1990 }
1991
1992 if (!error_msg.empty())
1993 logger().warn(error_msg);
1994
1995 for (size_t i = 0; i < bases.size(); ++i)
1996 {
1997 if (mesh.is_volume() || !mesh.is_polytope(i))
1998 continue;
1999
2000 const auto &mesh2d = static_cast<const mesh::Mesh2D &>(mesh);
2001 const int n_v = mesh2d.n_face_vertices(i);
2002
2003 for (int j = 0; j < n_v; ++j)
2004 {
2005 points.row(pts_index) = mesh2d.point(mesh2d.face_vertex(i, j));
2006 local_points.row(pts_index) = mesh2d.point(mesh2d.face_vertex(i, j));
2007 el_id(pts_index) = i;
2008 discr(pts_index) = output_orders(i);
2009 elements[i].vertices.push_back(pts_index);
2010
2011 pts_index++;
2012 }
2013 }
2014
2015 for (size_t i = 0; i < bases.size(); ++i)
2016 {
2017 if (!mesh.is_volume())
2018 {
2019 if (elements[i].vertices.size() == 1)
2020 elements[i].ctype = CellType::Vertex;
2021 else if (elements[i].vertices.size() == 2)
2022 elements[i].ctype = CellType::Line;
2023 else if (mesh.is_simplex(i))
2024 elements[i].ctype = CellType::Triangle;
2025 else if (mesh.is_cube(i))
2026 elements[i].ctype = CellType::Quadrilateral;
2027 else
2028 elements[i].ctype = CellType::Polygon;
2029 }
2030 else
2031 {
2032 if (mesh.is_simplex(i))
2033 elements[i].ctype = CellType::Tetrahedron;
2034 else if (mesh.is_cube(i))
2035 elements[i].ctype = CellType::Hexahedron;
2036 else if (mesh.is_prism(i))
2037 elements[i].ctype = CellType::Wedge;
2038 else if (mesh.is_pyramid(i))
2039 elements[i].ctype = CellType::Pyramid;
2040 }
2041 }
2042
2043 if (mesh.is_volume())
2044 {
2045 // ParaView does not reliably render high-order Lagrange pyramids;
2046 // tessellate only those cells while keeping linear pyramids and the
2047 // other element types in their native representation.
2048 std::vector<CellElement> expanded_elements;
2049 expanded_elements.reserve(elements.size());
2050 for (size_t i = 0; i < bases.size(); ++i)
2051 {
2052 if (!mesh.is_pyramid(i) || output_orders(i) == 1)
2053 {
2054 expanded_elements.push_back(std::move(elements[i]));
2055 continue;
2056 }
2057
2058 for (int t = 0; t < ref_element_sampler.pyramid_volume().rows(); ++t)
2059 {
2060 CellElement tet;
2061 tet.ctype = CellType::Tetrahedron;
2062 for (int j = 0; j < ref_element_sampler.pyramid_volume().cols(); ++j)
2063 tet.vertices.push_back(elements[i].vertices[ref_element_sampler.pyramid_volume()(t, j)]);
2064 expanded_elements.push_back(std::move(tet));
2065 }
2066 }
2067 elements.swap(expanded_elements);
2068 }
2069
2070 assert(pts_index == points.rows());
2071 }
2072
2073 void OutGeometryData::export_data(
2074 const OutputSpace &space,
2075 const OutputFieldFunction &output_fields,
2076 const bool is_time_dependent,
2077 const double tend_in,
2078 const double dt,
2079 const ExportOptions &opts,
2080 const std::string &vis_mesh_path) const
2081 {
2082 if (!space.mesh)
2083 {
2084 logger().error("Load the mesh first!");
2085 return;
2086 }
2087
2088 double tend = tend_in;
2089 if (tend <= 0)
2090 tend = 1;
2091
2092 if (!vis_mesh_path.empty() && !is_time_dependent)
2093 {
2094 save_vtu(
2095 vis_mesh_path, space, output_fields,
2096 tend, dt, opts);
2097 }
2098 }
2099
2100 bool OutGeometryData::ExportOptions::export_field(const std::string &field) const
2101 {
2102 return fields.empty() || std::find(fields.begin(), fields.end(), field) != fields.end();
2103 }
2104
2105 OutGeometryData::ExportOptions::ExportOptions(const json &args, const bool is_mesh_linear, const bool mesh_has_prisms, const bool is_problem_scalar)
2106 {
2107 fields = args["output"]["paraview"]["fields"];
2108
2109 volume = args["output"]["paraview"]["volume"];
2110 surface = args["output"]["paraview"]["surface"];
2111 wire = args["output"]["paraview"]["wireframe"];
2112 points = args["output"]["paraview"]["points"];
2113 contact_forces = args["output"]["paraview"]["options"]["contact_forces"] && !is_problem_scalar;
2114 friction_forces = args["output"]["paraview"]["options"]["friction_forces"] && !is_problem_scalar;
2115 normal_adhesion_forces = args["output"]["paraview"]["options"]["normal_adhesion_forces"] && !is_problem_scalar;
2116 tangential_adhesion_forces = args["output"]["paraview"]["options"]["tangential_adhesion_forces"] && !is_problem_scalar;
2117
2118 if (args["output"]["paraview"]["options"]["force_high_order"])
2119 use_sampler = false;
2120 else
2121 use_sampler = !(is_mesh_linear && args["output"]["paraview"]["high_order_mesh"]);
2122 boundary_only = use_sampler && args["output"]["advanced"]["vis_boundary_only"];
2123 sol_on_grid = args["output"]["advanced"]["sol_on_grid"] > 0;
2124
2125 discretization_order = args["output"]["paraview"]["options"]["discretization_order"];
2126
2127 reorder_output = args["output"]["data"]["advanced"]["reorder_nodes"];
2128
2129 use_hdf5 = args["output"]["paraview"]["options"]["use_hdf5"];
2130 }
2131
2132 void OutGeometryData::save_vtu(
2133 const std::string &path,
2134 const OutputSpace &space,
2135 const OutputFieldFunction &output_fields,
2136 const double t,
2137 const double dt,
2138 const ExportOptions &opts) const
2139 {
2140 if (!space.mesh)
2141 {
2142 logger().error("Load the mesh first!");
2143 return;
2144 }
2145
2146 const std::filesystem::path fs_path(path);
2147 const std::string path_stem = fs_path.stem().string();
2148 const std::string base_path = (fs_path.parent_path() / path_stem).string();
2149 paraviewo::VTMWriter vtm(t);
2150 save_vtu(path, space, output_fields, t, dt, opts, vtm, "");
2151 vtm.save(base_path + ".vtm");
2152 }
2153
2154 void OutGeometryData::save_vtu(
2155 const std::string &path,
2156 const OutputSpace &space,
2157 const OutputFieldFunction &output_fields,
2158 const double t,
2159 const double dt,
2160 const ExportOptions &opts,
2161 paraviewo::VTMWriter &vtm,
2162 const std::string &block_prefix) const
2163 {
2164 if (!space.mesh)
2165 {
2166 logger().error("Load the mesh first!");
2167 return;
2168 }
2169
2170 const bool save_contact =
2171 space.collision_mesh
2172 && (opts.contact_forces || opts.friction_forces || opts.normal_adhesion_forces || opts.tangential_adhesion_forces
2173 || (!opts.fields.empty() && opts.export_field("adaptive_dhat")));
2174
2175 logger().info("Saving vtu to {}; volume={}, surface={}, contact={}, points={}, wireframe={}",
2176 path, opts.volume, opts.surface, save_contact, opts.points, opts.wire);
2177
2178 const std::filesystem::path fs_path(path);
2179 const std::string path_stem = fs_path.stem().string();
2180 const std::string base_path = (fs_path.parent_path() / path_stem).string();
2181
2182 if (opts.volume)
2183 {
2184 save_volume(base_path + opts.file_extension(), space, output_fields, t, dt, opts);
2185 }
2186
2187 if (opts.surface)
2188 {
2189 save_surface(base_path + "_surf" + opts.file_extension(), space, output_fields, t, dt, opts);
2190 }
2191
2192 if (save_contact)
2193 {
2194 save_contact_surface(base_path + "_surf" + opts.file_extension(), space, output_fields, t, dt, opts);
2195 }
2196
2197 if (opts.wire)
2198 {
2199 save_wire(base_path + "_wire" + opts.file_extension(), space, output_fields, t, opts);
2200 }
2201
2202 if (opts.points)
2203 {
2204 save_points(base_path + "_points" + opts.file_extension(), space, output_fields, opts);
2205 }
2206
2207 const auto block_name = [&block_prefix](const std::string &name) {
2208 return block_prefix.empty() ? name : block_prefix + " " + name;
2209 };
2210 if (opts.volume)
2211 vtm.add_dataset(block_name("Volume"), "data", path_stem + opts.file_extension());
2212 if (opts.surface)
2213 vtm.add_dataset(block_name("Surface"), "data", path_stem + "_surf" + opts.file_extension());
2214 if (save_contact)
2215 vtm.add_dataset(block_name("Contact"), "data", path_stem + "_surf_contact" + opts.file_extension());
2216 if (opts.wire)
2217 vtm.add_dataset(block_name("Wireframe"), "data", path_stem + "_wire" + opts.file_extension());
2218 if (opts.points)
2219 vtm.add_dataset(block_name("Points"), "data", path_stem + "_points" + opts.file_extension());
2220 }
2221
2222 void OutGeometryData::save_volume(
2223 const std::string &path,
2224 const OutputSpace &space,
2225 const OutputFieldFunction &output_fields,
2226 const double t,
2227 const double dt,
2228 const ExportOptions &opts) const
2229 {
2230 if (!space.mesh || !space.geometry_bases)
2231 return;
2232
2233 static const std::map<int, Eigen::MatrixXd> empty_polys;
2234 static const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> empty_polys_3d;
2235
2236 const mesh::Mesh &mesh = *space.mesh;
2237 const std::vector<basis::ElementBases> &gbases = *space.geometry_bases;
2238 const std::map<int, Eigen::MatrixXd> &polys = space.polys ? *space.polys : empty_polys;
2239 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d = space.polys_3d ? *space.polys_3d : empty_polys_3d;
2240 const Eigen::VectorXi output_orders =
2241 space.output_orders.size() == mesh.n_elements()
2242 ? space.output_orders
2243 : Eigen::VectorXi::Ones(mesh.n_elements());
2244 const mesh::Obstacle *obstacle = space.obstacle;
2245
2246 Eigen::MatrixXd points;
2247 Eigen::MatrixXi tets;
2248 Eigen::MatrixXi el_id;
2249 Eigen::MatrixXd discr;
2250 Eigen::MatrixXd local_points;
2251 std::vector<CellElement> elements;
2252
2253 if (opts.use_sampler)
2254 build_vis_mesh(mesh, output_orders, gbases,
2255 polys, polys_3d, opts.boundary_only,
2256 points, tets, el_id, discr, local_points);
2257 else
2258 {
2259 build_high_order_vis_mesh(mesh, output_orders, gbases,
2260 points, elements, el_id, discr, local_points);
2261 }
2262
2263 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2264 if (opts.use_hdf5)
2265 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2266 else
2267 tmpw = std::make_shared<paraviewo::VTUWriter>();
2268 paraviewo::ParaviewWriter &writer = *tmpw;
2269
2270 if (obstacle && obstacle->n_vertices() > 0)
2271 {
2272 discr.conservativeResize(discr.size() + obstacle->n_vertices(), 1);
2273 discr.bottomRows(obstacle->n_vertices()).setZero();
2274 }
2275
2276 if (opts.discretization_order && opts.export_field("discr"))
2277 writer.add_field("discr", discr);
2278
2279 if (obstacle && obstacle->n_vertices() > 0)
2280 {
2281 const int orig_p = points.rows();
2282 points.conservativeResize(points.rows() + obstacle->n_vertices(), points.cols());
2283 points.bottomRows(obstacle->n_vertices()) = obstacle->v();
2284
2285 if (elements.empty())
2286 {
2287 for (int i = 0; i < tets.rows(); ++i)
2288 {
2289 elements.emplace_back();
2290 elements.back().ctype = CellType::Tetrahedron;
2291 for (int j = 0; j < tets.cols(); ++j)
2292 elements.back().vertices.push_back(tets(i, j));
2293 }
2294 }
2295
2296 for (int i = 0; i < obstacle->get_face_connectivity().rows(); ++i)
2297 {
2298 elements.emplace_back();
2299 elements.back().ctype = CellType::Triangle;
2300 for (int j = 0; j < obstacle->get_face_connectivity().cols(); ++j)
2301 elements.back().vertices.push_back(obstacle->get_face_connectivity()(i, j) + orig_p);
2302 }
2303
2304 for (int i = 0; i < obstacle->get_edge_connectivity().rows(); ++i)
2305 {
2306 elements.emplace_back();
2307 elements.back().ctype = CellType::Line;
2308 for (int j = 0; j < obstacle->get_edge_connectivity().cols(); ++j)
2309 elements.back().vertices.push_back(obstacle->get_edge_connectivity()(i, j) + orig_p);
2310 }
2311
2312 for (int i = 0; i < obstacle->get_vertex_connectivity().size(); ++i)
2313 {
2314 elements.emplace_back();
2315 elements.back().ctype = CellType::Vertex;
2316 elements.back().vertices.push_back(obstacle->get_vertex_connectivity()(i) + orig_p);
2317 }
2318 }
2319
2320 // Write the solution alias last so it is the default for warp-by-vector.
2321 OutputSample sample;
2322 sample.points = points;
2323 sample.local_points = local_points;
2324 sample.element_ids = el_id.col(0);
2325 sample.domain = OutputSample::Domain::Volume;
2326 sample.cell_count = elements.empty() ? tets.rows() : static_cast<int>(elements.size());
2327 sample.time = t;
2328 sample.dt = dt;
2329 add_output_fields(writer, sample, output_fields);
2330
2331 if (opts.sol_on_grid && output_fields && grid_points.rows() > 0)
2332 {
2333 OutputSample grid_sample;
2334 grid_sample.points = grid_points;
2335 grid_sample.element_ids = grid_points_to_elements.col(0);
2336 grid_sample.domain = OutputSample::Domain::Grid;
2337 grid_sample.local_points.resize(grid_points.rows(), mesh.dimension());
2338 grid_sample.local_points.setZero();
2339 for (int i = 0; i < grid_points.rows(); ++i)
2340 {
2341 if (grid_sample.element_ids(i) >= 0)
2342 grid_sample.local_points.row(i) = grid_points_bc.row(i).rightCols(mesh.dimension());
2343 }
2344 grid_sample.time = t;
2345 grid_sample.dt = dt;
2346 grid_sample.requested_fields = {
2347 "solution",
2348 "solution_gradient",
2349 "pressure",
2350 "pressure_gradient",
2351 };
2352
2353 io::write_matrix(path + "_grid.txt", grid_points);
2354 for (const OutputField &field : output_fields(grid_sample))
2355 {
2356 if (field.association != OutputField::Association::Point || field.values.rows() != grid_points.rows())
2357 continue;
2358 if (field.name == "solution")
2359 io::write_matrix(path + "_sol.txt", field.values);
2360 else if (field.name == "solution_gradient")
2361 io::write_matrix(path + "_grad.txt", field.values);
2362 else if (field.name == "pressure")
2363 io::write_matrix(path + "_p_sol.txt", field.values);
2364 else if (field.name == "pressure_gradient")
2365 io::write_matrix(path + "_p_grad.txt", field.values);
2366 }
2367 }
2368
2369 if (elements.empty())
2370 writer.write_mesh(path, points, tets, mesh.is_volume() ? CellType::Tetrahedron : CellType::Triangle);
2371 else
2372 writer.write_mesh(path, points, elements);
2373 }
2374
2375 void OutGeometryData::save_surface(
2376 const std::string &export_surface,
2377 const OutputSpace &space,
2378 const OutputFieldFunction &output_fields,
2379 const double t,
2380 const double dt_in,
2381 const ExportOptions &opts) const
2382 {
2383 if (!space.mesh || !space.geometry_bases || !space.total_local_boundary)
2384 return;
2385
2386 const mesh::Mesh &mesh = *space.mesh;
2387 const std::vector<basis::ElementBases> &gbases = *space.geometry_bases;
2388
2389 Eigen::MatrixXd boundary_vis_vertices;
2390 Eigen::MatrixXd boundary_vis_local_vertices;
2391 Eigen::MatrixXi boundary_vis_elements;
2392 Eigen::MatrixXi boundary_vis_elements_ids;
2393 Eigen::MatrixXi boundary_vis_primitive_ids;
2394 Eigen::MatrixXd boundary_vis_normals;
2395
2396 build_vis_boundary_mesh(mesh, gbases, *space.total_local_boundary,
2397 boundary_vis_vertices, boundary_vis_local_vertices, boundary_vis_elements,
2398 boundary_vis_elements_ids, boundary_vis_primitive_ids, boundary_vis_normals);
2399
2400 Eigen::MatrixXd discr, b_sidesets;
2401 discr.resize(boundary_vis_vertices.rows(), 1);
2402 b_sidesets.resize(boundary_vis_vertices.rows(), 1);
2403 b_sidesets.setZero();
2404
2405 for (int i = 0; i < boundary_vis_vertices.rows(); ++i)
2406 {
2407 const auto s_id = mesh.get_boundary_id(boundary_vis_primitive_ids(i));
2408 if (s_id > 0)
2409 {
2410 b_sidesets(i) = s_id;
2411 }
2412
2413 const int el_index = boundary_vis_elements_ids(i);
2414 discr(i) = space.output_orders.size() == mesh.n_elements() ? space.output_orders(el_index) : 1;
2415 }
2416
2417 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2418 if (opts.use_hdf5)
2419 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2420 else
2421 tmpw = std::make_shared<paraviewo::VTUWriter>();
2422 paraviewo::ParaviewWriter &writer = *tmpw;
2423
2424 if (opts.export_field("normals"))
2425 writer.add_field("normals", boundary_vis_normals);
2426 if (opts.export_field("discr"))
2427 writer.add_field("discr", discr);
2428 if (opts.export_field("sidesets"))
2429 writer.add_field("sidesets", b_sidesets);
2430
2431 // Write the solution alias last so it is the default for warp-by-vector.
2432 OutputSample sample;
2433 sample.points = boundary_vis_vertices;
2434 sample.local_points = boundary_vis_local_vertices;
2435 sample.element_ids = boundary_vis_elements_ids.col(0);
2436 sample.primitive_ids = boundary_vis_primitive_ids;
2437 sample.normals = boundary_vis_normals;
2438 sample.domain = OutputSample::Domain::Surface;
2439 sample.cell_count = boundary_vis_elements.rows();
2440 sample.time = t;
2441 sample.dt = dt_in;
2442 add_output_fields(writer, sample, output_fields);
2443 writer.write_mesh(export_surface, boundary_vis_vertices, boundary_vis_elements, mesh.is_volume() ? CellType::Triangle : CellType::Line);
2444 }
2445
2446 void OutGeometryData::save_contact_surface(
2447 const std::string &export_surface,
2448 const OutputSpace &space,
2449 const OutputFieldFunction &output_fields,
2450 const double t,
2451 const double dt_in,
2452 const ExportOptions &opts) const
2453 {
2454 if (!space.collision_mesh)
2455 return;
2456
2457 const ipc::CollisionMesh &collision_mesh = *space.collision_mesh;
2458
2459 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2460 if (opts.use_hdf5)
2461 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2462 else
2463 tmpw = std::make_shared<paraviewo::VTUWriter>();
2464 paraviewo::ParaviewWriter &writer = *tmpw;
2465
2466 // Write the solution alias last so it is the default for warp-by-vector.
2467 OutputSample sample;
2468 sample.points = collision_mesh.rest_positions();
2469 sample.domain = OutputSample::Domain::Contact;
2470 sample.cell_count = static_cast<int>(
2471 collision_mesh.dim() == 3 ? collision_mesh.num_faces() : collision_mesh.num_edges());
2472 sample.time = t;
2473 sample.dt = dt_in;
2474 add_output_fields(writer, sample, output_fields);
2475
2476 const std::filesystem::path surface_path(export_surface);
2477 const std::string contact_path =
2478 (surface_path.parent_path() / (surface_path.stem().string() + "_contact" + surface_path.extension().string())).string();
2479 writer.write_mesh(
2480 contact_path,
2481 collision_mesh.rest_positions(),
2482 collision_mesh.dim() == 3 ? collision_mesh.faces() : collision_mesh.edges(),
2483 collision_mesh.dim() == 3 ? CellType::Triangle : CellType::Line);
2484 }
2485
2486 void OutGeometryData::save_wire(
2487 const std::string &name,
2488 const OutputSpace &space,
2489 const OutputFieldFunction &output_fields,
2490 const double t,
2491 const ExportOptions &opts) const
2492 {
2493 if (!space.mesh || !space.geometry_bases)
2494 return;
2495
2496 static const std::map<int, Eigen::MatrixXd> empty_polys;
2497 static const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> empty_polys_3d;
2498
2499 const std::vector<basis::ElementBases> &gbases = *space.geometry_bases;
2500 const mesh::Mesh &mesh = *space.mesh;
2501 const std::map<int, Eigen::MatrixXd> &polys = space.polys ? *space.polys : empty_polys;
2502 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d = space.polys_3d ? *space.polys_3d : empty_polys_3d;
2503 const Eigen::VectorXi output_orders =
2504 space.output_orders.size() == mesh.n_elements()
2505 ? space.output_orders
2506 : Eigen::VectorXi::Ones(mesh.n_elements());
2507
2508 Eigen::MatrixXd points, discr, local_points;
2509 Eigen::MatrixXi cells, element_ids, edges;
2510 build_vis_mesh(mesh, output_orders, gbases, polys, polys_3d, /*boundary_only=*/false, points, cells, element_ids, discr, local_points);
2511 if (cells.size() > 0)
2512 igl::edges(cells, edges);
2513 else
2514 edges.resize(0, 2);
2515
2516 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2517 if (opts.use_hdf5)
2518 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2519 else
2520 tmpw = std::make_shared<paraviewo::VTUWriter>();
2521 paraviewo::ParaviewWriter &writer = *tmpw;
2522
2523 // Write the solution alias last so it is the default for warp-by-vector.
2524 OutputSample sample;
2525 sample.points = points;
2526 sample.local_points = local_points;
2527 if (element_ids.cols() > 0)
2528 sample.element_ids = element_ids.col(0);
2529 sample.domain = OutputSample::Domain::Wire;
2530 sample.cell_count = edges.rows();
2531 sample.time = t;
2532 add_output_fields(writer, sample, output_fields);
2533
2534 writer.write_mesh(name, points, edges, CellType::Line);
2535 }
2536
2537 void OutGeometryData::save_points(
2538 const std::string &path,
2539 const OutputSpace &space,
2540 const OutputFieldFunction &output_fields,
2541 const ExportOptions &opts) const
2542 {
2543 if (!space.mesh || !space.dirichlet_nodes || !space.dirichlet_nodes_position)
2544 return;
2545
2546 const auto &dirichlet_nodes = *space.dirichlet_nodes;
2547 const auto &dirichlet_nodes_position = *space.dirichlet_nodes_position;
2548 const mesh::Mesh &mesh = *space.mesh;
2549
2550 Eigen::MatrixXd b_sidesets(dirichlet_nodes_position.size(), 1);
2551 b_sidesets.setZero();
2552 Eigen::MatrixXd points(dirichlet_nodes_position.size(), mesh.dimension());
2553 std::vector<CellElement> cells(dirichlet_nodes_position.size());
2554
2555 for (int i = 0; i < dirichlet_nodes_position.size(); ++i)
2556 {
2557 const int n_id = dirichlet_nodes[i];
2558 const auto s_id = mesh.get_node_id(n_id);
2559 if (s_id > 0)
2560 {
2561 b_sidesets(i) = s_id;
2562 }
2563
2564 points.row(i) = dirichlet_nodes_position[i];
2565 cells[i].vertices.push_back(i);
2566 cells[i].ctype = CellType::Vertex;
2567 }
2568
2569 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2570 if (opts.use_hdf5)
2571 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2572 else
2573 tmpw = std::make_shared<paraviewo::VTUWriter>();
2574 paraviewo::ParaviewWriter &writer = *tmpw;
2575
2576 if (opts.export_field("sidesets"))
2577 writer.add_field("sidesets", b_sidesets);
2578
2579 // Write the solution alias last so it is the default for warp-by-vector.
2580 OutputSample sample;
2581 sample.points = points;
2582 sample.node_ids.resize(dirichlet_nodes.size());
2583 for (int i = 0; i < dirichlet_nodes.size(); ++i)
2584 sample.node_ids(i) = dirichlet_nodes[i];
2585 sample.domain = OutputSample::Domain::Points;
2586 sample.cell_count = static_cast<int>(cells.size());
2587 add_output_fields(writer, sample, output_fields);
2588 writer.write_mesh(path, points, cells);
2589 }
2590
2591 void OutGeometryData::save_pvd(
2592 const std::string &name,
2593 const std::function<std::string(int)> &vtu_names,
2594 int time_steps, double t0, double dt, int skip_frame) const
2595 {
2596 paraviewo::PVDWriter::save_pvd(name, vtu_names, time_steps, t0, dt, skip_frame);
2597 }
2598
2599 void OutGeometryData::init_sampler(const polyfem::mesh::Mesh &mesh, const double vismesh_rel_area)
2600 {
2601 ref_element_sampler.init(mesh.is_volume(), mesh.n_elements(), vismesh_rel_area);
2602 }
2603
2604 void OutGeometryData::build_grid(const polyfem::mesh::Mesh &mesh, const double spacing)
2605 {
2606 if (spacing <= 0)
2607 return;
2608
2609 RowVectorNd min, max;
2610 mesh.bounding_box(min, max);
2611 const RowVectorNd delta = max - min;
2612 const int nx = delta[0] / spacing + 1;
2613 const int ny = delta[1] / spacing + 1;
2614 const int nz = delta.cols() >= 3 ? (delta[2] / spacing + 1) : 1;
2615 const int n = nx * ny * nz;
2616
2617 grid_points.resize(n, delta.cols());
2618 int index = 0;
2619 for (int i = 0; i < nx; ++i)
2620 {
2621 const double x = (delta[0] / (nx - 1)) * i + min[0];
2622
2623 for (int j = 0; j < ny; ++j)
2624 {
2625 const double y = (delta[1] / (ny - 1)) * j + min[1];
2626
2627 if (delta.cols() <= 2)
2628 {
2629 grid_points.row(index++) << x, y;
2630 }
2631 else
2632 {
2633 for (int k = 0; k < nz; ++k)
2634 {
2635 const double z = (delta[2] / (nz - 1)) * k + min[2];
2636 grid_points.row(index++) << x, y, z;
2637 }
2638 }
2639 }
2640 }
2641
2642 assert(index == n);
2643
2644 std::vector<std::array<Eigen::Vector3d, 2>> boxes;
2645 mesh.elements_boxes(boxes);
2646
2647 SimpleBVH::BVH bvh;
2648 bvh.init(boxes);
2649
2650 const double eps = 1e-6;
2651
2652 grid_points_to_elements.resize(grid_points.rows(), 1);
2653 grid_points_to_elements.setConstant(-1);
2654
2655 grid_points_bc.resize(grid_points.rows(), mesh.is_volume() ? 4 : 3);
2656
2657 for (int i = 0; i < grid_points.rows(); ++i)
2658 {
2659 const Eigen::Vector3d min(
2660 grid_points(i, 0) - eps,
2661 grid_points(i, 1) - eps,
2662 (mesh.is_volume() ? grid_points(i, 2) : 0) - eps);
2663
2664 const Eigen::Vector3d max(
2665 grid_points(i, 0) + eps,
2666 grid_points(i, 1) + eps,
2667 (mesh.is_volume() ? grid_points(i, 2) : 0) + eps);
2668
2669 std::vector<unsigned int> candidates;
2670
2671 bvh.intersect_box(min, max, candidates);
2672
2673 for (const auto cand : candidates)
2674 {
2675 if (!mesh.is_simplex(cand))
2676 {
2677 logger().warn("Element {} is not simplex, skipping", cand);
2678 continue;
2679 }
2680
2681 Eigen::MatrixXd coords;
2682 mesh.barycentric_coords(grid_points.row(i), cand, coords);
2683
2684 for (int d = 0; d < coords.size(); ++d)
2685 {
2686 if (fabs(coords(d)) < 1e-8)
2687 coords(d) = 0;
2688 else if (fabs(coords(d) - 1) < 1e-8)
2689 coords(d) = 1;
2690 }
2691
2692 if (coords.array().minCoeff() >= 0 && coords.array().maxCoeff() <= 1)
2693 {
2694 grid_points_to_elements(i) = cand;
2695 grid_points_bc.row(i) = coords;
2696 break;
2697 }
2698 }
2699 }
2700 }
2701
2702 void OutStatsData::compute_mesh_size(const polyfem::mesh::Mesh &mesh_in, const std::vector<polyfem::basis::ElementBases> &bases_in, const int n_samples, const bool use_curved_mesh_size)
2703 {
2704 Eigen::MatrixXd samples_simplex, samples_cube, mapped, p0, p1, p;
2705
2706 mesh_size = 0;
2707 average_edge_length = 0;
2708 min_edge_length = std::numeric_limits<double>::max();
2709
2710 if (!use_curved_mesh_size)
2711 {
2712 mesh_in.get_edges(p0, p1);
2713 p = p0 - p1;
2714 min_edge_length = p.rowwise().norm().minCoeff();
2715 average_edge_length = p.rowwise().norm().mean();
2716 mesh_size = p.rowwise().norm().maxCoeff();
2717
2718 logger().info("hmin: {}", min_edge_length);
2719 logger().info("hmax: {}", mesh_size);
2720 logger().info("havg: {}", average_edge_length);
2721
2722 return;
2723 }
2724
2725 if (mesh_in.is_volume())
2726 {
2727 utils::EdgeSampler::sample_3d_simplex(n_samples, samples_simplex);
2728 utils::EdgeSampler::sample_3d_cube(n_samples, samples_cube);
2729 }
2730 else
2731 {
2732 utils::EdgeSampler::sample_2d_simplex(n_samples, samples_simplex);
2733 utils::EdgeSampler::sample_2d_cube(n_samples, samples_cube);
2734 }
2735
2736 int n = 0;
2737 for (size_t i = 0; i < bases_in.size(); ++i)
2738 {
2739 if (mesh_in.is_polytope(i))
2740 continue;
2741 int n_edges;
2742
2743 if (mesh_in.is_simplex(i))
2744 {
2745 n_edges = mesh_in.is_volume() ? 6 : 3;
2746 bases_in[i].eval_geom_mapping(samples_simplex, mapped);
2747 }
2748 else
2749 {
2750 n_edges = mesh_in.is_volume() ? 12 : 4;
2751 bases_in[i].eval_geom_mapping(samples_cube, mapped);
2752 }
2753
2754 for (int j = 0; j < n_edges; ++j)
2755 {
2756 double current_edge = 0;
2757 for (int k = 0; k < n_samples - 1; ++k)
2758 {
2759 p0 = mapped.row(j * n_samples + k);
2760 p1 = mapped.row(j * n_samples + k + 1);
2761 p = p0 - p1;
2762
2763 current_edge += p.norm();
2764 }
2765
2766 mesh_size = std::max(current_edge, mesh_size);
2767 min_edge_length = std::min(current_edge, min_edge_length);
2768 average_edge_length += current_edge;
2769 ++n;
2770 }
2771 }
2772
2773 average_edge_length /= n;
2774
2775 logger().info("hmin: {}", min_edge_length);
2776 logger().info("hmax: {}", mesh_size);
2777 logger().info("havg: {}", average_edge_length);
2778 }
2779
2780 void OutStatsData::reset()
2781 {
2782 *this = OutStatsData();
2783 }
2784
2785 void OutStatsData::count_flipped_elements(const polyfem::mesh::Mesh &mesh, const std::vector<polyfem::basis::ElementBases> &gbases)
2786 {
2787 using namespace mesh;
2788
2789 logger().info("Counting flipped elements...");
2790 const auto &els_tag = mesh.elements_tag();
2791
2792 // flipped_elements.clear();
2793 for (size_t i = 0; i < gbases.size(); ++i)
2794 {
2795 if (mesh.is_polytope(i))
2796 continue;
2797
2798 polyfem::assembler::ElementAssemblyValues vals;
2799 if (!vals.is_geom_mapping_positive(mesh.is_volume(), gbases[i]))
2800 {
2801 ++n_flipped;
2802
2803 static const std::vector<std::string> element_type_names{{
2804 "Simplex",
2805 "RegularInteriorCube",
2806 "RegularBoundaryCube",
2807 "SimpleSingularInteriorCube",
2808 "MultiSingularInteriorCube",
2809 "SimpleSingularBoundaryCube",
2810 "InterfaceCube",
2811 "MultiSingularBoundaryCube",
2812 "BoundaryPolytope",
2813 "InteriorPolytope",
2814 "Undefined",
2815 }};
2816
2817 log_and_throw_error("element {} is flipped, type {}", i, element_type_names[static_cast<int>(els_tag[i])]);
2818 }
2819 }
2820
2821 logger().info(" done");
2822
2823 // dynamic_cast<Mesh3D *>(mesh.get())->save({56}, 1, "mesh.HYBRID");
2824
2825 // std::sort(flipped_elements.begin(), flipped_elements.end());
2826 // auto it = std::unique(flipped_elements.begin(), flipped_elements.end());
2827 // flipped_elements.resize(std::distance(flipped_elements.begin(), it));
2828 }
2829
2830 void OutStatsData::compute_errors(
2831 const int n_bases,
2832 const std::vector<polyfem::basis::ElementBases> &bases,
2833 const std::vector<polyfem::basis::ElementBases> &gbases,
2834 const polyfem::mesh::Mesh &mesh,
2835 const assembler::Problem &problem,
2836 const double tend,
2837 const Eigen::MatrixXd &sol)
2838 {
2839 if (n_bases <= 0)
2840 {
2841 logger().error("Build the bases first!");
2842 return;
2843 }
2844 if (sol.size() <= 0)
2845 {
2846 logger().error("Solve the problem first!");
2847 return;
2848 }
2849
2850 int actual_dim = 1;
2851 if (!problem.is_scalar())
2852 actual_dim = mesh.dimension();
2853
2854 igl::Timer timer;
2855 timer.start();
2856 logger().info("Computing errors...");
2857 using std::max;
2858
2859 const int n_el = int(bases.size());
2860
2861 Eigen::MatrixXd v_exact, v_approx;
2862 Eigen::MatrixXd v_exact_grad(0, 0), v_approx_grad;
2863
2864 l2_err = 0;
2865 h1_err = 0;
2866 grad_max_err = 0;
2867 h1_semi_err = 0;
2868 linf_err = 0;
2869 lp_err = 0;
2870 // double pred_norm = 0;
2871
2872 static const int p = 8;
2873
2874 // Eigen::MatrixXd err_per_el(n_el, 5);
2875 polyfem::assembler::ElementAssemblyValues vals;
2876
2877 for (int e = 0; e < n_el; ++e)
2878 {
2879 vals.compute(e, mesh.is_volume(), bases[e], gbases[e]);
2880
2881 if (problem.has_exact_sol())
2882 {
2883 problem.exact(vals.val, tend, v_exact);
2884 problem.exact_grad(vals.val, tend, v_exact_grad);
2885 }
2886
2887 v_approx.resize(vals.val.rows(), actual_dim);
2888 v_approx.setZero();
2889
2890 v_approx_grad.resize(vals.val.rows(), mesh.dimension() * actual_dim);
2891 v_approx_grad.setZero();
2892
2893 const int n_loc_bases = int(vals.basis_values.size());
2894
2895 for (int i = 0; i < n_loc_bases; ++i)
2896 {
2897 const auto &val = vals.basis_values[i];
2898
2899 for (size_t ii = 0; ii < val.global.size(); ++ii)
2900 {
2901 for (int d = 0; d < actual_dim; ++d)
2902 {
2903 v_approx.col(d) += val.global[ii].val * sol(val.global[ii].index * actual_dim + d) * val.val;
2904 v_approx_grad.block(0, d * val.grad_t_m.cols(), v_approx_grad.rows(), val.grad_t_m.cols()) += val.global[ii].val * sol(val.global[ii].index * actual_dim + d) * val.grad_t_m;
2905 }
2906 }
2907 }
2908
2909 const auto err = problem.has_exact_sol() ? (v_exact - v_approx).eval().rowwise().norm().eval() : (v_approx).eval().rowwise().norm().eval();
2910 const auto err_grad = problem.has_exact_sol() ? (v_exact_grad - v_approx_grad).eval().rowwise().norm().eval() : (v_approx_grad).eval().rowwise().norm().eval();
2911
2912 // for(long i = 0; i < err.size(); ++i)
2913 // errors.push_back(err(i));
2914
2915 linf_err = std::max(linf_err, err.maxCoeff());
2916 grad_max_err = std::max(linf_err, err_grad.maxCoeff());
2917
2918 // {
2919 // const auto &mesh3d = *dynamic_cast<Mesh3D *>(mesh.get());
2920 // const auto v0 = mesh3d.point(mesh3d.cell_vertex(e, 0));
2921 // const auto v1 = mesh3d.point(mesh3d.cell_vertex(e, 1));
2922 // const auto v2 = mesh3d.point(mesh3d.cell_vertex(e, 2));
2923 // const auto v3 = mesh3d.point(mesh3d.cell_vertex(e, 3));
2924
2925 // Eigen::Matrix<double, 6, 3> ee;
2926 // ee.row(0) = v0 - v1;
2927 // ee.row(1) = v1 - v2;
2928 // ee.row(2) = v2 - v0;
2929
2930 // ee.row(3) = v0 - v3;
2931 // ee.row(4) = v1 - v3;
2932 // ee.row(5) = v2 - v3;
2933
2934 // Eigen::Matrix<double, 6, 1> en = ee.rowwise().norm();
2935
2936 // // Eigen::Matrix<double, 3*4, 1> alpha;
2937 // // alpha(0) = angle3(e.row(0), -e.row(1)); alpha(1) = angle3(e.row(1), -e.row(2)); alpha(2) = angle3(e.row(2), -e.row(0));
2938 // // alpha(3) = angle3(e.row(0), -e.row(4)); alpha(4) = angle3(e.row(4), e.row(3)); alpha(5) = angle3(-e.row(3), -e.row(0));
2939 // // alpha(6) = angle3(-e.row(4), -e.row(1)); alpha(7) = angle3(e.row(1), -e.row(5)); alpha(8) = angle3(e.row(5), e.row(4));
2940 // // alpha(9) = angle3(-e.row(2), -e.row(5)); alpha(10) = angle3(e.row(5), e.row(3)); alpha(11) = angle3(-e.row(3), e.row(2));
2941
2942 // const double S = (ee.row(0).cross(ee.row(1)).norm() + ee.row(0).cross(ee.row(4)).norm() + ee.row(4).cross(ee.row(1)).norm() + ee.row(2).cross(ee.row(5)).norm()) / 2;
2943 // const double V = std::abs(ee.row(3).dot(ee.row(2).cross(-ee.row(0))))/6;
2944 // const double rho = 3 * V / S;
2945 // const double hp = en.maxCoeff();
2946 // const int pp = disc_orders(e);
2947 // const int p_ref = args["space"]["discr_order"];
2948
2949 // err_per_el(e, 0) = err.mean();
2950 // err_per_el(e, 1) = err.maxCoeff();
2951 // err_per_el(e, 2) = std::pow(hp, pp+1)/(rho/hp); // /std::pow(average_edge_length, p_ref+1) * (sqrt(6)/12);
2952 // err_per_el(e, 3) = rho/hp;
2953 // err_per_el(e, 4) = (vals.det.array() * vals.quadrature.weights.array()).sum();
2954
2955 // // pred_norm += (pow(std::pow(hp, pp+1)/(rho/hp),p) * vals.det.array() * vals.quadrature.weights.array()).sum();
2956 // }
2957
2958 l2_err += (err.array() * err.array() * vals.det.array() * vals.quadrature.weights.array()).sum();
2959 h1_err += (err_grad.array() * err_grad.array() * vals.det.array() * vals.quadrature.weights.array()).sum();
2960 lp_err += (err.array().pow(p) * vals.det.array() * vals.quadrature.weights.array()).sum();
2961 }
2962
2963 h1_semi_err = sqrt(fabs(h1_err));
2964 h1_err = sqrt(fabs(l2_err) + fabs(h1_err));
2965 l2_err = sqrt(fabs(l2_err));
2966
2967 lp_err = pow(fabs(lp_err), 1. / p);
2968
2969 // pred_norm = pow(fabs(pred_norm), 1./p);
2970
2971 timer.stop();
2972 const double computing_errors_time = timer.getElapsedTime();
2973 logger().info(" took {}s", computing_errors_time);
2974
2975 logger().info("-- L2 error: {}", l2_err);
2976 logger().info("-- Lp error: {}", lp_err);
2977 logger().info("-- H1 error: {}", h1_err);
2978 logger().info("-- H1 semi error: {}", h1_semi_err);
2979 // logger().info("-- Perd norm: {}", pred_norm);
2980
2981 logger().info("-- Linf error: {}", linf_err);
2982 logger().info("-- grad max error: {}", grad_max_err);
2983
2984 // {
2985 // std::ofstream out("errs.txt");
2986 // out<<err_per_el;
2987 // out.close();
2988 // }
2989 }
2990
2991 void OutStatsData::compute_mesh_stats(const polyfem::mesh::Mesh &mesh)
2992 {
2993 using namespace polyfem::mesh;
2994
2995 simplex_count = 0;
2996 prism_count = 0;
2997 pyramid_count = 0;
2998 regular_count = 0;
2999 regular_boundary_count = 0;
3000 simple_singular_count = 0;
3001 multi_singular_count = 0;
3002 boundary_count = 0;
3003 non_regular_boundary_count = 0;
3004 non_regular_count = 0;
3005 undefined_count = 0;
3006 multi_singular_boundary_count = 0;
3007
3008 const auto &els_tag = mesh.elements_tag();
3009
3010 for (size_t i = 0; i < els_tag.size(); ++i)
3011 {
3012 const ElementType type = els_tag[i];
3013
3014 switch (type)
3015 {
3016 case ElementType::SIMPLEX:
3017 simplex_count++;
3018 break;
3019 case ElementType::PRISM:
3020 prism_count++;
3021 break;
3022 case ElementType::PYRAMID:
3023 pyramid_count++;
3024 break;
3025 case ElementType::REGULAR_INTERIOR_CUBE:
3026 regular_count++;
3027 break;
3028 case ElementType::REGULAR_BOUNDARY_CUBE:
3029 regular_boundary_count++;
3030 break;
3031 case ElementType::SIMPLE_SINGULAR_INTERIOR_CUBE:
3032 simple_singular_count++;
3033 break;
3034 case ElementType::MULTI_SINGULAR_INTERIOR_CUBE:
3035 multi_singular_count++;
3036 break;
3037 case ElementType::SIMPLE_SINGULAR_BOUNDARY_CUBE:
3038 boundary_count++;
3039 break;
3040 case ElementType::INTERFACE_CUBE:
3041 case ElementType::MULTI_SINGULAR_BOUNDARY_CUBE:
3042 multi_singular_boundary_count++;
3043 break;
3044 case ElementType::BOUNDARY_POLYTOPE:
3045 non_regular_boundary_count++;
3046 break;
3047 case ElementType::INTERIOR_POLYTOPE:
3048 non_regular_count++;
3049 break;
3050 case ElementType::UNDEFINED:
3051 undefined_count++;
3052 break;
3053 default:
3054 throw std::runtime_error("Unknown element type");
3055 }
3056 }
3057
3058 logger().info("simplex_count: \t{}", simplex_count);
3059 logger().info("prism_count: \t{}", prism_count);
3060 logger().info("pyramid_count: \t{}", pyramid_count);
3061 logger().info("regular_count: \t{}", regular_count);
3062 logger().info("regular_boundary_count: \t{}", regular_boundary_count);
3063 logger().info("simple_singular_count: \t{}", simple_singular_count);
3064 logger().info("multi_singular_count: \t{}", multi_singular_count);
3065 logger().info("boundary_count: \t{}", boundary_count);
3066 logger().info("multi_singular_boundary_count: \t{}", multi_singular_boundary_count);
3067 logger().info("non_regular_count: \t{}", non_regular_count);
3068 logger().info("non_regular_boundary_count: \t{}", non_regular_boundary_count);
3069 logger().info("undefined_count: \t{}", undefined_count);
3070 logger().info("total count:\t {}", mesh.n_elements());
3071 }
3072
3073 void OutStatsData::save_json(
3074 const nlohmann::json &args,
3075 const int n_bases, const int n_pressure_bases,
3076 const Eigen::MatrixXd &sol,
3077 const mesh::Mesh &mesh,
3078 const Eigen::VectorXi &disc_orders,
3079 const Eigen::VectorXi &disc_ordersq,
3080 const assembler::Problem &problem,
3081 const OutRuntimeData &runtime,
3082 const std::string &formulation,
3083 const bool isoparametric,
3084 const int sol_at_node_id,
3085 nlohmann::json &j) const
3086 {
3087
3088 j["args"] = args;
3089
3090 j["geom_order"] = mesh.orders().size() > 0 ? mesh.orders().maxCoeff() : 1;
3091 j["geom_order_min"] = mesh.orders().size() > 0 ? mesh.orders().minCoeff() : 1;
3092 j["discr_order_min"] = disc_orders.minCoeff();
3093 j["discr_order_max"] = disc_orders.maxCoeff();
3094 j["discr_orderq_min"] = disc_ordersq.minCoeff();
3095 j["discr_orderq_max"] = disc_ordersq.maxCoeff();
3096 j["iso_parametric"] = isoparametric;
3097 j["problem"] = problem.name();
3098 j["mat_size"] = mat_size;
3099 j["num_bases"] = n_bases;
3100 j["num_pressure_bases"] = n_pressure_bases;
3101 j["num_non_zero"] = nn_zero;
3102 j["num_flipped"] = n_flipped;
3103 j["num_dofs"] = num_dofs;
3104 j["num_vertices"] = mesh.n_vertices();
3105 j["num_elements"] = mesh.n_elements();
3106
3107 j["num_p1"] = (disc_orders.array() == 1).count();
3108 j["num_p2"] = (disc_orders.array() == 2).count();
3109 j["num_p3"] = (disc_orders.array() == 3).count();
3110 j["num_p4"] = (disc_orders.array() == 4).count();
3111 j["num_p5"] = (disc_orders.array() == 5).count();
3112
3113 j["mesh_size"] = mesh_size;
3114 j["max_angle"] = max_angle;
3115
3116 j["sigma_max"] = sigma_max;
3117 j["sigma_min"] = sigma_min;
3118 j["sigma_avg"] = sigma_avg;
3119
3120 j["min_edge_length"] = min_edge_length;
3121 j["average_edge_length"] = average_edge_length;
3122
3123 j["err_l2"] = l2_err;
3124 j["err_h1"] = h1_err;
3125 j["err_h1_semi"] = h1_semi_err;
3126 j["err_linf"] = linf_err;
3127 j["err_linf_grad"] = grad_max_err;
3128 j["err_lp"] = lp_err;
3129
3130 j["spectrum"] = {spectrum(0), spectrum(1), spectrum(2), spectrum(3)};
3131 j["spectrum_condest"] = std::abs(spectrum(3)) / std::abs(spectrum(0));
3132
3133 // j["errors"] = errors;
3134
3135 j["time_building_basis"] = runtime.building_basis_time;
3136 j["time_loading_mesh"] = runtime.loading_mesh_time;
3137 j["time_computing_poly_basis"] = runtime.computing_poly_basis_time;
3138 j["time_assembling_stiffness_mat"] = runtime.assembling_stiffness_mat_time;
3139 j["time_assembling_mass_mat"] = runtime.assembling_mass_mat_time;
3140 j["time_assigning_rhs"] = runtime.assigning_rhs_time;
3141 j["time_solving"] = runtime.solving_time;
3142 // j["time_computing_errors"] = runtime.computing_errors_time;
3143
3144 j["solver_info"] = solver_info;
3145
3146 j["count_simplex"] = simplex_count;
3147 j["count_prism"] = prism_count;
3148 j["count_pyramid"] = pyramid_count;
3149 j["count_regular"] = regular_count;
3150 j["count_regular_boundary"] = regular_boundary_count;
3151 j["count_simple_singular"] = simple_singular_count;
3152 j["count_multi_singular"] = multi_singular_count;
3153 j["count_boundary"] = boundary_count;
3154 j["count_non_regular_boundary"] = non_regular_boundary_count;
3155 j["count_non_regular"] = non_regular_count;
3156 j["count_undefined"] = undefined_count;
3157 j["count_multi_singular_boundary"] = multi_singular_boundary_count;
3158
3159 j["is_simplicial"] = mesh.n_elements() == simplex_count;
3160
3161 j["peak_memory"] = getPeakRSS() / (1024 * 1024);
3162
3163 const int actual_dim = problem.is_scalar() ? 1 : mesh.dimension();
3164
3165 std::vector<double> mmin(actual_dim);
3166 std::vector<double> mmax(actual_dim);
3167
3168 for (int d = 0; d < actual_dim; ++d)
3169 {
3170 mmin[d] = std::numeric_limits<double>::max();
3171 mmax[d] = -std::numeric_limits<double>::max();
3172 }
3173
3174 for (int i = 0; i < sol.size(); i += actual_dim)
3175 {
3176 for (int d = 0; d < actual_dim; ++d)
3177 {
3178 mmin[d] = std::min(mmin[d], sol(i + d));
3179 mmax[d] = std::max(mmax[d], sol(i + d));
3180 }
3181 }
3182
3183 std::vector<double> sol_at_node(actual_dim);
3184
3185 if (sol_at_node_id >= 0)
3186 {
3187 const int node_id = sol_at_node_id;
3188
3189 for (int d = 0; d < actual_dim; ++d)
3190 {
3191 sol_at_node[d] = sol(node_id * actual_dim + d);
3192 }
3193 }
3194
3195 j["sol_at_node"] = sol_at_node;
3196 j["sol_min"] = mmin;
3197 j["sol_max"] = mmax;
3198
3199#if defined(POLYFEM_WITH_CPP_THREADS)
3200 j["num_threads"] = utils::get_n_threads();
3201#elif defined(POLYFEM_WITH_TBB)
3202 j["num_threads"] = utils::get_n_threads();
3203#else
3204 j["num_threads"] = 1;
3205#endif
3206
3207 j["formulation"] = formulation;
3208
3209 logger().info("done");
3210 }
3211
3212} // namespace polyfem::io
val
double val
Definition Assembler.cpp:89
vals
ElementAssemblyValues vals
Definition Assembler.cpp:25
weights
std::vector< std::pair< int, double > > weights
Definition PeriodicBoundaryLagrangianForm.cpp:19
y
int y
Definition SplineBasis3d.cpp:55
z
int z
Definition SplineBasis3d.cpp:55
x
int x
Definition SplineBasis3d.cpp:55
auto_p_bases.hpp
auto_pyramid_bases.hpp
auto_q_bases.hpp
polyfem::assembler::ElementAssemblyValues
stores per element basis values at given quadrature points and geometric mapping
Definition ElementAssemblyValues.hpp:14
polyfem::assembler::ElementAssemblyValues::compute
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
Definition ElementAssemblyValues.cpp:164
polyfem::assembler::Problem::name
const std::string & name() const
Definition Problem.hpp:29
polyfem::assembler::Problem::exact_grad
virtual void exact_grad(const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const
Definition Problem.hpp:58
polyfem::assembler::Problem::is_scalar
virtual bool is_scalar() const =0
polyfem::assembler::Problem::has_exact_sol
virtual bool has_exact_sol() const =0
polyfem::assembler::Problem::exact
virtual void exact(const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const
Definition Problem.hpp:57
polyfem::basis::Basis
Represents one basis function and its gradient.
Definition Basis.hpp:44
polyfem::basis::Basis::global
const std::vector< Local2Global > & global() const
Definition Basis.hpp:104
polyfem::basis::ElementBases
Stores the basis functions for a given element in a mesh (facet in 2d, cell in 3d).
Definition ElementBases.hpp:17
polyfem::io::OutGeometryData::build_vis_boundary_mesh
void build_vis_boundary_mesh(const mesh::Mesh &mesh, const std::vector< basis::ElementBases > &gbases, const std::vector< mesh::LocalBoundary > &total_local_boundary, Eigen::MatrixXd &boundary_vis_vertices, Eigen::MatrixXd &boundary_vis_local_vertices, Eigen::MatrixXi &boundary_vis_elements, Eigen::MatrixXi &boundary_vis_elements_ids, Eigen::MatrixXi &boundary_vis_primitive_ids, Eigen::MatrixXd &boundary_vis_normals) const
builds the boundary mesh for visualization
Definition OutData.cpp:1408
polyfem::io::OutGeometryData::grid_points_bc
Eigen::MatrixXd grid_points_bc
grid mesh boundaries
Definition OutData.hpp:249
polyfem::io::OutGeometryData::build_high_order_vis_mesh
void build_high_order_vis_mesh(const mesh::Mesh &mesh, const Eigen::VectorXi &output_orders, const std::vector< basis::ElementBases > &bases, Eigen::MatrixXd &points, std::vector< paraviewo::CellElement > &elements, Eigen::MatrixXi &el_id, Eigen::MatrixXd &discr, Eigen::MatrixXd &local_points) const
builds high-der visualzation mesh per element all disconnected it also retuns the mapping to element ...
polyfem::io::OutGeometryData::grid_points
Eigen::MatrixXd grid_points
grid mesh points to export solution sampled on a grid
Definition OutData.hpp:245
polyfem::io::OutGeometryData::extract_boundary_mesh_sampled
static void extract_boundary_mesh_sampled(const mesh::Mesh &mesh, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< mesh::LocalBoundary > &total_local_boundary, Eigen::MatrixXd &node_positions, Eigen::MatrixXi &boundary_edges, Eigen::MatrixXi &boundary_triangles, std::vector< Eigen::Triplet< double > > &displacement_map_entries, const int sampling_order=0)
extracts a collision proxy sampling every boundary face on a uniform lattice of the globally maximal ...
Definition OutData.cpp:434
polyfem::io::OutGeometryData::save_volume
void save_volume(const std::string &path, const OutputSpace &space, const OutputFieldFunction &output_fields, const double t, const double dt, const ExportOptions &opts) const
saves the volume vtu file
Definition OutData.cpp:2222
polyfem::io::OutGeometryData::build_vis_mesh
void build_vis_mesh(const mesh::Mesh &mesh, const Eigen::VectorXi &disc_orders, const std::vector< basis::ElementBases > &gbases, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const bool boundary_only, Eigen::MatrixXd &points, Eigen::MatrixXi &tets, Eigen::MatrixXi &el_id, Eigen::MatrixXd &discr, Eigen::MatrixXd &local_points) const
builds visualzation mesh, upsampled mesh used for visualization the visualization mesh is a dense mes...
Definition OutData.cpp:1634
polyfem::io::OutGeometryData::build_grid
void build_grid(const polyfem::mesh::Mesh &mesh, const double spacing)
builds the grid to export the solution
Definition OutData.cpp:2604
polyfem::io::OutGeometryData::save_wire
void save_wire(const std::string &name, const OutputSpace &space, const OutputFieldFunction &output_fields, const double t, const ExportOptions &opts) const
saves the wireframe
Definition OutData.cpp:2486
polyfem::io::OutGeometryData::extract_boundary_mesh
static void extract_boundary_mesh(const mesh::Mesh &mesh, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< mesh::LocalBoundary > &total_local_boundary, Eigen::MatrixXd &node_positions, Eigen::MatrixXi &boundary_edges, Eigen::MatrixXi &boundary_triangles, std::vector< Eigen::Triplet< double > > &displacement_map_entries)
extracts the boundary mesh
Definition OutData.cpp:707
polyfem::io::OutGeometryData::save_pvd
void save_pvd(const std::string &name, const std::function< std::string(int)> &vtu_names, int time_steps, double t0, double dt, int skip_frame=1) const
save a PVD of a time dependent simulation
Definition OutData.cpp:2591
polyfem::io::OutGeometryData::save_contact_surface
void save_contact_surface(const std::string &export_surface, const OutputSpace &space, const OutputFieldFunction &output_fields, const double t, const double dt_in, const ExportOptions &opts) const
saves the surface vtu file for for constact quantites, eg contact or friction forces
Definition OutData.cpp:2446
polyfem::io::OutGeometryData::export_data
void export_data(const OutputSpace &space, const OutputFieldFunction &output_fields, const bool is_time_dependent, const double tend_in, const double dt, const ExportOptions &opts, const std::string &vis_mesh_path) const
exports everytihng, txt, vtu, etc
Definition OutData.cpp:2073
polyfem::io::OutGeometryData::save_points
void save_points(const std::string &path, const OutputSpace &space, const OutputFieldFunction &output_fields, const ExportOptions &opts) const
saves the nodal values
Definition OutData.cpp:2537
polyfem::io::OutGeometryData::save_vtu
void save_vtu(const std::string &path, const OutputSpace &space, const OutputFieldFunction &output_fields, const double t, const double dt, const ExportOptions &opts) const
saves the vtu file for time t
Definition OutData.cpp:2132
polyfem::io::OutGeometryData::init_sampler
void init_sampler(const polyfem::mesh::Mesh &mesh, const double vismesh_rel_area)
unitalize the ref element sampler
Definition OutData.cpp:2599
polyfem::io::OutGeometryData::save_surface
void save_surface(const std::string &export_surface, const OutputSpace &space, const OutputFieldFunction &output_fields, const double t, const double dt_in, const ExportOptions &opts) const
saves the surface vtu file for for surface quantites, eg traction forces
Definition OutData.cpp:2375
polyfem::io::OutGeometryData::grid_points_to_elements
Eigen::MatrixXi grid_points_to_elements
grid mesh mapping to fe elements
Definition OutData.hpp:247
polyfem::io::OutGeometryData::ref_element_sampler
utils::RefElementSampler ref_element_sampler
used to sample the solution
Definition OutData.hpp:242
polyfem::io::OutRuntimeData
timers from polyfem.
Definition OutStatsData.hpp:18
polyfem::io::OutRuntimeData::loading_mesh_time
double loading_mesh_time
time to load the mesh
Definition OutStatsData.hpp:23
polyfem::io::OutRuntimeData::assembling_stiffness_mat_time
double assembling_stiffness_mat_time
time to assembly
Definition OutStatsData.hpp:27
polyfem::io::OutRuntimeData::assigning_rhs_time
double assigning_rhs_time
time to computing the rhs
Definition OutStatsData.hpp:31
polyfem::io::OutRuntimeData::assembling_mass_mat_time
double assembling_mass_mat_time
time to assembly mass
Definition OutStatsData.hpp:29
polyfem::io::OutRuntimeData::building_basis_time
double building_basis_time
time to construct the basis
Definition OutStatsData.hpp:21
polyfem::io::OutRuntimeData::solving_time
double solving_time
time to solve
Definition OutStatsData.hpp:33
polyfem::io::OutRuntimeData::computing_poly_basis_time
double computing_poly_basis_time
time to build the polygonal/polyhedral bases
Definition OutStatsData.hpp:25
polyfem::io::OutStatsData
all stats from polyfem
Definition OutStatsData.hpp:51
polyfem::io::OutStatsData::count_flipped_elements
void count_flipped_elements(const polyfem::mesh::Mesh &mesh, const std::vector< polyfem::basis::ElementBases > &gbases)
counts the number of flipped elements
Definition OutData.cpp:2785
polyfem::io::OutStatsData::compute_errors
void compute_errors(const int n_bases, const std::vector< polyfem::basis::ElementBases > &bases, const std::vector< polyfem::basis::ElementBases > &gbases, const polyfem::mesh::Mesh &mesh, const assembler::Problem &problem, const double tend, const Eigen::MatrixXd &sol)
compute errors
Definition OutData.cpp:2830
polyfem::io::OutStatsData::compute_mesh_size
void compute_mesh_size(const polyfem::mesh::Mesh &mesh_in, const std::vector< polyfem::basis::ElementBases > &bases_in, const int n_samples, const bool use_curved_mesh_size)
computes the mesh size, it samples every edges n_samples times uses curved_mesh_size (false by defaul...
Definition OutData.cpp:2702
polyfem::io::OutStatsData::reset
void reset()
clears all stats
Definition OutData.cpp:2780
polyfem::io::OutStatsData::save_json
void save_json(const nlohmann::json &args, const int n_bases, const int n_pressure_bases, const Eigen::MatrixXd &sol, const mesh::Mesh &mesh, const Eigen::VectorXi &disc_orders, const Eigen::VectorXi &disc_ordersq, const assembler::Problem &problem, const OutRuntimeData &runtime, const std::string &formulation, const bool isoparametric, const int sol_at_node_id, nlohmann::json &j) const
saves the output statistic to a json object
Definition OutData.cpp:3073
polyfem::io::OutStatsData::compute_mesh_stats
void compute_mesh_stats(const polyfem::mesh::Mesh &mesh)
compute stats (counts els type, mesh lenght, etc), step 1 of solve
Definition OutData.cpp:2991
polyfem::mesh::LocalBoundary
Boundary primitive IDs for a single element.
Definition LocalBoundary.hpp:27
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::n_cell_faces
virtual int n_cell_faces(const int c_id) const =0
polyfem::mesh::Mesh3D::next_around_face
virtual Navigation3D::Index next_around_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_elements
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
Definition Mesh.hpp:174
polyfem::mesh::Mesh::n_vertices
virtual int n_vertices() const =0
number of vertices
polyfem::mesh::Mesh::is_polytope
bool is_polytope(const int el_id) const
checks if element is polygon compatible
Definition Mesh.cpp:450
polyfem::mesh::Mesh::get_edges
virtual void get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1) const =0
Get all the edges.
polyfem::mesh::Mesh::is_simplicial
bool is_simplicial() const
checks if the mesh is simplicial
Definition Mesh.hpp:645
polyfem::mesh::Mesh::is_conforming
virtual bool is_conforming() const =0
if the mesh is conforming
polyfem::mesh::Mesh::bounding_box
virtual void bounding_box(RowVectorNd &min, RowVectorNd &max) const =0
computes the bbox of the mesh
polyfem::mesh::Mesh::barycentric_coords
virtual void barycentric_coords(const RowVectorNd &p, const int el_id, Eigen::MatrixXd &coord) const =0
constructs barycentric coodiantes for a point p.
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::orders
const Eigen::MatrixXi & orders() const
order of each element
Definition Mesh.hpp:296
polyfem::mesh::Mesh::get_boundary_id
virtual int get_boundary_id(const int primitive) const
Get the boundary selection of an element (face in 3d, edge in 2d)
Definition Mesh.hpp:499
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::has_poly
bool has_poly() const
checks if the mesh has polytopes
Definition Mesh.hpp:603
polyfem::mesh::Mesh::dimension
int dimension() const
utily for dimension
Definition Mesh.hpp:164
polyfem::mesh::Mesh::n_faces
virtual int n_faces() const =0
number of faces
polyfem::mesh::Mesh::elements_tag
const std::vector< ElementType > & elements_tag() const
Returns the elements types.
Definition Mesh.hpp:439
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::n_face_vertices
virtual int n_face_vertices(const int f_id) const =0
number of vertices of a face
polyfem::mesh::Mesh::elements_boxes
virtual void elements_boxes(std::vector< std::array< Eigen::Vector3d, 2 > > &boxes) const =0
constructs a box around every element (3d cell, 2d face)
polyfem::mesh::Mesh::is_boundary_element
virtual bool is_boundary_element(const int element_global_id) const =0
is cell boundary
polyfem::mesh::Mesh::get_node_id
virtual int get_node_id(const int node_id) const
Get the boundary selection of a node.
Definition Mesh.hpp:508
polyfem::mesh::Obstacle::get_edge_connectivity
const Eigen::MatrixXi & get_edge_connectivity() const
Definition Obstacle.hpp:48
polyfem::mesh::Obstacle::get_face_connectivity
const Eigen::MatrixXi & get_face_connectivity() const
Definition Obstacle.hpp:47
polyfem::mesh::Obstacle::v
const Eigen::MatrixXd & v() const
Definition Obstacle.hpp:42
polyfem::mesh::Obstacle::get_vertex_connectivity
const Eigen::VectorXi & get_vertex_connectivity() const
Definition Obstacle.hpp:49
polyfem::utils::BoundarySampler::sample_parametric_prism_face
static void sample_parametric_prism_face(int index, int n_samples, Eigen::MatrixXd &uv, Eigen::MatrixXd &samples)
Definition BoundarySampler.cpp:512
polyfem::utils::BoundarySampler::normal_for_quad_edge
static void normal_for_quad_edge(int index, Eigen::MatrixXd &normal)
Definition BoundarySampler.cpp:709
polyfem::utils::BoundarySampler::normal_for_tri_edge
static void normal_for_tri_edge(int index, Eigen::MatrixXd &normal)
Definition BoundarySampler.cpp:719
polyfem::utils::BoundarySampler::normal_for_quad_face
static void normal_for_quad_face(int index, Eigen::MatrixXd &normal)
Definition BoundarySampler.cpp:729
polyfem::utils::BoundarySampler::sample_parametric_pyramid_face
static void sample_parametric_pyramid_face(int index, int n_samples, Eigen::MatrixXd &uv, Eigen::MatrixXd &samples)
Definition BoundarySampler.cpp:570
polyfem::utils::BoundarySampler::sample_parametric_tri_face
static void sample_parametric_tri_face(int index, int n_samples, Eigen::MatrixXd &uv, Eigen::MatrixXd &samples)
Definition BoundarySampler.cpp:481
polyfem::utils::BoundarySampler::normal_for_prism_face
static void normal_for_prism_face(int index, Eigen::MatrixXd &normal)
Definition BoundarySampler.cpp:753
polyfem::utils::BoundarySampler::normal_for_tri_face
static void normal_for_tri_face(int index, Eigen::MatrixXd &normal)
Definition BoundarySampler.cpp:741
polyfem::utils::BoundarySampler::sample_parametric_quad_face
static void sample_parametric_quad_face(int index, int n_samples, Eigen::MatrixXd &uv, Eigen::MatrixXd &samples)
Definition BoundarySampler.cpp:456
polyfem::utils::BoundarySampler::normal_for_polygon_edge
static void normal_for_polygon_edge(int face_id, int edge_id, const mesh::Mesh &mesh, Eigen::MatrixXd &normal)
Definition BoundarySampler.cpp:777
polyfem::utils::BoundarySampler::normal_for_pyramid_face
static void normal_for_pyramid_face(int index, Eigen::MatrixXd &normal)
Definition BoundarySampler.cpp:765
polyfem::utils::BoundarySampler::sample_parametric_quad_edge
static void sample_parametric_quad_edge(int index, int n_samples, Eigen::MatrixXd &uv, Eigen::MatrixXd &samples)
Definition BoundarySampler.cpp:422
polyfem::utils::BoundarySampler::sample_polygon_edge
static void sample_polygon_edge(int face_id, int edge_id, int n_samples, const mesh::Mesh &mesh, Eigen::MatrixXd &uv, Eigen::MatrixXd &samples)
Definition BoundarySampler.cpp:628
polyfem::utils::BoundarySampler::sample_parametric_tri_edge
static void sample_parametric_tri_edge(int index, int n_samples, Eigen::MatrixXd &uv, Eigen::MatrixXd &samples)
Definition BoundarySampler.cpp:439
polyfem::utils::EdgeSampler::sample_3d_simplex
static void sample_3d_simplex(const int resolution, Eigen::MatrixXd &samples)
Definition EdgeSampler.cpp:63
polyfem::utils::EdgeSampler::sample_3d_cube
static void sample_3d_cube(const int resolution, Eigen::MatrixXd &samples)
Definition EdgeSampler.cpp:103
polyfem::utils::EdgeSampler::sample_2d_cube
static void sample_2d_cube(const int resolution, Eigen::MatrixXd &samples)
Definition EdgeSampler.cpp:9
polyfem::utils::EdgeSampler::sample_2d_simplex
static void sample_2d_simplex(const int resolution, Eigen::MatrixXd &samples)
Definition EdgeSampler.cpp:38
polyfem::utils::RefElementSampler::init
void init(const bool is_volume, const int n_elements, const double target_rel_area)
Definition RefElementSampler.cpp:330
polyfem::utils::RefElementSampler::pyramid_volume
const Eigen::MatrixXi & pyramid_volume() const
Definition RefElementSampler.hpp:65
polyfem::utils::RefElementSampler::pyramid_points
const Eigen::MatrixXd & pyramid_points() const
Definition RefElementSampler.hpp:63
polyfem::utils::RefElementSampler::simplex_volume
const Eigen::MatrixXi & simplex_volume() const
Definition RefElementSampler.hpp:53
getPeakRSS
size_t getPeakRSS(void)
Returns the peak (maximum so far) resident set size (physical memory use) measured in bytes,...
Definition getRSS.c:37
p_bases.vertices
list vertices
Definition p_bases.py:238
p_bases.points
points
Definition p_bases.py:252
polyfem::autogen::q_nodes_2d
void q_nodes_2d(const int q, Eigen::MatrixXd &val)
Definition auto_q_bases_2d_nodes.cpp:67
polyfem::autogen::pyramid_nodes_3d
void pyramid_nodes_3d(const int pyramid, Eigen::MatrixXd &val)
Definition auto_pyramid_bases.cpp:8341
polyfem::autogen::prism_nodes_3d
void prism_nodes_3d(const int p, const int q, Eigen::MatrixXd &val)
Definition prism_bases.cpp:92
polyfem::autogen::p_nodes_2d
void p_nodes_2d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:969
polyfem::autogen::p_nodes_3d
void p_nodes_3d(const int p, Eigen::MatrixXd &val)
Definition auto_p_bases.cpp:3060
polyfem::autogen::q_nodes_3d
void q_nodes_3d(const int q, Eigen::MatrixXd &val)
Definition auto_q_bases_3d_nodes.cpp:149
polyfem::io::OutputFieldFunction
std::function< std::vector< OutputField >(const OutputSample &)> OutputFieldFunction
Definition OutputData.hpp:73
polyfem::io::CellElement
paraviewo::CellElement CellElement
Definition OutData.cpp:62
polyfem::io::CellType
paraviewo::CellType CellType
Definition OutData.cpp:61
polyfem::io::write_matrix
bool write_matrix(const std::string &path, const Mat &mat)
Writes a matrix to a file. Determines the file format based on the path's extension.
Definition MatrixIO.cpp:42
polyfem::mesh::ElementType
ElementType
Type of Element, check [Poly-Spline Finite Element Method] for a complete description.
Definition Mesh.hpp:31
polyfem::utils::get_n_threads
size_t get_n_threads()
Definition par_for.hpp:51
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
polyfem::json
nlohmann::json json
Definition Common.hpp:9
polyfem::RowVectorNd
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13
polyfem::log_and_throw_error
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
par_for.hpp
prism_bases.hpp
polyfem::io::OutGeometryData::ExportOptions::file_extension
std::string file_extension() const
return the extension of the output paraview files depending on use_hdf5
Definition OutData.hpp:65
polyfem::io::OutGeometryData::ExportOptions::fields
std::vector< std::string > fields
Definition OutData.hpp:34
polyfem::io::OutGeometryData::ExportOptions::ExportOptions
ExportOptions(const json &args, const bool is_mesh_linear, const bool mesh_has_prisms, const bool is_problem_scalar)
initialize the flags based on the input args
Definition OutData.cpp:2105
polyfem::io::OutGeometryData::ExportOptions::export_field
bool export_field(const std::string &field) const
Definition OutData.cpp:2100
polyfem::io::OutputFieldOptions::export_field
bool export_field(const std::string &field) const
Definition OutData.cpp:56
polyfem::io::OutputFieldOptions::fields
std::vector< std::string > fields
Definition OutputData.hpp:29
polyfem::io::OutputSample::node_ids
Eigen::VectorXi node_ids
Definition OutputData.hpp:51
polyfem::io::OutputSample::normals
Eigen::MatrixXd normals
Definition OutputData.hpp:52
polyfem::io::OutputSample::primitive_ids
Eigen::VectorXi primitive_ids
Definition OutputData.hpp:50
polyfem::io::OutputSample::requested_fields
std::vector< std::string > requested_fields
Definition OutputData.hpp:53
polyfem::io::OutputSample::element_ids
Eigen::VectorXi element_ids
Definition OutputData.hpp:49
polyfem::io::OutputSample::local_points
Eigen::MatrixXd local_points
Definition OutputData.hpp:48
polyfem::io::OutputSpace::mesh
const mesh::Mesh * mesh
Definition OutputData.hpp:77
polyfem::io::OutputSpace::output_orders
Eigen::VectorXi output_orders
Definition OutputData.hpp:79
polyfem::io::OutputSpace::total_local_boundary
const std::vector< mesh::LocalBoundary > * total_local_boundary
Definition OutputData.hpp:82
polyfem::io::OutputSpace::geometry_bases
const std::vector< basis::ElementBases > * geometry_bases
Definition OutputData.hpp:78
polyfem::io::OutputSpace::dirichlet_nodes_position
const std::vector< RowVectorNd > * dirichlet_nodes_position
Definition OutputData.hpp:86
polyfem::io::OutputSpace::dirichlet_nodes
const std::vector< int > * dirichlet_nodes
Definition OutputData.hpp:85
polyfem::io::OutputSpace::polys
const std::map< int, Eigen::MatrixXd > * polys
Definition OutputData.hpp:80
polyfem::io::OutputSpace::obstacle
const mesh::Obstacle * obstacle
Definition OutputData.hpp:83
polyfem::io::OutputSpace::collision_mesh
const ipc::CollisionMesh * collision_mesh
Definition OutputData.hpp:84
polyfem::io::OutputSpace::polys_3d
const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > * polys_3d
Definition OutputData.hpp:81