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