PolyFEM
Loading...
Searching...
No Matches
OutData.cpp
Go to the documentation of this file.
1#include "OutData.hpp"
2
3#include <polyfem/io/Evaluator.hpp>
4
5#include <polyfem/legacy/State.hpp>
6
7#include <polyfem/assembler/ElementAssemblyValues.hpp>
8#include <polyfem/assembler/AssemblyValues.hpp>
9#include <polyfem/assembler/MatParams.hpp>
10#include <polyfem/assembler/Mass.hpp>
11
12#include <polyfem/basis/ElementBases.hpp>
13
14#include <polyfem/mesh/MeshUtils.hpp>
15#include <polyfem/mesh/mesh2D/Mesh2D.hpp>
16#include <polyfem/mesh/mesh3D/Mesh3D.hpp>
17
18#include <polyfem/time_integrator/ImplicitTimeIntegrator.hpp>
19
20#include <polyfem/solver/forms/SmoothContactForm.hpp>
21#include <polyfem/solver/forms/FrictionForm.hpp>
22#include <polyfem/solver/NLProblem.hpp>
23#include <polyfem/solver/forms/BodyForm.hpp>
24#include <polyfem/solver/forms/BodyForm.hpp>
25#include <polyfem/solver/forms/BarrierContactForm.hpp>
26#include <polyfem/solver/forms/ElasticForm.hpp>
27#include <polyfem/solver/forms/FrictionForm.hpp>
28#include <polyfem/solver/forms/InertiaForm.hpp>
29#include <polyfem/solver/forms/LaggedRegForm.hpp>
30#include <polyfem/solver/forms/RayleighDampingForm.hpp>
31#include <polyfem/solver/forms/lagrangian/AugmentedLagrangianForm.hpp>
32
33#include <polyfem/utils/EdgeSampler.hpp>
34#include <polyfem/utils/Logger.hpp>
35#include <polyfem/utils/par_for.hpp>
36#include <polyfem/utils/BoundarySampler.hpp>
37#include <polyfem/utils/Timer.hpp>
38#include <polyfem/utils/MaybeParallelFor.hpp>
39#include <polyfem/utils/getRSS.h>
40
41#include <polyfem/autogen/auto_p_bases.hpp>
42#include <polyfem/autogen/auto_q_bases.hpp>
43#include <polyfem/autogen/prism_bases.hpp>
44#include <polyfem/autogen/auto_pyramid_bases.hpp>
45
46#include <paraviewo/VTMWriter.hpp>
47#include <paraviewo/PVDWriter.hpp>
48
49#include <SimpleBVH/BVH.hpp>
50
51#include <igl/write_triangle_mesh.h>
52#include <igl/edges.h>
53#include <igl/facet_adjacency_matrix.h>
54#include <igl/connected_components.h>
55
56#include <ipc/ipc.hpp>
57
58#include <filesystem>
59
60namespace polyfem::legacy::io
61{
62 using polyfem::io::Evaluator;
63 using CellType = paraviewo::CellType;
64 using CellElement = paraviewo::CellElement;
65
66 namespace
67 {
68 void compute_traction_forces(const State &state, const Eigen::MatrixXd &solution, const double t, Eigen::MatrixXd &traction_forces, bool skip_dirichlet = true)
69 {
70 int actual_dim = 1;
71 if (!state.problem->is_scalar())
72 actual_dim = state.mesh->dimension();
73 else
74 return;
75
76 const std::vector<basis::ElementBases> &bases = state.bases;
77 const std::vector<basis::ElementBases> &gbases = state.geom_bases();
78
79 Eigen::MatrixXd uv, samples, gtmp, rhs_fun, deform_mat, trafo;
80 Eigen::VectorXi global_primitive_ids;
81 Eigen::MatrixXd points, normals;
82 Eigen::VectorXd weights;
83 polyfem::assembler::ElementAssemblyValues vals;
84
85 traction_forces.setZero(state.n_bases * actual_dim, 1);
86
87 for (const auto &lb : state.total_local_boundary)
88 {
89 const int e = lb.element_id();
90 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, false, uv, points, normals, weights, global_primitive_ids);
91
92 if (!has_samples)
93 continue;
94
95 const basis::ElementBases &gbs = gbases[e];
96 const basis::ElementBases &bs = bases[e];
97
98 vals.compute(e, state.mesh->is_volume(), points, bs, gbs);
99
100 for (int n = 0; n < normals.rows(); ++n)
101 {
102 trafo = vals.jac_it[n].inverse();
103
104 if (solution.size() > 0)
105 {
106 assert(actual_dim == 2 || actual_dim == 3);
107 deform_mat.resize(actual_dim, actual_dim);
108 deform_mat.setZero();
109 for (const auto &b : vals.basis_values)
110 {
111 for (const auto &g : b.global)
112 {
113 for (int d = 0; d < actual_dim; ++d)
114 {
115 deform_mat.row(d) += solution(g.index * actual_dim + d) * b.grad.row(n);
116 }
117 }
118 }
119
120 trafo += deform_mat;
121 }
122
123 normals.row(n) = normals.row(n) * trafo.inverse();
124 normals.row(n).normalize();
125 }
126
127 std::vector<assembler::Assembler::NamedMatrix> tensor_flat;
128 state.assembler->compute_tensor_value(assembler::OutputData(t, e, bs, gbs, points, solution), tensor_flat);
129
130 for (long n = 0; n < vals.basis_values.size(); ++n)
131 {
132 const polyfem::assembler::AssemblyValues &v = vals.basis_values[n];
133
134 const int g_index = v.global[0].index * actual_dim;
135
136 for (int q = 0; q < points.rows(); ++q)
137 {
138 // TF computed only from cauchy stress
139 assert(tensor_flat[0].first == "cauchy_stess");
140 assert(tensor_flat[0].second.row(q).size() == actual_dim * actual_dim);
141
142 Eigen::MatrixXd stress_tensor = utils::unflatten(tensor_flat[0].second.row(q), actual_dim);
143
144 traction_forces.block(g_index, 0, actual_dim, 1) += stress_tensor * normals.row(q).transpose() * v.val(q) * weights(q);
145 }
146 }
147 }
148 }
149 } // namespace
150
151 void OutGeometryData::extract_boundary_mesh(
152 const mesh::Mesh &mesh,
153 const int n_bases,
154 const std::vector<basis::ElementBases> &bases,
155 const std::vector<mesh::LocalBoundary> &total_local_boundary,
156 Eigen::MatrixXd &node_positions,
157 Eigen::MatrixXi &boundary_edges,
158 Eigen::MatrixXi &boundary_triangles,
159 std::vector<Eigen::Triplet<double>> &displacement_map_entries)
160 {
161 using namespace polyfem::mesh;
162
163 displacement_map_entries.clear();
164
165 if (mesh.is_volume())
166 {
167 if (mesh.has_poly())
168 {
169 logger().warn("Skipping as the mesh has polygons");
170 return;
171 }
172
173 const bool is_simplicial = mesh.is_simplicial();
174
175 node_positions.resize(n_bases + (is_simplicial ? 0 : mesh.n_faces()), 3);
176 node_positions.setZero();
177 const Mesh3D &mesh3d = dynamic_cast<const Mesh3D &>(mesh);
178
179 std::vector<std::tuple<int, int, int>> tris;
180
181 std::vector<bool> visited_node(n_bases, false);
182
183 std::stringstream print_warning;
184
185 for (const LocalBoundary &lb : total_local_boundary)
186 {
187 const basis::ElementBases &b = bases[lb.element_id()];
188
189 for (int j = 0; j < lb.size(); ++j)
190 {
191 const int eid = lb.global_primitive_id(j);
192 const int lid = lb[j];
193 const Eigen::VectorXi nodes = b.local_nodes_for_primitive(eid, mesh3d);
194
195 if (mesh.is_cube(lb.element_id()))
196 {
197 assert(!is_simplicial);
198 assert(!mesh.has_poly());
199 std::vector<int> loc_nodes;
200 RowVectorNd bary = RowVectorNd::Zero(3);
201
202 for (long n = 0; n < nodes.size(); ++n)
203 {
204 auto &bs = b.bases[nodes(n)];
205 const auto &glob = bs.global();
206 if (glob.size() != 1)
207 continue;
208
209 int gindex = glob.front().index;
210 node_positions.row(gindex) = glob.front().node;
211 bary += glob.front().node;
212 loc_nodes.push_back(gindex);
213 }
214
215 if (loc_nodes.size() != 4)
216 {
217 logger().trace("skipping element {} since it is not Q1", eid);
218 continue;
219 }
220
221 bary /= 4;
222
223 const int new_node = n_bases + eid;
224 node_positions.row(new_node) = bary;
225 tris.emplace_back(loc_nodes[1], loc_nodes[0], new_node);
226 tris.emplace_back(loc_nodes[2], loc_nodes[1], new_node);
227 tris.emplace_back(loc_nodes[3], loc_nodes[2], new_node);
228 tris.emplace_back(loc_nodes[0], loc_nodes[3], new_node);
229
230 for (int q = 0; q < 4; ++q)
231 {
232 if (!visited_node[loc_nodes[q]])
233 displacement_map_entries.emplace_back(loc_nodes[q], loc_nodes[q], 1);
234
235 visited_node[loc_nodes[q]] = true;
236 displacement_map_entries.emplace_back(new_node, loc_nodes[q], 0.25);
237 }
238
239 continue;
240 }
241 else if (mesh.is_prism(lb.element_id()))
242 {
243 assert(!is_simplicial);
244 assert(!mesh.has_poly());
245 std::vector<int> loc_nodes;
246 RowVectorNd bary = RowVectorNd::Zero(3);
247
248 for (long n = 0; n < nodes.size(); ++n)
249 {
250 auto &bs = b.bases[nodes(n)];
251 const auto &glob = bs.global();
252 if (glob.size() != 1)
253 continue;
254
255 int gindex = glob.front().index;
256 node_positions.row(gindex) = glob.front().node;
257 bary += glob.front().node;
258 loc_nodes.push_back(gindex);
259 }
260
261 auto update_mapping = [&displacement_map_entries, &visited_node](const std::vector<int> &loc_nodes) {
262 for (int k = 0; k < loc_nodes.size(); ++k)
263 {
264 if (!visited_node[loc_nodes[k]])
265 displacement_map_entries.emplace_back(loc_nodes[k], loc_nodes[k], 1);
266
267 visited_node[loc_nodes[k]] = true;
268 }
269 };
270
271 if (loc_nodes.size() == 3)
272 {
273 tris.emplace_back(loc_nodes[0], loc_nodes[1], loc_nodes[2]);
274
275 update_mapping(loc_nodes);
276 }
277 else if (loc_nodes.size() == 6)
278 {
279 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[5]);
280 tris.emplace_back(loc_nodes[3], loc_nodes[1], loc_nodes[4]);
281 tris.emplace_back(loc_nodes[4], loc_nodes[2], loc_nodes[5]);
282 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[5]);
283
284 update_mapping(loc_nodes);
285 }
286 else if (loc_nodes.size() == 10)
287 {
288 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[8]);
289 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[9]);
290 tris.emplace_back(loc_nodes[4], loc_nodes[1], loc_nodes[5]);
291 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[9]);
292 tris.emplace_back(loc_nodes[6], loc_nodes[2], loc_nodes[7]);
293 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[9]);
294 tris.emplace_back(loc_nodes[8], loc_nodes[3], loc_nodes[9]);
295 tris.emplace_back(loc_nodes[9], loc_nodes[4], loc_nodes[5]);
296 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[9]);
297 update_mapping(loc_nodes);
298 }
299 else if (loc_nodes.size() == 15)
300 {
301 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[11]);
302 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[12]);
303 tris.emplace_back(loc_nodes[3], loc_nodes[12], loc_nodes[11]);
304 tris.emplace_back(loc_nodes[12], loc_nodes[10], loc_nodes[11]);
305 tris.emplace_back(loc_nodes[4], loc_nodes[5], loc_nodes[13]);
306 tris.emplace_back(loc_nodes[4], loc_nodes[13], loc_nodes[12]);
307 tris.emplace_back(loc_nodes[12], loc_nodes[13], loc_nodes[14]);
308 tris.emplace_back(loc_nodes[12], loc_nodes[14], loc_nodes[10]);
309 tris.emplace_back(loc_nodes[14], loc_nodes[9], loc_nodes[10]);
310 tris.emplace_back(loc_nodes[5], loc_nodes[1], loc_nodes[6]);
311 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[13]);
312 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[13]);
313 tris.emplace_back(loc_nodes[13], loc_nodes[7], loc_nodes[14]);
314 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[14]);
315 tris.emplace_back(loc_nodes[14], loc_nodes[8], loc_nodes[9]);
316 tris.emplace_back(loc_nodes[8], loc_nodes[2], loc_nodes[9]);
317 update_mapping(loc_nodes);
318 }
319 else if (loc_nodes.size() == 4)
320 {
321 bary /= 4;
322
323 const int new_node = n_bases + eid;
324 node_positions.row(new_node) = bary;
325 tris.emplace_back(loc_nodes[1], loc_nodes[0], new_node);
326 tris.emplace_back(loc_nodes[2], loc_nodes[1], new_node);
327 tris.emplace_back(loc_nodes[3], loc_nodes[2], new_node);
328 tris.emplace_back(loc_nodes[0], loc_nodes[3], new_node);
329
330 update_mapping(loc_nodes);
331 }
332 else
333 {
334 logger().trace("skipping element {} since it is not linear, it has {} nodes", eid, loc_nodes.size());
335 continue;
336 }
337
338 continue;
339 }
340
341 else if (mesh.is_pyramid(lb.element_id()))
342 {
343 assert(!is_simplicial);
344 assert(!mesh.has_poly());
345 std::vector<int> loc_nodes;
346
347 for (long n = 0; n < nodes.size(); ++n)
348 {
349 auto &bs = b.bases[nodes(n)];
350 const auto &glob = bs.global();
351 if (glob.size() != 1)
352 continue;
353
354 int gindex = glob.front().index;
355 node_positions.row(gindex) = glob.front().node;
356 loc_nodes.push_back(gindex);
357 }
358
359 auto update_mapping = [&displacement_map_entries, &visited_node](const std::vector<int> &loc_nodes) {
360 for (int k = 0; k < loc_nodes.size(); ++k)
361 {
362 if (!visited_node[loc_nodes[k]])
363 displacement_map_entries.emplace_back(loc_nodes[k], loc_nodes[k], 1);
364
365 visited_node[loc_nodes[k]] = true;
366 }
367 };
368
369 if (loc_nodes.size() == 3)
370 {
371 tris.emplace_back(loc_nodes[0], loc_nodes[1], loc_nodes[2]);
372 update_mapping(loc_nodes);
373 }
374 else if (loc_nodes.size() == 4)
375 {
376 tris.emplace_back(loc_nodes[0], loc_nodes[1], loc_nodes[2]);
377 tris.emplace_back(loc_nodes[0], loc_nodes[2], loc_nodes[3]);
378 update_mapping(loc_nodes);
379 }
380 else
381 {
382 logger().trace("skipping element {} since it is not linear, it has {} nodes", eid, loc_nodes.size());
383 continue;
384 }
385
386 continue;
387 }
388
389 if (!mesh.is_simplex(lb.element_id()))
390 {
391 logger().trace("skipping element {} since it is not a simplex or hex", eid);
392 continue;
393 }
394
395 assert(mesh.is_simplex(lb.element_id()));
396
397 std::vector<int> loc_nodes;
398
399 bool is_follower = false;
400 if (!mesh3d.is_conforming())
401 {
402 for (long n = 0; n < nodes.size(); ++n)
403 {
404 auto &bs = b.bases[nodes(n)];
405 const auto &glob = bs.global();
406 if (glob.size() != 1)
407 {
408 is_follower = true;
409 break;
410 }
411 }
412 }
413
414 if (is_follower)
415 continue;
416
417 for (long n = 0; n < nodes.size(); ++n)
418 {
419 const basis::Basis &bs = b.bases[nodes(n)];
420 const std::vector<basis::Local2Global> &glob = bs.global();
421 if (glob.size() != 1)
422 continue;
423
424 int gindex = glob.front().index;
425 node_positions.row(gindex) = glob.front().node;
426 loc_nodes.push_back(gindex);
427 }
428
429 if (loc_nodes.size() == 3)
430 {
431 tris.emplace_back(loc_nodes[0], loc_nodes[1], loc_nodes[2]);
432 }
433 else if (loc_nodes.size() == 6)
434 {
435 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[5]);
436 tris.emplace_back(loc_nodes[3], loc_nodes[1], loc_nodes[4]);
437 tris.emplace_back(loc_nodes[4], loc_nodes[2], loc_nodes[5]);
438 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[5]);
439 }
440 else if (loc_nodes.size() == 10)
441 {
442 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[8]);
443 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[9]);
444 tris.emplace_back(loc_nodes[4], loc_nodes[1], loc_nodes[5]);
445 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[9]);
446 tris.emplace_back(loc_nodes[6], loc_nodes[2], loc_nodes[7]);
447 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[9]);
448 tris.emplace_back(loc_nodes[8], loc_nodes[3], loc_nodes[9]);
449 tris.emplace_back(loc_nodes[9], loc_nodes[4], loc_nodes[5]);
450 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[9]);
451 }
452 else if (loc_nodes.size() == 15)
453 {
454 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[11]);
455 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[12]);
456 tris.emplace_back(loc_nodes[3], loc_nodes[12], loc_nodes[11]);
457 tris.emplace_back(loc_nodes[12], loc_nodes[10], loc_nodes[11]);
458 tris.emplace_back(loc_nodes[4], loc_nodes[5], loc_nodes[13]);
459 tris.emplace_back(loc_nodes[4], loc_nodes[13], loc_nodes[12]);
460 tris.emplace_back(loc_nodes[12], loc_nodes[13], loc_nodes[14]);
461 tris.emplace_back(loc_nodes[12], loc_nodes[14], loc_nodes[10]);
462 tris.emplace_back(loc_nodes[14], loc_nodes[9], loc_nodes[10]);
463 tris.emplace_back(loc_nodes[5], loc_nodes[1], loc_nodes[6]);
464 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[13]);
465 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[13]);
466 tris.emplace_back(loc_nodes[13], loc_nodes[7], loc_nodes[14]);
467 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[14]);
468 tris.emplace_back(loc_nodes[14], loc_nodes[8], loc_nodes[9]);
469 tris.emplace_back(loc_nodes[8], loc_nodes[2], loc_nodes[9]);
470 }
471 else
472 {
473 print_warning << loc_nodes.size() << " ";
474 // assert(false);
475 }
476
477 if (!is_simplicial)
478 {
479 for (int k = 0; k < loc_nodes.size(); ++k)
480 {
481 if (!visited_node[loc_nodes[k]])
482 displacement_map_entries.emplace_back(loc_nodes[k], loc_nodes[k], 1);
483
484 visited_node[loc_nodes[k]] = true;
485 }
486 }
487 }
488 }
489
490 if (print_warning.str().size() > 0)
491 logger().warn("Skipping faces as theys have {} nodes, boundary export supported up to p4", print_warning.str());
492
493 boundary_triangles.resize(tris.size(), 3);
494 for (int i = 0; i < tris.size(); ++i)
495 {
496 boundary_triangles.row(i) << std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]);
497 }
498
499 if (boundary_triangles.rows() > 0)
500 {
501 igl::edges(boundary_triangles, boundary_edges);
502 }
503 }
504 else
505 {
506 node_positions.resize(n_bases, 2);
507 node_positions.setZero();
508 const Mesh2D &mesh2d = dynamic_cast<const Mesh2D &>(mesh);
509
510 std::vector<std::pair<int, int>> edges;
511
512 for (const LocalBoundary &lb : total_local_boundary)
513 {
514 const basis::ElementBases &b = bases[lb.element_id()];
515
516 for (int j = 0; j < lb.size(); ++j)
517 {
518 const int eid = lb.global_primitive_id(j);
519 const int lid = lb[j];
520 const Eigen::VectorXi nodes = b.local_nodes_for_primitive(eid, mesh2d);
521
522 int prev_node = -1;
523
524 for (long n = 0; n < nodes.size(); ++n)
525 {
526 const basis::Basis &bs = b.bases[nodes(n)];
527 const std::vector<basis::Local2Global> &glob = bs.global();
528 if (glob.size() != 1)
529 continue;
530
531 int gindex = glob.front().index;
532 node_positions.row(gindex) = glob.front().node.head<2>();
533
534 if (prev_node >= 0)
535 edges.emplace_back(prev_node, gindex);
536
537 prev_node = gindex;
538 }
539 }
540 }
541
542 boundary_triangles.resize(0, 0);
543 boundary_edges.resize(edges.size(), 2);
544 for (int i = 0; i < edges.size(); ++i)
545 {
546 boundary_edges.row(i) << edges[i].first, edges[i].second;
547 }
548 }
549 }
550
551 void OutGeometryData::build_vis_boundary_mesh(
552 const mesh::Mesh &mesh,
553 const std::vector<basis::ElementBases> &bases,
554 const std::vector<basis::ElementBases> &gbases,
555 const std::vector<mesh::LocalBoundary> &total_local_boundary,
556 const Eigen::MatrixXd &solution,
557 const int problem_dim,
558 Eigen::MatrixXd &boundary_vis_vertices,
559 Eigen::MatrixXd &boundary_vis_local_vertices,
560 Eigen::MatrixXi &boundary_vis_elements,
561 Eigen::MatrixXi &boundary_vis_elements_ids,
562 Eigen::MatrixXi &boundary_vis_primitive_ids,
563 Eigen::MatrixXd &boundary_vis_normals,
564 Eigen::MatrixXd &displaced_boundary_vis_normals) const
565 {
566 using namespace polyfem::mesh;
567
568 std::vector<Eigen::MatrixXd> lv, vertices, allnormals, displaced_allnormals;
569 std::vector<int> el_ids, global_primitive_ids;
570 Eigen::MatrixXd uv, local_pts, tmp_n, normals, displaced_normals, trafo, deform_mat;
571 assembler::ElementAssemblyValues vals;
572 const auto &sampler = ref_element_sampler;
573 const int n_samples = sampler.num_samples();
574 int size = 0;
575
576 std::vector<std::pair<int, int>> edges;
577 std::vector<std::tuple<int, int, int>> tris;
578
579 for (auto it = total_local_boundary.begin(); it != total_local_boundary.end(); ++it)
580 {
581 const auto &lb = *it;
582 const auto &gbs = gbases[lb.element_id()];
583 const auto &bs = bases[lb.element_id()];
584
585 for (int k = 0; k < lb.size(); ++k)
586 {
587 switch (lb.type())
588 {
589 case BoundaryType::TRI_LINE:
590 utils::BoundarySampler::normal_for_tri_edge(lb[k], tmp_n);
591 utils::BoundarySampler::sample_parametric_tri_edge(lb[k], n_samples, uv, local_pts);
592 break;
593 case BoundaryType::QUAD_LINE:
594 utils::BoundarySampler::normal_for_quad_edge(lb[k], tmp_n);
595 utils::BoundarySampler::sample_parametric_quad_edge(lb[k], n_samples, uv, local_pts);
596 break;
597 case BoundaryType::QUAD:
598 utils::BoundarySampler::normal_for_quad_face(lb[k], tmp_n);
599 utils::BoundarySampler::sample_parametric_quad_face(lb[k], n_samples, uv, local_pts);
600 break;
601 case BoundaryType::TRI:
602 utils::BoundarySampler::normal_for_tri_face(lb[k], tmp_n);
603 utils::BoundarySampler::sample_parametric_tri_face(lb[k], n_samples, uv, local_pts);
604 break;
605 case BoundaryType::PRISM:
606 utils::BoundarySampler::normal_for_prism_face(lb[k], tmp_n);
607 utils::BoundarySampler::sample_parametric_prism_face(lb[k], n_samples, uv, local_pts);
608 break;
609 case BoundaryType::PYRAMID:
610 utils::BoundarySampler::normal_for_pyramid_face(lb[k], tmp_n);
611 utils::BoundarySampler::sample_parametric_pyramid_face(lb[k], n_samples, uv, local_pts);
612 break;
613 case BoundaryType::POLYGON:
614 utils::BoundarySampler::normal_for_polygon_edge(lb.element_id(), lb.global_primitive_id(k), mesh, tmp_n);
615 utils::BoundarySampler::sample_polygon_edge(lb.element_id(), lb.global_primitive_id(k), n_samples, mesh, uv, local_pts);
616 break;
617 case BoundaryType::POLYHEDRON:
618 assert(false);
619 break;
620 case BoundaryType::INVALID:
621 assert(false);
622 break;
623 default:
624 assert(false);
625 }
626
627 vertices.emplace_back();
628 lv.emplace_back(local_pts);
629 el_ids.push_back(lb.element_id());
630 global_primitive_ids.push_back(lb.global_primitive_id(k));
631 gbs.eval_geom_mapping(local_pts, vertices.back());
632 vals.compute(lb.element_id(), mesh.is_volume(), local_pts, bs, gbs);
633 const int tris_start = tris.size();
634
635 if (mesh.is_volume())
636 {
637 const bool prism_quad = lb.type() == BoundaryType::PRISM && lb[k] >= 2;
638 const bool prism_tri = lb.type() == BoundaryType::PRISM && lb[k] < 2;
639
640 const bool pyramid_quad = lb.type() == BoundaryType::PYRAMID && lb[k] == 0;
641 const bool pyramid_tri = lb.type() == BoundaryType::PYRAMID && lb[k] > 0;
642
643 if (lb.type() == BoundaryType::QUAD || prism_quad || pyramid_quad)
644 {
645 const auto map = [n_samples, size](int i, int j) { return j * n_samples + i + size; };
646
647 for (int j = 0; j < n_samples - 1; ++j)
648 {
649 for (int i = 0; i < n_samples - 1; ++i)
650 {
651 tris.emplace_back(map(i, j), map(i + 1, j), map(i, j + 1));
652 tris.emplace_back(map(i + 1, j + 1), map(i, j + 1), map(i + 1, j));
653 }
654 }
655 }
656 else if (lb.type() == BoundaryType::TRI || prism_tri || pyramid_tri)
657 {
658 int index = 0;
659 std::vector<int> mapp(n_samples * n_samples, -1);
660 for (int j = 0; j < n_samples; ++j)
661 {
662 for (int i = 0; i < n_samples - j; ++i)
663 {
664 mapp[j * n_samples + i] = index;
665 ++index;
666 }
667 }
668 const auto map = [mapp, n_samples](int i, int j) {
669 if (j * n_samples + i >= mapp.size())
670 return -1;
671 return mapp[j * n_samples + i];
672 };
673
674 for (int j = 0; j < n_samples - 1; ++j)
675 {
676 for (int i = 0; i < n_samples - j; ++i)
677 {
678 if (map(i, j) >= 0 && map(i + 1, j) >= 0 && map(i, j + 1) >= 0)
679 tris.emplace_back(map(i, j) + size, map(i + 1, j) + size, map(i, j + 1) + size);
680
681 if (map(i + 1, j + 1) >= 0 && map(i, j + 1) >= 0 && map(i + 1, j) >= 0)
682 tris.emplace_back(map(i + 1, j + 1) + size, map(i, j + 1) + size, map(i + 1, j) + size);
683 }
684 }
685 }
686 else
687 {
688 assert(false);
689 }
690 }
691 else
692 {
693 for (int i = 0; i < vertices.back().rows() - 1; ++i)
694 edges.emplace_back(i + size, i + size + 1);
695 }
696
697 normals.resize(vals.jac_it.size(), tmp_n.cols());
698 displaced_normals.resize(vals.jac_it.size(), tmp_n.cols());
699
700 for (int n = 0; n < vals.jac_it.size(); ++n)
701 {
702 trafo = vals.jac_it[n].inverse();
703
704 if (problem_dim == 2 || problem_dim == 3)
705 {
706
707 if (solution.size() > 0)
708 {
709 deform_mat.resize(problem_dim, problem_dim);
710 deform_mat.setZero();
711 for (const auto &b : vals.basis_values)
712 for (const auto &g : b.global)
713 for (int d = 0; d < problem_dim; ++d)
714 deform_mat.row(d) += solution(g.index * problem_dim + d) * b.grad.row(n);
715
716 trafo += deform_mat;
717 }
718 }
719
720 normals.row(n) = tmp_n * vals.jac_it[n];
721 normals.row(n).normalize();
722
723 displaced_normals.row(n) = tmp_n * trafo.inverse();
724 displaced_normals.row(n).normalize();
725 }
726
727 allnormals.push_back(normals);
728 displaced_allnormals.push_back(displaced_normals);
729
730 tmp_n.setZero();
731 for (int n = 0; n < vals.jac_it.size(); ++n)
732 {
733 tmp_n += normals.row(n);
734 }
735
736 if (mesh.is_volume())
737 {
738 Eigen::Vector3d e1 = vertices.back().row(std::get<1>(tris.back()) - size) - vertices.back().row(std::get<0>(tris.back()) - size);
739 Eigen::Vector3d e2 = vertices.back().row(std::get<2>(tris.back()) - size) - vertices.back().row(std::get<0>(tris.back()) - size);
740
741 Eigen::Vector3d n = e1.cross(e2);
742 Eigen::Vector3d nn = tmp_n.transpose();
743
744 if (n.dot(nn) < 0)
745 {
746 for (int i = tris_start; i < tris.size(); ++i)
747 {
748 tris[i] = std::tuple<int, int, int>(std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]));
749 }
750 }
751 }
752
753 size += vertices.back().rows();
754 }
755 }
756
757 boundary_vis_vertices.resize(size, vertices.front().cols());
758 boundary_vis_local_vertices.resize(size, vertices.front().cols());
759 boundary_vis_elements_ids.resize(size, 1);
760 boundary_vis_primitive_ids.resize(size, 1);
761 boundary_vis_normals.resize(size, vertices.front().cols());
762 displaced_boundary_vis_normals.resize(size, vertices.front().cols());
763
764 if (mesh.is_volume())
765 boundary_vis_elements.resize(tris.size(), 3);
766 else
767 boundary_vis_elements.resize(edges.size(), 2);
768
769 int index = 0;
770 int ii = 0;
771 for (const auto &v : vertices)
772 {
773 boundary_vis_vertices.block(index, 0, v.rows(), v.cols()) = v;
774 boundary_vis_local_vertices.block(index, 0, v.rows(), v.cols()) = lv[ii];
775 boundary_vis_elements_ids.block(index, 0, v.rows(), 1).setConstant(el_ids[ii]);
776 boundary_vis_primitive_ids.block(index, 0, v.rows(), 1).setConstant(global_primitive_ids[ii++]);
777 index += v.rows();
778 }
779
780 index = 0;
781 for (const auto &n : allnormals)
782 {
783 boundary_vis_normals.block(index, 0, n.rows(), n.cols()) = n;
784 index += n.rows();
785 }
786
787 index = 0;
788 for (const auto &n : displaced_allnormals)
789 {
790 displaced_boundary_vis_normals.block(index, 0, n.rows(), n.cols()) = n;
791 index += n.rows();
792 }
793
794 index = 0;
795 if (mesh.is_volume())
796 {
797 for (const auto &t : tris)
798 {
799 boundary_vis_elements.row(index) << std::get<0>(t), std::get<1>(t), std::get<2>(t);
800 ++index;
801 }
802 }
803 else
804 {
805 for (const auto &e : edges)
806 {
807 boundary_vis_elements.row(index) << e.first, e.second;
808 ++index;
809 }
810 }
811 }
812
813 void OutGeometryData::build_vis_mesh(
814 const mesh::Mesh &mesh,
815 const Eigen::VectorXi &disc_orders,
816 const std::vector<basis::ElementBases> &gbases,
817 const std::map<int, Eigen::MatrixXd> &polys,
818 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
819 const bool boundary_only,
820 Eigen::MatrixXd &points,
821 Eigen::MatrixXi &tets,
822 Eigen::MatrixXi &el_id,
823 Eigen::MatrixXd &discr) const
824 {
825 const auto &sampler = ref_element_sampler;
826
827 const auto &current_bases = gbases;
828 int tet_total_size = 0;
829 int pts_total_size = 0;
830
831 Eigen::MatrixXd vis_pts_poly;
832 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
833
834 for (size_t i = 0; i < current_bases.size(); ++i)
835 {
836 const auto &bs = current_bases[i];
837
838 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
839 continue;
840
841 if (mesh.is_simplex(i))
842 {
843 tet_total_size += sampler.simplex_volume().rows();
844 pts_total_size += sampler.simplex_points().rows();
845 }
846 else if (mesh.is_cube(i))
847 {
848 tet_total_size += sampler.cube_volume().rows();
849 pts_total_size += sampler.cube_points().rows();
850 }
851 else if (mesh.is_prism(i))
852 {
853 tet_total_size += sampler.prism_volume().rows();
854 pts_total_size += sampler.prism_points().rows();
855 }
856 else if (mesh.is_pyramid(i))
857 {
858 tet_total_size += sampler.pyramid_volume().rows();
859 pts_total_size += sampler.pyramid_points().rows();
860 }
861 else
862 {
863 if (mesh.is_volume())
864 {
865 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
866
867 tet_total_size += vis_faces_poly.rows();
868 pts_total_size += vis_pts_poly.rows();
869 }
870 else
871 {
872 sampler.sample_polygon(polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
873
874 tet_total_size += vis_faces_poly.rows();
875 pts_total_size += vis_pts_poly.rows();
876 }
877 }
878 }
879
880 points.resize(pts_total_size, mesh.dimension());
881 tets.resize(tet_total_size, mesh.is_volume() ? 4 : 3);
882
883 el_id.resize(pts_total_size, 1);
884 discr.resize(pts_total_size, 1);
885
886 Eigen::MatrixXd mapped, tmp;
887 int tet_index = 0, pts_index = 0;
888
889 for (size_t i = 0; i < current_bases.size(); ++i)
890 {
891 const auto &bs = current_bases[i];
892
893 if (boundary_only && mesh.is_volume() && !mesh.is_boundary_element(i))
894 continue;
895
896 if (mesh.is_simplex(i))
897 {
898 bs.eval_geom_mapping(sampler.simplex_points(), mapped);
899
900 tets.block(tet_index, 0, sampler.simplex_volume().rows(), tets.cols()) = sampler.simplex_volume().array() + pts_index;
901 tet_index += sampler.simplex_volume().rows();
902
903 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
904 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
905 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
906 pts_index += mapped.rows();
907 }
908 else if (mesh.is_cube(i))
909 {
910 bs.eval_geom_mapping(sampler.cube_points(), mapped);
911
912 tets.block(tet_index, 0, sampler.cube_volume().rows(), tets.cols()) = sampler.cube_volume().array() + pts_index;
913 tet_index += sampler.cube_volume().rows();
914
915 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
916 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
917 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
918 pts_index += mapped.rows();
919 }
920 else if (mesh.is_prism(i))
921 {
922 bs.eval_geom_mapping(sampler.prism_points(), mapped);
923
924 tets.block(tet_index, 0, sampler.prism_volume().rows(), tets.cols()) = sampler.prism_volume().array() + pts_index;
925 tet_index += sampler.prism_volume().rows();
926
927 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
928 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
929 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
930 pts_index += mapped.rows();
931 }
932 else if (mesh.is_pyramid(i))
933 {
934 bs.eval_geom_mapping(sampler.pyramid_points(), mapped);
935
936 tets.block(tet_index, 0, sampler.pyramid_volume().rows(), tets.cols()) = sampler.pyramid_volume().array() + pts_index;
937 tet_index += sampler.pyramid_volume().rows();
938
939 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
940 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
941 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
942 pts_index += mapped.rows();
943 }
944 else
945 {
946 if (mesh.is_volume())
947 {
948 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
949 bs.eval_geom_mapping(vis_pts_poly, mapped);
950
951 tets.block(tet_index, 0, vis_faces_poly.rows(), tets.cols()) = vis_faces_poly.array() + pts_index;
952 tet_index += vis_faces_poly.rows();
953
954 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
955 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(-1);
956 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
957 pts_index += mapped.rows();
958 }
959 else
960 {
961 sampler.sample_polygon(polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
962 bs.eval_geom_mapping(vis_pts_poly, mapped);
963
964 tets.block(tet_index, 0, vis_faces_poly.rows(), tets.cols()) = vis_faces_poly.array() + pts_index;
965 tet_index += vis_faces_poly.rows();
966
967 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
968 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(-1);
969 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
970 pts_index += mapped.rows();
971 }
972 }
973 }
974
975 assert(pts_index == points.rows());
976 assert(tet_index == tets.rows());
977 }
978
979 void OutGeometryData::build_high_order_vis_mesh(
980 const mesh::Mesh &mesh,
981 const Eigen::VectorXi &disc_orders,
982 const Eigen::VectorXi &disc_ordersq,
983 const std::vector<basis::ElementBases> &bases,
984 Eigen::MatrixXd &points,
985 std::vector<CellElement> &elements,
986 Eigen::MatrixXi &el_id,
987 Eigen::MatrixXd &discr) const
988 {
989 // if (!mesh)
990 // {
991 // logger().error("Load the mesh first!");
992 // return;
993 // }
994 // if (n_bases <= 0)
995 // {
996 // logger().error("Build the bases first!");
997 // return;
998 // }
999 // assert(mesh.is_linear());
1000
1001 std::vector<RowVectorNd> nodes;
1002 int pts_total_size = 0;
1003 elements.resize(bases.size());
1004 Eigen::MatrixXd ref_pts;
1005
1006 for (size_t i = 0; i < bases.size(); ++i)
1007 {
1008 const auto &bs = bases[i];
1009 if (mesh.is_volume())
1010 {
1011 if (mesh.is_simplex(i))
1012 autogen::p_nodes_3d(disc_orders(i), ref_pts);
1013 else if (mesh.is_cube(i))
1014 autogen::q_nodes_3d(disc_orders(i), ref_pts);
1015 else if (mesh.is_prism(i))
1016 {
1017 int max_order = std::max(disc_orders(i), disc_ordersq(i));
1018 autogen::prism_nodes_3d(max_order, max_order, ref_pts);
1019 }
1020 else if (mesh.is_pyramid(i))
1021 {
1022 autogen::pyramid_nodes_3d(disc_orders(i) == 2 ? -1 : disc_orders(i), ref_pts);
1023 }
1024 else
1025 continue;
1026 }
1027 else
1028 {
1029 if (mesh.is_simplex(i))
1030 autogen::p_nodes_2d(disc_orders(i), ref_pts);
1031 else if (mesh.is_cube(i))
1032 autogen::q_nodes_2d(disc_orders(i), ref_pts);
1033 else
1034 {
1035 const int n_v = static_cast<const mesh::Mesh2D &>(mesh).n_face_vertices(i);
1036 ref_pts.resize(n_v, 2);
1037 }
1038 }
1039
1040 pts_total_size += ref_pts.rows();
1041 }
1042
1043 points.resize(pts_total_size, mesh.dimension());
1044
1045 el_id.resize(pts_total_size, 1);
1046 discr.resize(pts_total_size, 1);
1047
1048 Eigen::MatrixXd mapped;
1049 int pts_index = 0;
1050
1051 std::string error_msg = "";
1052
1053 for (size_t i = 0; i < bases.size(); ++i)
1054 {
1055 const auto &bs = bases[i];
1056 if (mesh.is_volume())
1057 {
1058 if (mesh.is_simplex(i))
1059 autogen::p_nodes_3d(disc_orders(i), ref_pts);
1060 else if (mesh.is_cube(i))
1061 autogen::q_nodes_3d(disc_orders(i), ref_pts);
1062 else if (mesh.is_prism(i))
1063 {
1064 int max_order = std::max(disc_orders(i), disc_ordersq(i));
1065 autogen::prism_nodes_3d(max_order, max_order, ref_pts);
1066 }
1067 else if (mesh.is_pyramid(i))
1068 {
1069 autogen::pyramid_nodes_3d(disc_orders(i) == 2 ? -1 : disc_orders(i), ref_pts);
1070 }
1071 else
1072 continue;
1073 }
1074 else
1075 {
1076 if (mesh.is_simplex(i))
1077 autogen::p_nodes_2d(disc_orders(i), ref_pts);
1078 else if (mesh.is_cube(i))
1079 autogen::q_nodes_2d(disc_orders(i), ref_pts);
1080 else
1081 continue;
1082 }
1083
1084 bs.eval_geom_mapping(ref_pts, mapped);
1085
1086 for (int j = 0; j < mapped.rows(); ++j)
1087 {
1088 points.row(pts_index) = mapped.row(j);
1089 el_id(pts_index) = i;
1090 discr(pts_index) = disc_orders(i);
1091 elements[i].vertices.push_back(pts_index);
1092
1093 pts_index++;
1094 }
1095
1096 if (mesh.is_simplex(i))
1097 {
1098 if (mesh.is_volume())
1099 {
1100 const int n_nodes = elements[i].vertices.size();
1101 if (disc_orders(i) >= 3)
1102 {
1103 std::swap(elements[i].vertices[16], elements[i].vertices[17]);
1104 std::swap(elements[i].vertices[17], elements[i].vertices[18]);
1105 std::swap(elements[i].vertices[18], elements[i].vertices[19]);
1106 }
1107 if (disc_orders(i) > 4)
1108 error_msg = "Saving high-order meshes not implemented for P5+ elements!";
1109 }
1110 else
1111 {
1112 if (disc_orders(i) == 4)
1113 {
1114 const int n_nodes = elements[i].vertices.size();
1115 std::swap(elements[i].vertices[n_nodes - 1], elements[i].vertices[n_nodes - 2]);
1116 }
1117 if (disc_orders(i) > 4)
1118 error_msg = "Saving high-order meshes not implemented for P5+ elements!";
1119 }
1120 }
1121 else if (mesh.is_cube(i) && mesh.is_volume())
1122 {
1123 const int n_nodes = elements[i].vertices.size();
1124 if (disc_orders(i) == 2) // Lagrange hex, order=2
1125 {
1126 std::swap(elements[i].vertices[12], elements[i].vertices[16]);
1127 std::swap(elements[i].vertices[13], elements[i].vertices[17]);
1128 std::swap(elements[i].vertices[14], elements[i].vertices[18]);
1129 std::swap(elements[i].vertices[15], elements[i].vertices[19]);
1130 std::swap(elements[i].vertices[18], elements[i].vertices[19]); // a hack fix
1131 }
1132 // if (disc_orders(i) == 3) // Incomplete fix, need to fix order on the edge
1133 // {
1134 // std::swap(elements[i].vertices[24], elements[i].vertices[16]);
1135 // std::swap(elements[i].vertices[25], elements[i].vertices[17]);
1136 // std::swap(elements[i].vertices[26], elements[i].vertices[18]);
1137 // std::swap(elements[i].vertices[27], elements[i].vertices[19]);
1138 // std::swap(elements[i].vertices[28], elements[i].vertices[20]);
1139 // std::swap(elements[i].vertices[29], elements[i].vertices[21]);
1140 // std::swap(elements[i].vertices[30], elements[i].vertices[22]);
1141 // std::swap(elements[i].vertices[31], elements[i].vertices[23]);
1142 // std::swap(elements[i].vertices[28], elements[i].vertices[30]); // hack
1143 // std::swap(elements[i].vertices[29], elements[i].vertices[31]); // hack
1144 // }
1145 if (disc_orders(i) > 2)
1146 error_msg = "Saving high-order meshes not implemented for P2+ elements!";
1147 }
1148 else if (disc_orders(i) > 1)
1149 error_msg = "Saving high-order meshes not implemented for Q2+ elements!";
1150 }
1151
1152 if (!error_msg.empty())
1153 logger().warn(error_msg);
1154
1155 for (size_t i = 0; i < bases.size(); ++i)
1156 {
1157 if (mesh.is_volume() || !mesh.is_polytope(i))
1158 continue;
1159
1160 const auto &mesh2d = static_cast<const mesh::Mesh2D &>(mesh);
1161 const int n_v = mesh2d.n_face_vertices(i);
1162
1163 for (int j = 0; j < n_v; ++j)
1164 {
1165 points.row(pts_index) = mesh2d.point(mesh2d.face_vertex(i, j));
1166 el_id(pts_index) = i;
1167 discr(pts_index) = disc_orders(i);
1168 elements[i].vertices.push_back(pts_index);
1169
1170 pts_index++;
1171 }
1172 }
1173
1174 for (size_t i = 0; i < bases.size(); ++i)
1175 {
1176 if (!mesh.is_volume())
1177 {
1178 if (elements[i].vertices.size() == 1)
1179 elements[i].ctype = CellType::Vertex;
1180 else if (elements[i].vertices.size() == 2)
1181 elements[i].ctype = CellType::Line;
1182 else if (mesh.is_simplex(i))
1183 elements[i].ctype = CellType::Triangle;
1184 else if (mesh.is_cube(i))
1185 elements[i].ctype = CellType::Quadrilateral;
1186 else
1187 elements[i].ctype = CellType::Polygon;
1188 }
1189 else
1190 {
1191 if (mesh.is_simplex(i))
1192 elements[i].ctype = CellType::Tetrahedron;
1193 else if (mesh.is_cube(i))
1194 elements[i].ctype = CellType::Hexahedron;
1195 else if (mesh.is_prism(i))
1196 elements[i].ctype = CellType::Wedge;
1197 else if (mesh.is_pyramid(i))
1198 elements[i].ctype = CellType::Pyramid;
1199 }
1200 }
1201
1202 assert(pts_index == points.rows());
1203 }
1204
1205 void OutGeometryData::export_data(
1206 const State &state,
1207 const Eigen::MatrixXd &sol,
1208 const Eigen::MatrixXd &pressure,
1209 const bool is_time_dependent,
1210 const double tend_in,
1211 const double dt,
1212 const ExportOptions &opts,
1213 const std::string &vis_mesh_path,
1214 const std::string &nodes_path,
1215 const std::string &solution_path,
1216 const std::string &stress_path,
1217 const std::string &mises_path,
1218 const bool is_contact_enabled) const
1219 {
1220 if (!state.mesh)
1221 {
1222 logger().error("Load the mesh first!");
1223 return;
1224 }
1225 const int n_bases = state.n_bases;
1226 const std::vector<basis::ElementBases> &bases = state.bases;
1227 const std::vector<basis::ElementBases> &gbases = state.geom_bases();
1228 const mesh::Mesh &mesh = *state.mesh;
1229 const Eigen::VectorXi &in_node_to_node = state.in_node_to_node;
1230 const Eigen::MatrixXd &rhs = state.rhs;
1231 const assembler::Problem &problem = *state.problem;
1232
1233 if (n_bases <= 0)
1234 {
1235 logger().error("Build the bases first!");
1236 return;
1237 }
1238 // if (rhs.size() <= 0)
1239 // {
1240 // logger().error("Assemble the rhs first!");
1241 // return;
1242 // }
1243 if (sol.size() <= 0)
1244 {
1245 logger().error("Solve the problem first!");
1246 return;
1247 }
1248
1249 if (!solution_path.empty())
1250 {
1251 logger().info("Saving solution to {}", solution_path);
1252
1253 std::ofstream out(solution_path);
1254 out.precision(100);
1255 out << std::scientific;
1256 if (opts.reorder_output)
1257 {
1258 int problem_dim = (problem.is_scalar() ? 1 : mesh.dimension());
1259 Eigen::VectorXi reordering(n_bases);
1260 reordering.setConstant(-1);
1261
1262 for (int i = 0; i < in_node_to_node.size(); ++i)
1263 {
1264 reordering[in_node_to_node[i]] = i;
1265 }
1266 Eigen::MatrixXd tmp_sol = utils::unflatten(sol, problem_dim);
1267 Eigen::MatrixXd tmp(tmp_sol.rows(), tmp_sol.cols());
1268
1269 for (int i = 0; i < reordering.size(); ++i)
1270 {
1271 if (reordering[i] < 0)
1272 continue;
1273
1274 tmp.row(reordering[i]) = tmp_sol.row(i);
1275 }
1276
1277 for (int i = 0; i < tmp.rows(); ++i)
1278 {
1279 for (int j = 0; j < tmp.cols(); ++j)
1280 out << tmp(i, j) << " ";
1281
1282 out << std::endl;
1283 }
1284 }
1285 else
1286 out << sol << std::endl;
1287 out.close();
1288 }
1289
1290 double tend = tend_in;
1291 if (tend <= 0)
1292 tend = 1;
1293
1294 if (!vis_mesh_path.empty() && !is_time_dependent)
1295 {
1296 save_vtu(
1297 vis_mesh_path, state, sol, pressure,
1298 tend, dt, opts,
1299 is_contact_enabled);
1300 }
1301 if (!nodes_path.empty())
1302 {
1303 Eigen::MatrixXd nodes(n_bases, mesh.dimension());
1304 for (const basis::ElementBases &eb : bases)
1305 {
1306 for (const basis::Basis &b : eb.bases)
1307 {
1308 // for(const auto &lg : b.global())
1309 for (size_t ii = 0; ii < b.global().size(); ++ii)
1310 {
1311 const auto &lg = b.global()[ii];
1312 nodes.row(lg.index) = lg.node;
1313 }
1314 }
1315 }
1316 std::ofstream out(nodes_path);
1317 out.precision(100);
1318 out << nodes;
1319 out.close();
1320 }
1321 if (!stress_path.empty())
1322 {
1323 Eigen::MatrixXd result;
1324 Eigen::VectorXd mises;
1325 polyfem::io::Evaluator::compute_stress_at_quadrature_points(
1326 mesh, problem.is_scalar(),
1327 bases, gbases, state.disc_orders, state.disc_ordersq, *state.assembler,
1328 sol, tend, result, mises);
1329 std::ofstream out(stress_path);
1330 out.precision(20);
1331 out << result;
1332 }
1333 if (!mises_path.empty())
1334 {
1335 Eigen::MatrixXd result;
1336 Eigen::VectorXd mises;
1337 polyfem::io::Evaluator::compute_stress_at_quadrature_points(
1338 mesh, problem.is_scalar(),
1339 bases, gbases, state.disc_orders, state.disc_ordersq, *state.assembler,
1340 sol, tend, result, mises);
1341 std::ofstream out(mises_path);
1342 out.precision(20);
1343 out << mises;
1344 }
1345 }
1346
1347 bool OutGeometryData::ExportOptions::export_field(const std::string &field) const
1348 {
1349 return fields.empty() || std::find(fields.begin(), fields.end(), field) != fields.end();
1350 }
1351
1352 OutGeometryData::ExportOptions::ExportOptions(const json &args, const bool is_mesh_linear, const bool mesh_has_prisms, const bool is_problem_scalar)
1353 {
1354 fields = args["output"]["paraview"]["fields"];
1355
1356 volume = args["output"]["paraview"]["volume"];
1357 surface = args["output"]["paraview"]["surface"];
1358 wire = args["output"]["paraview"]["wireframe"];
1359 points = args["output"]["paraview"]["points"];
1360 contact_forces = args["output"]["paraview"]["options"]["contact_forces"] && !is_problem_scalar;
1361 friction_forces = args["output"]["paraview"]["options"]["friction_forces"] && !is_problem_scalar;
1362 normal_adhesion_forces = args["output"]["paraview"]["options"]["normal_adhesion_forces"] && !is_problem_scalar;
1363 tangential_adhesion_forces = args["output"]["paraview"]["options"]["tangential_adhesion_forces"] && !is_problem_scalar;
1364
1365 if (args["output"]["paraview"]["options"]["force_high_order"])
1366 use_sampler = false;
1367 else
1368 use_sampler = !(is_mesh_linear && args["output"]["paraview"]["high_order_mesh"]);
1369 boundary_only = use_sampler && args["output"]["advanced"]["vis_boundary_only"];
1370 material_params = args["output"]["paraview"]["options"]["material"];
1371 body_ids = args["output"]["paraview"]["options"]["body_ids"];
1372 sol_on_grid = args["output"]["advanced"]["sol_on_grid"] > 0;
1373 velocity = args["output"]["paraview"]["options"]["velocity"];
1374 acceleration = args["output"]["paraview"]["options"]["acceleration"];
1375 forces = args["output"]["paraview"]["options"]["forces"] && !is_problem_scalar;
1376 jacobian_validity = args["output"]["paraview"]["options"]["jacobian_validity"] && !is_problem_scalar;
1377
1378 scalar_values = args["output"]["paraview"]["options"]["scalar_values"];
1379 tensor_values = args["output"]["paraview"]["options"]["tensor_values"] && !is_problem_scalar;
1380 discretization_order = args["output"]["paraview"]["options"]["discretization_order"];
1381 nodes = args["output"]["paraview"]["options"]["nodes"] && !is_problem_scalar;
1382
1383 use_spline = args["space"]["basis_type"] == "Spline";
1384
1385 reorder_output = args["output"]["data"]["advanced"]["reorder_nodes"];
1386
1387 use_hdf5 = args["output"]["paraview"]["options"]["use_hdf5"];
1388 }
1389
1390 void OutGeometryData::save_vtu(
1391 const std::string &path,
1392 const State &state,
1393 const Eigen::MatrixXd &sol,
1394 const Eigen::MatrixXd &pressure,
1395 const double t,
1396 const double dt,
1397 const ExportOptions &opts,
1398 const bool is_contact_enabled) const
1399 {
1400 if (!state.mesh)
1401 {
1402 logger().error("Load the mesh first!");
1403 return;
1404 }
1405 const mesh::Mesh &mesh = *state.mesh;
1406 const Eigen::MatrixXd &rhs = state.rhs;
1407
1408 if (state.n_bases <= 0)
1409 {
1410 logger().error("Build the bases first!");
1411 return;
1412 }
1413 // if (rhs.size() <= 0)
1414 // {
1415 // logger().error("Assemble the rhs first!");
1416 // return;
1417 // }
1418 if (sol.size() <= 0)
1419 {
1420 logger().error("Solve the problem first!");
1421 return;
1422 }
1423
1424 const bool save_contact = is_contact_enabled && (opts.contact_forces || opts.friction_forces);
1425
1426 logger().info("Saving vtu to {}; volume={}, surface={}, contact={}, points={}, wireframe={}",
1427 path, opts.volume, opts.surface, save_contact, opts.points, opts.wire);
1428
1429 const std::filesystem::path fs_path(path);
1430 const std::string path_stem = fs_path.stem().string();
1431 const std::string base_path = (fs_path.parent_path() / path_stem).string();
1432
1433 if (opts.volume)
1434 {
1435 save_volume(base_path + opts.file_extension(), state, sol, pressure, t, dt, opts);
1436 }
1437
1438 if (opts.surface)
1439 {
1440 save_surface(base_path + "_surf" + opts.file_extension(), state, sol, pressure, t, dt, opts,
1441 is_contact_enabled);
1442 }
1443
1444 if (is_contact_enabled && (opts.contact_forces || opts.friction_forces || opts.normal_adhesion_forces || opts.tangential_adhesion_forces))
1445 {
1446 save_contact_surface(base_path + "_surf" + opts.file_extension(), state, sol, pressure, t, dt, opts,
1447 is_contact_enabled);
1448 }
1449
1450 if (opts.wire)
1451 {
1452 save_wire(base_path + "_wire" + opts.file_extension(), state, sol, t, opts);
1453 }
1454
1455 if (opts.points)
1456 {
1457 save_points(base_path + "_points" + opts.file_extension(), state, sol, opts);
1458 }
1459
1460 paraviewo::VTMWriter vtm(t);
1461 if (opts.volume)
1462 vtm.add_dataset("Volume", "data", path_stem + opts.file_extension());
1463 if (opts.surface)
1464 vtm.add_dataset("Surface", "data", path_stem + "_surf" + opts.file_extension());
1465 if (is_contact_enabled && (opts.contact_forces || opts.friction_forces || opts.normal_adhesion_forces || opts.tangential_adhesion_forces))
1466 vtm.add_dataset("Contact", "data", path_stem + "_surf_contact" + opts.file_extension());
1467 if (opts.wire)
1468 vtm.add_dataset("Wireframe", "data", path_stem + "_wire" + opts.file_extension());
1469 if (opts.points)
1470 vtm.add_dataset("Points", "data", path_stem + "_points" + opts.file_extension());
1471 vtm.save(base_path + ".vtm");
1472 }
1473
1474 void OutGeometryData::save_volume(
1475 const std::string &path,
1476 const State &state,
1477 const Eigen::MatrixXd &sol,
1478 const Eigen::MatrixXd &pressure,
1479 const double t,
1480 const double dt,
1481 const ExportOptions &opts) const
1482 {
1483 const Eigen::VectorXi &disc_orders = state.disc_orders;
1484 const Eigen::VectorXi &disc_ordersq = state.disc_ordersq;
1485 const auto &density = state.mass_matrix_assembler->density();
1486 const std::vector<basis::ElementBases> &bases = state.bases;
1487 const std::vector<basis::ElementBases> &pressure_bases = state.pressure_bases;
1488 const std::vector<basis::ElementBases> &gbases = state.geom_bases();
1489 const std::map<int, Eigen::MatrixXd> &polys = state.polys;
1490 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d = state.polys_3d;
1491 const assembler::Assembler &assembler = *state.assembler;
1492 const std::shared_ptr<time_integrator::ImplicitTimeIntegrator> time_integrator = state.solve_data.time_integrator;
1493 const mesh::Mesh &mesh = *state.mesh;
1494 const mesh::Obstacle &obstacle = state.obstacle;
1495 const assembler::Problem &problem = *state.problem;
1496
1497 Eigen::MatrixXd points;
1498 Eigen::MatrixXi tets;
1499 Eigen::MatrixXi el_id;
1500 Eigen::MatrixXd discr;
1501 std::vector<CellElement> elements;
1502
1503 if (opts.use_sampler)
1504 build_vis_mesh(mesh, disc_orders, gbases,
1505 state.polys, state.polys_3d, opts.boundary_only,
1506 points, tets, el_id, discr);
1507 else
1508 build_high_order_vis_mesh(mesh, disc_orders, disc_ordersq, bases,
1509 points, elements, el_id, discr);
1510
1511 Eigen::MatrixXd fun, exact_fun, err, node_fun;
1512
1513 if (opts.sol_on_grid)
1514 {
1515 const int problem_dim = problem.is_scalar() ? 1 : mesh.dimension();
1516 Eigen::MatrixXd tmp, tmp_grad;
1517 Eigen::MatrixXd tmp_p, tmp_grad_p;
1518 Eigen::MatrixXd res(grid_points_to_elements.size(), problem_dim);
1519 res.setConstant(std::numeric_limits<double>::quiet_NaN());
1520 Eigen::MatrixXd res_grad(grid_points_to_elements.size(), problem_dim * problem_dim);
1521 res_grad.setConstant(std::numeric_limits<double>::quiet_NaN());
1522
1523 Eigen::MatrixXd res_p(grid_points_to_elements.size(), 1);
1524 res_p.setConstant(std::numeric_limits<double>::quiet_NaN());
1525 Eigen::MatrixXd res_grad_p(grid_points_to_elements.size(), problem_dim);
1526 res_grad_p.setConstant(std::numeric_limits<double>::quiet_NaN());
1527
1528 for (int i = 0; i < grid_points_to_elements.size(); ++i)
1529 {
1530 const int el_id = grid_points_to_elements(i);
1531 if (el_id < 0)
1532 continue;
1533 assert(mesh.is_simplex(el_id));
1534 const Eigen::MatrixXd bc = grid_points_bc.row(i);
1535 Eigen::MatrixXd pt(1, bc.cols() - 1);
1536 for (int d = 1; d < bc.cols(); ++d)
1537 pt(d - 1) = bc(d);
1538 polyfem::io::Evaluator::interpolate_at_local_vals(
1539 mesh, problem.is_scalar(), bases, gbases,
1540 el_id, pt, sol, tmp, tmp_grad);
1541
1542 res.row(i) = tmp;
1543 res_grad.row(i) = tmp_grad;
1544
1545 if (state.mixed_assembler != nullptr)
1546 {
1547 polyfem::io::Evaluator::interpolate_at_local_vals(
1548 mesh, 1, pressure_bases, gbases,
1549 el_id, pt, pressure, tmp_p, tmp_grad_p);
1550 res_p.row(i) = tmp_p;
1551 res_grad_p.row(i) = tmp_grad_p;
1552 }
1553 }
1554
1555 std::ofstream os(path + "_sol.txt");
1556 os << res;
1557
1558 std::ofstream osg(path + "_grad.txt");
1559 osg << res_grad;
1560
1561 std::ofstream osgg(path + "_grid.txt");
1562 osgg << grid_points;
1563
1564 if (state.mixed_assembler != nullptr)
1565 {
1566 std::ofstream osp(path + "_p_sol.txt");
1567 osp << res_p;
1568
1569 std::ofstream osgp(path + "_p_grad.txt");
1570 osgp << res_grad_p;
1571 }
1572 }
1573
1574 Eigen::Vector<bool, -1> validity;
1575 if (opts.jacobian_validity)
1576 polyfem::io::Evaluator::mark_flipped_cells(
1577 mesh, gbases, bases, state.disc_orders,
1578 state.polys, state.polys_3d, ref_element_sampler,
1579 points.rows(), sol, validity, opts.use_sampler, opts.boundary_only);
1580
1581 polyfem::io::Evaluator::interpolate_function(
1582 mesh, problem.is_scalar(), bases, disc_orders, disc_ordersq,
1583 state.polys, state.polys_3d, ref_element_sampler,
1584 points.rows(), sol, fun, opts.use_sampler, opts.boundary_only);
1585
1586 {
1587 Eigen::MatrixXd tmp = Eigen::VectorXd::LinSpaced(sol.size(), 0, sol.size() - 1);
1588
1589 polyfem::io::Evaluator::interpolate_function(
1590 mesh, problem.is_scalar(), bases, disc_orders, disc_ordersq,
1591 state.polys, state.polys_3d, ref_element_sampler,
1592 points.rows(), tmp, node_fun, opts.use_sampler, opts.boundary_only);
1593 }
1594
1595 if (obstacle.n_vertices() > 0)
1596 {
1597 fun.conservativeResize(fun.rows() + obstacle.n_vertices(), fun.cols());
1598 node_fun.conservativeResize(node_fun.rows() + obstacle.n_vertices(), node_fun.cols());
1599 node_fun.bottomRows(obstacle.n_vertices()).setZero();
1600 // obstacle.update_displacement(t, fun);
1601 // NOTE: Assuming the obstacle displacement is the last part of the solution
1602 fun.bottomRows(obstacle.n_vertices()) = utils::unflatten(sol.bottomRows(obstacle.ndof()), fun.cols());
1603 }
1604
1605 if (problem.has_exact_sol())
1606 {
1607 problem.exact(points, t, exact_fun);
1608 err = (fun - exact_fun).eval().rowwise().norm();
1609
1610 if (obstacle.n_vertices() > 0)
1611 {
1612 exact_fun.conservativeResize(exact_fun.rows() + obstacle.n_vertices(), exact_fun.cols());
1613 // obstacle.update_displacement(t, exact_fun);
1614 exact_fun.bottomRows(obstacle.n_vertices()) = utils::unflatten(sol.bottomRows(obstacle.ndof()), fun.cols());
1615
1616 err.conservativeResize(err.rows() + obstacle.n_vertices(), 1);
1617 err.bottomRows(obstacle.n_vertices()).setZero();
1618 }
1619 }
1620
1621 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
1622 if (opts.use_hdf5)
1623 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
1624 else
1625 tmpw = std::make_shared<paraviewo::VTUWriter>();
1626 paraviewo::ParaviewWriter &writer = *tmpw;
1627
1628 if (validity.size() && opts.export_field("validity"))
1629 writer.add_field("validity", validity.cast<double>());
1630
1631 if (opts.nodes && opts.export_field("nodes"))
1632 writer.add_field("nodes", node_fun);
1633
1634 if (problem.is_time_dependent())
1635 {
1636 bool is_time_integrator_valid = time_integrator != nullptr;
1637
1638 if (opts.velocity || opts.export_field("velocity"))
1639 {
1640 const Eigen::VectorXd velocity =
1641 is_time_integrator_valid ? (time_integrator->v_prev()) : Eigen::VectorXd::Zero(sol.size());
1642 save_volume_vector_field(state, points, opts, "velocity", velocity, writer);
1643 }
1644
1645 if (opts.acceleration || opts.export_field("acceleration"))
1646 {
1647 const Eigen::VectorXd acceleration =
1648 is_time_integrator_valid ? (time_integrator->a_prev()) : Eigen::VectorXd::Zero(sol.size());
1649 save_volume_vector_field(state, points, opts, "acceleration", acceleration, writer);
1650 }
1651 }
1652
1653 if (opts.forces)
1654 {
1655 const double s = state.solve_data.time_integrator
1656 ? state.solve_data.time_integrator->acceleration_scaling()
1657 : 1;
1658
1659 for (const auto &[name, form] : state.solve_data.named_forms())
1660 {
1661 // NOTE: Assumes this form will be null for the entire sim
1662 if (form == nullptr)
1663 continue;
1664
1665 Eigen::VectorXd force;
1666 if (form->enabled())
1667 {
1668 form->first_derivative(sol, force);
1669 force *= -1.0 / s; // Divide by acceleration scaling to get units of force
1670 }
1671 else
1672 {
1673 force.setZero(sol.size());
1674 }
1675 if (opts.export_field(name + "_forces"))
1676 save_volume_vector_field(state, points, opts, name + "_forces", force, writer);
1677 }
1678 }
1679
1680 // if(problem->is_mixed())
1681 if (state.mixed_assembler != nullptr && opts.export_field("pressure"))
1682 {
1683 Eigen::MatrixXd interp_p;
1684 polyfem::io::Evaluator::interpolate_function(
1685 mesh, 1, // FIXME: disc_orders should use pressure discr orders, works only with sampler
1686 pressure_bases, disc_orders, disc_ordersq, state.polys, state.polys_3d, ref_element_sampler,
1687 points.rows(), pressure, interp_p, opts.use_sampler, opts.boundary_only);
1688
1689 if (obstacle.n_vertices() > 0)
1690 {
1691 interp_p.conservativeResize(interp_p.size() + obstacle.n_vertices(), 1);
1692 interp_p.bottomRows(obstacle.n_vertices()).setZero();
1693 }
1694
1695 writer.add_field("pressure", interp_p);
1696 }
1697
1698 if (obstacle.n_vertices() > 0)
1699 {
1700 discr.conservativeResize(discr.size() + obstacle.n_vertices(), 1);
1701 discr.bottomRows(obstacle.n_vertices()).setZero();
1702 }
1703
1704 if (opts.discretization_order && opts.export_field("discr"))
1705 writer.add_field("discr", discr);
1706
1707 if (problem.has_exact_sol())
1708 {
1709 if (opts.export_field("exact"))
1710 writer.add_field("exact", exact_fun);
1711 if (opts.export_field("error"))
1712 writer.add_field("error", err);
1713 }
1714
1715 if (fun.cols() != 1 && (opts.scalar_values || opts.tensor_values || (!opts.use_spline && (opts.scalar_values || opts.tensor_values))))
1716 {
1717 std::vector<assembler::Assembler::NamedMatrix> vals, tvals;
1718 polyfem::io::Evaluator::compute_scalar_value(
1719 mesh, problem.is_scalar(), bases, gbases,
1720 disc_orders, disc_ordersq, state.polys, state.polys_3d,
1721 *state.assembler,
1722 ref_element_sampler, points.rows(), sol, t, vals, opts.use_sampler, opts.boundary_only);
1723
1724 for (auto &[_, v] : vals)
1725 utils::append_rows_of_zeros(v, obstacle.n_vertices());
1726
1727 if (opts.scalar_values)
1728 {
1729 for (const auto &[name, v] : vals)
1730 {
1731 if (opts.export_field(name))
1732 writer.add_field(name, v);
1733 }
1734 }
1735
1736 if (opts.tensor_values)
1737 {
1738 polyfem::io::Evaluator::compute_tensor_value(
1739 mesh, problem.is_scalar(), bases, gbases, disc_orders, disc_ordersq,
1740 state.polys, state.polys_3d, *state.assembler, ref_element_sampler,
1741 points.rows(), sol, t, tvals, opts.use_sampler, opts.boundary_only);
1742
1743 for (auto &[_, v] : tvals)
1744 utils::append_rows_of_zeros(v, obstacle.n_vertices());
1745
1746 for (const auto &[name, v] : tvals)
1747 {
1748 const int stride = mesh.dimension();
1749 assert(v.cols() % stride == 0);
1750
1751 if (!opts.export_field(name))
1752 continue;
1753
1754 for (int i = 0; i < v.cols(); i += stride)
1755 {
1756 const Eigen::MatrixXd tmp = v.middleCols(i, stride);
1757 assert(tmp.cols() == stride);
1758
1759 const int ii = (i / stride) + 1;
1760 writer.add_field(fmt::format("{:s}_{:d}", name, ii), tmp);
1761 }
1762 }
1763 }
1764
1765 if (!opts.use_spline && (opts.scalar_values || opts.tensor_values))
1766 {
1767 polyfem::io::Evaluator::average_grad_based_function(
1768 mesh, problem.is_scalar(), state.n_bases, bases, gbases,
1769 disc_orders, disc_ordersq, state.polys, state.polys_3d, *state.assembler,
1770 ref_element_sampler, t, points.rows(), sol, vals, tvals,
1771 opts.use_sampler, opts.boundary_only);
1772
1773 if (obstacle.n_vertices() > 0)
1774 {
1775 for (auto &v : vals)
1776 {
1777 v.second.conservativeResize(v.second.size() + obstacle.n_vertices(), 1);
1778 v.second.bottomRows(obstacle.n_vertices()).setZero();
1779 }
1780 }
1781
1782 if (opts.scalar_values)
1783 {
1784 for (const auto &v : vals)
1785 {
1786 if (opts.export_field(fmt::format("{:s}_avg", v.first)))
1787 writer.add_field(fmt::format("{:s}_avg", v.first), v.second);
1788 }
1789 }
1790 if (opts.tensor_values)
1791 {
1792 for (const auto &v : tvals)
1793 {
1794 const int stride = mesh.dimension();
1795 assert(v.second.cols() % stride == 0);
1796
1797 if (!opts.export_field(fmt::format("{:s}_avg", v.first)))
1798 continue;
1799
1800 for (int i = 0; i < v.second.cols(); i += stride)
1801 {
1802 const Eigen::MatrixXd tmp = v.second.middleCols(i, stride);
1803 assert(tmp.cols() == stride);
1804
1805 const int ii = (i / stride) + 1;
1806 writer.add_field(
1807 fmt::format("{:s}_avg_{:d}", v.first, ii), tmp);
1808 }
1809 }
1810 }
1811 }
1812 }
1813
1814 if (opts.material_params)
1815 {
1816 const auto &params = assembler.parameters();
1817
1818 std::map<std::string, Eigen::MatrixXd> param_val;
1819 for (const auto &[p, _] : params)
1820 param_val[p] = Eigen::MatrixXd(points.rows(), 1);
1821 Eigen::MatrixXd rhos(points.rows(), 1);
1822
1823 Eigen::MatrixXd local_pts;
1824 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
1825
1826 int index = 0;
1827 const auto &sampler = ref_element_sampler;
1828 for (int e = 0; e < int(bases.size()); ++e)
1829 {
1830 const basis::ElementBases &gbs = gbases[e];
1831 const basis::ElementBases &bs = bases[e];
1832
1833 if (opts.use_sampler)
1834 {
1835 if (mesh.is_simplex(e))
1836 local_pts = sampler.simplex_points();
1837 else if (mesh.is_cube(e))
1838 local_pts = sampler.cube_points();
1839 else if (mesh.is_prism(e))
1840 local_pts = sampler.prism_points();
1841 else if (mesh.is_pyramid(e))
1842 local_pts = sampler.pyramid_points();
1843 else
1844 {
1845 if (mesh.is_volume())
1846 sampler.sample_polyhedron(polys_3d.at(e).first, polys_3d.at(e).second, local_pts, vis_faces_poly, vis_edges_poly);
1847 else
1848 sampler.sample_polygon(polys.at(e), local_pts, vis_faces_poly, vis_edges_poly);
1849 }
1850 }
1851 else
1852 {
1853 if (mesh.is_volume())
1854 {
1855 if (mesh.is_simplex(e))
1856 autogen::p_nodes_3d(disc_orders(e), local_pts);
1857 else if (mesh.is_cube(e))
1858 autogen::q_nodes_3d(disc_orders(e), local_pts);
1859 else if (mesh.is_prism(e))
1860 {
1861 const auto o = std::max(disc_orders(e), disc_ordersq(e));
1862 autogen::prism_nodes_3d(o, o, local_pts);
1863 }
1864 else if (mesh.is_pyramid(e))
1865 autogen::pyramid_nodes_3d(disc_orders(e) == 2 ? -1 : disc_orders(e), local_pts);
1866
1867 else
1868 continue;
1869 }
1870 else
1871 {
1872 if (mesh.is_simplex(e))
1873 autogen::p_nodes_2d(disc_orders(e), local_pts);
1874 else if (mesh.is_cube(e))
1875 autogen::q_nodes_2d(disc_orders(e), local_pts);
1876 else
1877 {
1878 const auto &mesh2d = static_cast<const mesh::Mesh2D &>(mesh);
1879 const int n_v = mesh2d.n_face_vertices(e);
1880 local_pts.resize(n_v, 2);
1881
1882 for (int j = 0; j < n_v; ++j)
1883 {
1884 local_pts.row(j) = mesh2d.point(mesh2d.face_vertex(e, j));
1885 }
1886 }
1887 }
1888 }
1889
1890 assembler::ElementAssemblyValues vals;
1891 vals.compute(e, mesh.is_volume(), local_pts, bs, gbs);
1892
1893 for (int j = 0; j < vals.val.rows(); ++j)
1894 {
1895 for (const auto &[p, func] : params)
1896 param_val.at(p)(index) = func(local_pts.row(j), vals.val.row(j), t, e);
1897
1898 rhos(index) = density(local_pts.row(j), vals.val.row(j), t, e);
1899
1900 ++index;
1901 }
1902 }
1903
1904 assert(index == points.rows());
1905
1906 if (obstacle.n_vertices() > 0)
1907 {
1908 for (auto &[_, tmp] : param_val)
1909 {
1910 tmp.conservativeResize(tmp.size() + obstacle.n_vertices(), 1);
1911 tmp.bottomRows(obstacle.n_vertices()).setZero();
1912 }
1913
1914 rhos.conservativeResize(rhos.size() + obstacle.n_vertices(), 1);
1915 rhos.bottomRows(obstacle.n_vertices()).setZero();
1916 }
1917 for (const auto &[p, tmp] : param_val)
1918 {
1919 if (opts.export_field(p))
1920 writer.add_field(p, tmp);
1921 }
1922 if (opts.export_field("rho"))
1923 writer.add_field("rho", rhos);
1924 }
1925
1926 if (opts.body_ids || opts.export_field("body_ids"))
1927 {
1928
1929 Eigen::MatrixXd ids(points.rows(), 1);
1930
1931 for (int i = 0; i < points.rows(); ++i)
1932 {
1933 ids(i) = mesh.get_body_id(el_id(i));
1934 }
1935
1936 if (obstacle.n_vertices() > 0)
1937 {
1938 ids.conservativeResize(ids.size() + obstacle.n_vertices(), 1);
1939 ids.bottomRows(obstacle.n_vertices()).setZero();
1940 }
1941
1942 writer.add_field("body_ids", ids);
1943 }
1944
1945 // if (opts.export_field("rhs"))
1946 // {
1947 // interpolate_function(pts_index, rhs, fun, opts.boundary_only);
1948 // writer.add_field("rhs", fun);
1949 // }
1950
1951 if (fun.cols() != 1 && state.mixed_assembler == nullptr && opts.export_field("traction_force"))
1952 {
1953 Eigen::MatrixXd traction_forces, traction_forces_fun;
1954 compute_traction_forces(state, sol, t, traction_forces, false);
1955
1956 polyfem::io::Evaluator::interpolate_function(
1957 mesh, problem.is_scalar(), bases, disc_orders, disc_ordersq,
1958 state.polys, state.polys_3d, ref_element_sampler,
1959 points.rows(), traction_forces, traction_forces_fun, opts.use_sampler, opts.boundary_only);
1960
1961 if (obstacle.n_vertices() > 0)
1962 {
1963 traction_forces_fun.conservativeResize(traction_forces_fun.rows() + obstacle.n_vertices(), traction_forces_fun.cols());
1964 traction_forces_fun.bottomRows(obstacle.n_vertices()).setZero();
1965 }
1966
1967 writer.add_field("traction_force", traction_forces_fun);
1968 }
1969
1970 if (fun.cols() != 1 && state.mixed_assembler == nullptr && opts.export_field("gradient_of_elastic_potential"))
1971 {
1972 try
1973 {
1974 Eigen::VectorXd potential_grad;
1975 Eigen::MatrixXd potential_grad_fun;
1976 if (state.solve_data.elastic_form)
1977 state.solve_data.elastic_form->first_derivative(sol, potential_grad);
1978
1979 polyfem::io::Evaluator::interpolate_function(
1980 mesh, problem.is_scalar(), bases, disc_orders, disc_ordersq,
1981 state.polys, state.polys_3d, ref_element_sampler,
1982 points.rows(), potential_grad, potential_grad_fun, opts.use_sampler, opts.boundary_only);
1983
1984 if (obstacle.n_vertices() > 0)
1985 {
1986 potential_grad_fun.conservativeResize(potential_grad_fun.rows() + obstacle.n_vertices(), potential_grad_fun.cols());
1987 potential_grad_fun.bottomRows(obstacle.n_vertices()).setZero();
1988 }
1989
1990 writer.add_field("gradient_of_elastic_potential", potential_grad_fun);
1991 }
1992 catch (std::exception &)
1993 {
1994 }
1995 }
1996
1997 if (fun.cols() != 1 && state.mixed_assembler == nullptr && opts.export_field("gradient_of_contact_potential"))
1998 {
1999 try
2000 {
2001 Eigen::VectorXd potential_grad;
2002 Eigen::MatrixXd potential_grad_fun;
2003 if (state.solve_data.contact_form && state.solve_data.contact_form->weight() > 0)
2004 {
2005 state.solve_data.contact_form->first_derivative(sol, potential_grad);
2006 potential_grad *= -state.solve_data.contact_form->barrier_stiffness() / state.solve_data.contact_form->weight();
2007
2008 polyfem::io::Evaluator::interpolate_function(
2009 mesh, problem.is_scalar(), bases, state.disc_orders, state.disc_ordersq,
2010 state.polys, state.polys_3d, ref_element_sampler,
2011 points.rows(), potential_grad, potential_grad_fun, opts.use_sampler, opts.boundary_only);
2012
2013 if (obstacle.n_vertices() > 0)
2014 {
2015 potential_grad_fun.conservativeResize(potential_grad_fun.rows() + obstacle.n_vertices(), potential_grad_fun.cols());
2016 potential_grad_fun.bottomRows(obstacle.n_vertices()).setZero();
2017 }
2018
2019 writer.add_field("gradient_of_contact_potential", potential_grad_fun);
2020 }
2021 }
2022 catch (std::exception &)
2023 {
2024 }
2025 }
2026
2027 // Write the solution last so it is the default for warp-by-vector
2028 writer.add_field("solution", fun);
2029
2030 if (obstacle.n_vertices() > 0)
2031 {
2032 const int orig_p = points.rows();
2033 points.conservativeResize(points.rows() + obstacle.n_vertices(), points.cols());
2034 points.bottomRows(obstacle.n_vertices()) = obstacle.v();
2035
2036 if (elements.empty())
2037 {
2038 for (int i = 0; i < tets.rows(); ++i)
2039 {
2040 elements.emplace_back();
2041 elements.back().ctype = CellType::Tetrahedron;
2042 for (int j = 0; j < tets.cols(); ++j)
2043 elements.back().vertices.push_back(tets(i, j));
2044 }
2045 }
2046
2047 for (int i = 0; i < obstacle.get_face_connectivity().rows(); ++i)
2048 {
2049 elements.emplace_back();
2050 elements.back().ctype = CellType::Triangle;
2051 for (int j = 0; j < obstacle.get_face_connectivity().cols(); ++j)
2052 elements.back().vertices.push_back(obstacle.get_face_connectivity()(i, j) + orig_p);
2053 }
2054
2055 for (int i = 0; i < obstacle.get_edge_connectivity().rows(); ++i)
2056 {
2057 elements.emplace_back();
2058 elements.back().ctype = CellType::Line;
2059 for (int j = 0; j < obstacle.get_edge_connectivity().cols(); ++j)
2060 elements.back().vertices.push_back(obstacle.get_edge_connectivity()(i, j) + orig_p);
2061 }
2062
2063 for (int i = 0; i < obstacle.get_vertex_connectivity().size(); ++i)
2064 {
2065 elements.emplace_back();
2066 elements.back().ctype = CellType::Vertex;
2067 elements.back().vertices.push_back(obstacle.get_vertex_connectivity()(i) + orig_p);
2068 }
2069 }
2070
2071 if (elements.empty())
2072 writer.write_mesh(path, points, tets, mesh.is_volume() ? CellType::Tetrahedron : CellType::Triangle);
2073 else
2074 writer.write_mesh(path, points, elements);
2075 }
2076
2077 void OutGeometryData::save_volume_vector_field(
2078 const State &state,
2079 const Eigen::MatrixXd &points,
2080 const ExportOptions &opts,
2081 const std::string &name,
2082 const Eigen::VectorXd &field,
2083 paraviewo::ParaviewWriter &writer) const
2084 {
2085 Eigen::MatrixXd inerpolated_field;
2086 polyfem::io::Evaluator::interpolate_function(
2087 *state.mesh, state.problem->is_scalar(), state.bases, state.disc_orders, state.disc_ordersq,
2088 state.polys, state.polys_3d, ref_element_sampler,
2089 points.rows(), field, inerpolated_field, opts.use_sampler, opts.boundary_only);
2090
2091 if (state.obstacle.n_vertices() > 0)
2092 {
2093 inerpolated_field.conservativeResize(
2094 inerpolated_field.rows() + state.obstacle.n_vertices(), inerpolated_field.cols());
2095 inerpolated_field.bottomRows(state.obstacle.n_vertices()) =
2096 utils::unflatten(field.tail(state.obstacle.ndof()), inerpolated_field.cols());
2097 }
2098 if (opts.export_field(name))
2099 writer.add_field(name, inerpolated_field);
2100 }
2101
2102 void OutGeometryData::save_surface(
2103 const std::string &export_surface,
2104 const State &state,
2105 const Eigen::MatrixXd &sol,
2106 const Eigen::MatrixXd &pressure,
2107 const double t,
2108 const double dt_in,
2109 const ExportOptions &opts,
2110 const bool is_contact_enabled) const
2111 {
2112
2113 const Eigen::VectorXi &disc_orders = state.disc_orders;
2114 const auto &density = state.mass_matrix_assembler->density();
2115 const std::vector<basis::ElementBases> &bases = state.bases;
2116 const std::vector<basis::ElementBases> &pressure_bases = state.pressure_bases;
2117 const std::vector<basis::ElementBases> &gbases = state.geom_bases();
2118 const assembler::Assembler &assembler = *state.assembler;
2119 const assembler::Problem &problem = *state.problem;
2120 const mesh::Mesh &mesh = *state.mesh;
2121 int problem_dim = (problem.is_scalar() ? 1 : mesh.dimension());
2122
2123 Eigen::MatrixXd boundary_vis_vertices;
2124 Eigen::MatrixXd boundary_vis_local_vertices;
2125 Eigen::MatrixXi boundary_vis_elements;
2126 Eigen::MatrixXi boundary_vis_elements_ids;
2127 Eigen::MatrixXi boundary_vis_primitive_ids;
2128 Eigen::MatrixXd boundary_vis_normals;
2129 Eigen::MatrixXd displaced_boundary_vis_normals;
2130
2131 build_vis_boundary_mesh(mesh, bases, gbases, state.total_local_boundary, sol, problem_dim,
2132 boundary_vis_vertices, boundary_vis_local_vertices, boundary_vis_elements,
2133 boundary_vis_elements_ids, boundary_vis_primitive_ids, boundary_vis_normals,
2134 displaced_boundary_vis_normals);
2135
2136 Eigen::MatrixXd fun, interp_p, discr, vect, b_sidesets;
2137
2138 Eigen::MatrixXd lsol, lp, lgrad, lpgrad;
2139
2140 int actual_dim = 1;
2141 if (!problem.is_scalar())
2142 actual_dim = mesh.dimension();
2143
2144 discr.resize(boundary_vis_vertices.rows(), 1);
2145 fun.resize(boundary_vis_vertices.rows(), actual_dim);
2146 interp_p.resize(boundary_vis_vertices.rows(), 1);
2147 vect.resize(boundary_vis_vertices.rows(), mesh.dimension());
2148
2149 b_sidesets.resize(boundary_vis_vertices.rows(), 1);
2150 b_sidesets.setZero();
2151
2152 for (int i = 0; i < boundary_vis_vertices.rows(); ++i)
2153 {
2154 const auto s_id = mesh.get_boundary_id(boundary_vis_primitive_ids(i));
2155 if (s_id > 0)
2156 {
2157 b_sidesets(i) = s_id;
2158 }
2159
2160 const int el_index = boundary_vis_elements_ids(i);
2161 polyfem::io::Evaluator::interpolate_at_local_vals(
2162 mesh, problem.is_scalar(), bases, gbases,
2163 el_index, boundary_vis_local_vertices.row(i), sol, lsol, lgrad);
2164 assert(lsol.size() == actual_dim);
2165 if (state.mixed_assembler != nullptr)
2166 {
2167 polyfem::io::Evaluator::interpolate_at_local_vals(
2168 mesh, 1, pressure_bases, gbases,
2169 el_index, boundary_vis_local_vertices.row(i), pressure, lp, lpgrad);
2170 assert(lp.size() == 1);
2171 interp_p(i) = lp(0);
2172 }
2173
2174 discr(i) = disc_orders(el_index);
2175 for (int j = 0; j < actual_dim; ++j)
2176 {
2177 fun(i, j) = lsol(j);
2178 }
2179
2180 if (actual_dim == 1)
2181 {
2182 assert(lgrad.size() == mesh.dimension());
2183 for (int j = 0; j < mesh.dimension(); ++j)
2184 {
2185 vect(i, j) = lgrad(j);
2186 }
2187 }
2188 else
2189 {
2190 assert(lgrad.size() == actual_dim * actual_dim);
2191 std::vector<assembler::Assembler::NamedMatrix> tensor_flat;
2192 const basis::ElementBases &gbs = gbases[el_index];
2193 const basis::ElementBases &bs = bases[el_index];
2194 assembler.compute_tensor_value(assembler::OutputData(t, el_index, bs, gbs, boundary_vis_local_vertices.row(i), sol), tensor_flat);
2195 // TF computed only from cauchy stress
2196 assert(tensor_flat[0].first == "cauchy_stess");
2197 assert(tensor_flat[0].second.size() == actual_dim * actual_dim);
2198
2199 Eigen::Map<Eigen::MatrixXd> tensor(tensor_flat[0].second.data(), actual_dim, actual_dim);
2200 vect.row(i) = displaced_boundary_vis_normals.row(i) * tensor;
2201
2202 double area = 0;
2203 if (mesh.is_volume())
2204 {
2205 if (mesh.is_simplex(el_index))
2206 area = mesh.tri_area(boundary_vis_primitive_ids(i));
2207 else if (mesh.is_cube(el_index))
2208 area = mesh.quad_area(boundary_vis_primitive_ids(i));
2209 else if (mesh.is_prism(el_index))
2210 {
2211 const int tmp = boundary_vis_primitive_ids(i);
2212 area = mesh.n_face_vertices(tmp) == 4 ? mesh.quad_area(tmp) : mesh.tri_area(tmp);
2213 }
2214 else if (mesh.is_pyramid(el_index))
2215 {
2216 const int tmp = boundary_vis_primitive_ids(i);
2217 area = mesh.n_face_vertices(tmp) == 4 ? mesh.quad_area(tmp) : mesh.tri_area(tmp);
2218 }
2219 }
2220 else
2221 area = mesh.edge_length(boundary_vis_primitive_ids(i));
2222
2223 vect.row(i) *= area;
2224 }
2225 }
2226
2227 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2228 if (opts.use_hdf5)
2229 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2230 else
2231 tmpw = std::make_shared<paraviewo::VTUWriter>();
2232 paraviewo::ParaviewWriter &writer = *tmpw;
2233
2234 if (opts.export_field("normals"))
2235 writer.add_field("normals", boundary_vis_normals);
2236 if (opts.export_field("displaced_normals"))
2237 writer.add_field("displaced_normals", displaced_boundary_vis_normals);
2238 if (state.mixed_assembler != nullptr && opts.export_field("pressure"))
2239 writer.add_field("pressure", interp_p);
2240 if (opts.export_field("discr"))
2241 writer.add_field("discr", discr);
2242 if (opts.export_field("sidesets"))
2243 writer.add_field("sidesets", b_sidesets);
2244
2245 if (actual_dim == 1 && opts.export_field("solution_grad"))
2246 writer.add_field("solution_grad", vect);
2247 else if (opts.export_field("traction_force"))
2248 {
2249 writer.add_field("traction_force", vect);
2250 }
2251
2252 if (opts.material_params)
2253 {
2254 const auto &params = assembler.parameters();
2255
2256 std::map<std::string, Eigen::MatrixXd> param_val;
2257 for (const auto &[p, _] : params)
2258 param_val[p] = Eigen::MatrixXd(boundary_vis_vertices.rows(), 1);
2259 Eigen::MatrixXd rhos(boundary_vis_vertices.rows(), 1);
2260
2261 for (int i = 0; i < boundary_vis_vertices.rows(); ++i)
2262 {
2263 double lambda, mu;
2264
2265 for (const auto &[p, func] : params)
2266 param_val.at(p)(i) = func(boundary_vis_local_vertices.row(i), boundary_vis_vertices.row(i), t, boundary_vis_elements_ids(i));
2267
2268 rhos(i) = density(boundary_vis_local_vertices.row(i), boundary_vis_vertices.row(i), t, boundary_vis_elements_ids(i));
2269 }
2270
2271 for (const auto &[p, tmp] : param_val)
2272 {
2273 if (opts.export_field(p))
2274 writer.add_field(p, tmp);
2275 }
2276 if (opts.export_field("rho"))
2277 writer.add_field("rho", rhos);
2278 }
2279
2280 if (opts.body_ids || opts.export_field("body_ids"))
2281 {
2282
2283 Eigen::MatrixXd ids(boundary_vis_vertices.rows(), 1);
2284
2285 for (int i = 0; i < boundary_vis_vertices.rows(); ++i)
2286 {
2287 ids(i) = mesh.get_body_id(boundary_vis_elements_ids(i));
2288 }
2289
2290 writer.add_field("body_ids", ids);
2291 }
2292
2293 // Write the solution last so it is the default for warp-by-vector
2294 writer.add_field("solution", fun);
2295 writer.write_mesh(export_surface, boundary_vis_vertices, boundary_vis_elements, mesh.is_volume() ? CellType::Triangle : CellType::Line);
2296 }
2297
2298 void OutGeometryData::save_contact_surface(
2299 const std::string &export_surface,
2300 const State &state,
2301 const Eigen::MatrixXd &sol,
2302 const Eigen::MatrixXd &pressure,
2303 const double t,
2304 const double dt_in,
2305 const ExportOptions &opts,
2306 const bool is_contact_enabled) const
2307 {
2308 const mesh::Mesh &mesh = *state.mesh;
2309 const ipc::CollisionMesh &collision_mesh = state.collision_mesh;
2310 const double dhat = state.args["contact"]["dhat"];
2311 const double friction_coefficient = state.args["contact"]["friction_coefficient"];
2312 const double epsv = state.args["contact"]["epsv"];
2313 const double dhat_a = state.args["contact"]["adhesion"]["dhat_a"];
2314 const double dhat_p = state.args["contact"]["adhesion"]["dhat_p"];
2315 const double Y = state.args["contact"]["adhesion"]["adhesion_strength"];
2316 const double epsa = state.args["contact"]["adhesion"]["epsa"];
2317 const double tangential_adhesion_coefficient = state.args["contact"]["adhesion"]["tangential_adhesion_coefficient"];
2318 const std::shared_ptr<solver::ContactForm> &contact_form = state.solve_data.contact_form;
2319 const std::shared_ptr<solver::FrictionForm> &friction_form = state.solve_data.friction_form;
2320 const std::shared_ptr<solver::NormalAdhesionForm> &normal_adhesion_form = state.solve_data.normal_adhesion_form;
2321 const std::shared_ptr<solver::TangentialAdhesionForm> &tangential_adhesion_form = state.solve_data.tangential_adhesion_form;
2322
2323 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2324 if (opts.use_hdf5)
2325 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2326 else
2327 tmpw = std::make_shared<paraviewo::VTUWriter>();
2328 paraviewo::ParaviewWriter &writer = *tmpw;
2329
2330 const int problem_dim = mesh.dimension();
2331 const Eigen::MatrixXd full_displacements = utils::unflatten(sol, problem_dim);
2332 const Eigen::MatrixXd surface_displacements = collision_mesh.map_displacements(full_displacements);
2333
2334 const Eigen::MatrixXd displaced_surface = collision_mesh.displace_vertices(full_displacements);
2335
2336 ipc::NormalCollisions collision_set;
2337 // collision_set.set_use_convergent_formulation(state.args["contact"]["use_convergent_formulation"]);
2338 if (state.args["contact"]["use_convergent_formulation"])
2339 {
2340 collision_set.set_use_area_weighting(state.args["contact"]["use_area_weighting"]);
2341 collision_set.set_use_improved_max_approximator(state.args["contact"]["use_improved_max_operator"]);
2342 }
2343
2344 collision_set.build(
2345 collision_mesh, displaced_surface, dhat,
2346 /*dmin=*/0, ipc::create_broad_phase(state.args["solver"]["contact"]["CCD"]["broad_phase"]).get());
2347
2348 const double barrier_stiffness = contact_form != nullptr ? contact_form->barrier_stiffness() : 1;
2349 ipc::BarrierPotential barrier_potential(dhat, barrier_stiffness);
2350 if (state.args["contact"]["use_convergent_formulation"])
2351 {
2352 barrier_potential.set_use_physical_barrier(state.args["contact"]["use_physical_barrier"]);
2353 }
2354
2355 if (opts.contact_forces || opts.export_field("contact_forces"))
2356 {
2357 Eigen::MatrixXd forces = -barrier_potential.gradient(collision_set, collision_mesh, displaced_surface);
2358
2359 Eigen::MatrixXd forces_reshaped = utils::unflatten(forces, problem_dim);
2360
2361 assert(forces_reshaped.rows() == surface_displacements.rows());
2362 assert(forces_reshaped.cols() == surface_displacements.cols());
2363 writer.add_field("contact_forces", forces_reshaped);
2364 }
2365
2366 if (contact_form && state.args["contact"]["use_gcp_formulation"] && state.args["contact"]["use_adaptive_dhat"] && opts.export_field("adaptive_dhat"))
2367 {
2368 const auto form = std::dynamic_pointer_cast<solver::SmoothContactForm>(contact_form);
2369 assert(form);
2370 const auto &set = form->collision_set();
2371
2372 if (problem_dim == 2)
2373 {
2374 Eigen::VectorXd dhats(collision_mesh.num_edges());
2375 dhats.setConstant(dhat);
2376 for (int e = 0; e < dhats.size(); e++)
2377 dhats(e) = set.get_edge_dhat(e);
2378
2379 writer.add_cell_field("dhat", dhats);
2380 }
2381 else
2382 {
2383 Eigen::VectorXd fdhats(collision_mesh.num_faces());
2384 fdhats.setConstant(dhat);
2385 for (int e = 0; e < fdhats.size(); e++)
2386 fdhats(e) = set.get_face_dhat(e);
2387
2388 writer.add_cell_field("dhat_face", fdhats);
2389
2390 Eigen::VectorXd vdhats(collision_mesh.num_vertices());
2391 vdhats.setConstant(dhat);
2392 for (int i = 0; i < vdhats.size(); i++)
2393 vdhats(i) = set.get_vert_dhat(i);
2394
2395 writer.add_field("dhat_vert", vdhats);
2396 }
2397 }
2398
2399 if (opts.friction_forces || opts.export_field("friction_forces"))
2400 {
2401 ipc::TangentialCollisions friction_collision_set;
2402 friction_collision_set.build(
2403 collision_mesh, displaced_surface, collision_set,
2404 barrier_potential, friction_coefficient);
2405
2406 ipc::FrictionPotential friction_potential(epsv);
2407
2408 Eigen::MatrixXd velocities;
2409 if (state.solve_data.time_integrator != nullptr)
2410 velocities = state.solve_data.time_integrator->v_prev();
2411 else
2412 velocities = sol;
2413 velocities = collision_mesh.map_displacements(utils::unflatten(velocities, collision_mesh.dim()));
2414
2415 Eigen::MatrixXd forces = -friction_potential.gradient(
2416 friction_collision_set, collision_mesh, velocities);
2417
2418 Eigen::MatrixXd forces_reshaped = utils::unflatten(forces, problem_dim);
2419
2420 assert(forces_reshaped.rows() == surface_displacements.rows());
2421 assert(forces_reshaped.cols() == surface_displacements.cols());
2422 writer.add_field("friction_forces", forces_reshaped);
2423 }
2424
2425 const bool adhesion_enabled = state.is_adhesion_enabled();
2426 const bool need_normal_adhesion = adhesion_enabled && (opts.normal_adhesion_forces || opts.export_field("normal_adhesion_forces"));
2427 const bool need_tangential_adhesion = adhesion_enabled && (opts.tangential_adhesion_forces || opts.export_field("tangential_adhesion_forces"));
2428
2429 ipc::NormalCollisions adhesion_collision_set;
2430 ipc::NormalAdhesionPotential normal_adhesion_potential(dhat_p, dhat_a, Y, 1);
2431
2432 if (need_normal_adhesion || need_tangential_adhesion)
2433 {
2434 adhesion_collision_set.build(
2435 collision_mesh, displaced_surface, dhat_a,
2436 /*dmin=*/0, ipc::create_broad_phase(state.args["solver"]["contact"]["CCD"]["broad_phase"]).get());
2437 }
2438
2439 if (need_normal_adhesion)
2440 {
2441 Eigen::MatrixXd forces = -1 * normal_adhesion_potential.gradient(adhesion_collision_set, collision_mesh, displaced_surface);
2442
2443 Eigen::MatrixXd forces_reshaped = utils::unflatten(forces, problem_dim);
2444
2445 assert(forces_reshaped.rows() == surface_displacements.rows());
2446 assert(forces_reshaped.cols() == surface_displacements.cols());
2447 writer.add_field("normal_adhesion_forces", forces_reshaped);
2448 }
2449
2450 if (need_tangential_adhesion)
2451 {
2452 ipc::TangentialCollisions tangential_collision_set;
2453 tangential_collision_set.build(
2454 collision_mesh, displaced_surface, adhesion_collision_set,
2455 normal_adhesion_potential, tangential_adhesion_coefficient);
2456
2457 ipc::TangentialAdhesionPotential tangential_adhesion_potential(epsa);
2458
2459 Eigen::MatrixXd velocities;
2460 if (state.solve_data.time_integrator != nullptr)
2461 velocities = state.solve_data.time_integrator->v_prev();
2462 else
2463 velocities = sol;
2464 velocities = collision_mesh.map_displacements(utils::unflatten(velocities, collision_mesh.dim()));
2465
2466 Eigen::MatrixXd forces = -tangential_adhesion_potential.gradient(
2467 tangential_collision_set, collision_mesh, velocities);
2468
2469 Eigen::MatrixXd forces_reshaped = utils::unflatten(forces, problem_dim);
2470
2471 assert(forces_reshaped.rows() == surface_displacements.rows());
2472 assert(forces_reshaped.cols() == surface_displacements.cols());
2473 writer.add_field("tangential_adhesion_forces", forces_reshaped);
2474 }
2475
2476 assert(collision_mesh.rest_positions().rows() == surface_displacements.rows());
2477 assert(collision_mesh.rest_positions().cols() == surface_displacements.cols());
2478
2479 // Write the solution last so it is the default for warp-by-vector
2480 writer.add_field("solution", surface_displacements);
2481
2482 writer.write_mesh(
2483 export_surface.substr(0, export_surface.length() - 4) + "_contact.vtu",
2484 collision_mesh.rest_positions(),
2485 problem_dim == 3 ? collision_mesh.faces() : collision_mesh.edges(),
2486 problem_dim == 3 ? CellType::Triangle : CellType::Line);
2487 }
2488
2489 void OutGeometryData::save_wire(
2490 const std::string &name,
2491 const State &state,
2492 const Eigen::MatrixXd &sol,
2493 const double t,
2494 const ExportOptions &opts) const
2495 {
2496 const std::vector<basis::ElementBases> &gbases = state.geom_bases();
2497 const mesh::Mesh &mesh = *state.mesh;
2498 const assembler::Problem &problem = *state.problem;
2499
2500 const auto &sampler = ref_element_sampler;
2501
2502 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
2503 Eigen::MatrixXd vis_pts_poly;
2504
2505 const auto &current_bases = gbases;
2506 int seg_total_size = 0;
2507 int pts_total_size = 0;
2508 int faces_total_size = 0;
2509
2510 for (size_t i = 0; i < current_bases.size(); ++i)
2511 {
2512 const auto &bs = current_bases[i];
2513
2514 if (mesh.is_simplex(i))
2515 {
2516 pts_total_size += sampler.simplex_points().rows();
2517 seg_total_size += sampler.simplex_edges().rows();
2518 faces_total_size += sampler.simplex_faces().rows();
2519 }
2520 else if (mesh.is_cube(i))
2521 {
2522 pts_total_size += sampler.cube_points().rows();
2523 seg_total_size += sampler.cube_edges().rows();
2524 faces_total_size += sampler.cube_faces().rows();
2525 }
2526 else if (mesh.is_prism(i))
2527 {
2528 pts_total_size += sampler.prism_points().rows();
2529 seg_total_size += sampler.prism_edges().rows();
2530 faces_total_size += sampler.prism_faces().rows();
2531 }
2532 else if (mesh.is_pyramid(i))
2533 {
2534 pts_total_size += sampler.pyramid_points().rows();
2535 seg_total_size += sampler.pyramid_edges().rows();
2536 faces_total_size += sampler.pyramid_faces().rows();
2537 }
2538 else
2539 {
2540 if (mesh.is_volume())
2541 sampler.sample_polyhedron(state.polys_3d.at(i).first, state.polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
2542 else
2543 sampler.sample_polygon(state.polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
2544
2545 pts_total_size += vis_pts_poly.rows();
2546 seg_total_size += vis_edges_poly.rows();
2547 faces_total_size += vis_faces_poly.rows();
2548 }
2549 }
2550
2551 Eigen::MatrixXd points(pts_total_size, mesh.dimension());
2552 Eigen::MatrixXi edges(seg_total_size, 2);
2553 Eigen::MatrixXi faces(faces_total_size, 3);
2554 points.setZero();
2555
2556 Eigen::MatrixXd mapped, tmp;
2557 int seg_index = 0, pts_index = 0, face_index = 0;
2558 for (size_t i = 0; i < current_bases.size(); ++i)
2559 {
2560 const auto &bs = current_bases[i];
2561
2562 if (mesh.is_simplex(i))
2563 {
2564 bs.eval_geom_mapping(sampler.simplex_points(), mapped);
2565 edges.block(seg_index, 0, sampler.simplex_edges().rows(), edges.cols()) = sampler.simplex_edges().array() + pts_index;
2566 seg_index += sampler.simplex_edges().rows();
2567
2568 faces.block(face_index, 0, sampler.simplex_faces().rows(), 3) = sampler.simplex_faces().array() + pts_index;
2569 face_index += sampler.simplex_faces().rows();
2570
2571 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2572 pts_index += mapped.rows();
2573 }
2574 else if (mesh.is_cube(i))
2575 {
2576 bs.eval_geom_mapping(sampler.cube_points(), mapped);
2577 edges.block(seg_index, 0, sampler.cube_edges().rows(), edges.cols()) = sampler.cube_edges().array() + pts_index;
2578 seg_index += sampler.cube_edges().rows();
2579
2580 faces.block(face_index, 0, sampler.cube_faces().rows(), 3) = sampler.cube_faces().array() + pts_index;
2581 face_index += sampler.cube_faces().rows();
2582
2583 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2584 pts_index += mapped.rows();
2585 }
2586 else if (mesh.is_prism(i))
2587 {
2588 bs.eval_geom_mapping(sampler.prism_points(), mapped);
2589 edges.block(seg_index, 0, sampler.prism_edges().rows(), edges.cols()) = sampler.prism_edges().array() + pts_index;
2590 seg_index += sampler.prism_edges().rows();
2591
2592 faces.block(face_index, 0, sampler.prism_faces().rows(), 3) = sampler.prism_faces().array() + pts_index;
2593 face_index += sampler.prism_faces().rows();
2594
2595 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2596 pts_index += mapped.rows();
2597 }
2598 else if (mesh.is_pyramid(i))
2599 {
2600 bs.eval_geom_mapping(sampler.pyramid_points(), mapped);
2601 edges.block(seg_index, 0, sampler.pyramid_edges().rows(), edges.cols()) = sampler.pyramid_edges().array() + pts_index;
2602 seg_index += sampler.pyramid_edges().rows();
2603
2604 faces.block(face_index, 0, sampler.pyramid_faces().rows(), 3) = sampler.pyramid_faces().array() + pts_index;
2605 face_index += sampler.pyramid_faces().rows();
2606
2607 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2608 pts_index += mapped.rows();
2609 }
2610 else
2611 {
2612 if (mesh.is_volume())
2613 sampler.sample_polyhedron(state.polys_3d.at(i).first, state.polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
2614 else
2615 sampler.sample_polygon(state.polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
2616
2617 edges.block(seg_index, 0, vis_edges_poly.rows(), edges.cols()) = vis_edges_poly.array() + pts_index;
2618 seg_index += vis_edges_poly.rows();
2619
2620 faces.block(face_index, 0, vis_faces_poly.rows(), 3) = vis_faces_poly.array() + pts_index;
2621 face_index += vis_faces_poly.rows();
2622
2623 points.block(pts_index, 0, vis_pts_poly.rows(), points.cols()) = vis_pts_poly;
2624 pts_index += vis_pts_poly.rows();
2625 }
2626 }
2627
2628 assert(pts_index == points.rows());
2629 assert(face_index == faces.rows());
2630
2631 if (mesh.is_volume())
2632 {
2633 // reverse all faces
2634 for (long i = 0; i < faces.rows(); ++i)
2635 {
2636 const int v0 = faces(i, 0);
2637 const int v1 = faces(i, 1);
2638 const int v2 = faces(i, 2);
2639
2640 int tmpc = faces(i, 2);
2641 faces(i, 2) = faces(i, 1);
2642 faces(i, 1) = tmpc;
2643 }
2644 }
2645 else
2646 {
2647 Eigen::Matrix2d mmat;
2648 for (long i = 0; i < faces.rows(); ++i)
2649 {
2650 const int v0 = faces(i, 0);
2651 const int v1 = faces(i, 1);
2652 const int v2 = faces(i, 2);
2653
2654 mmat.row(0) = points.row(v2) - points.row(v0);
2655 mmat.row(1) = points.row(v1) - points.row(v0);
2656
2657 if (mmat.determinant() > 0)
2658 {
2659 int tmpc = faces(i, 2);
2660 faces(i, 2) = faces(i, 1);
2661 faces(i, 1) = tmpc;
2662 }
2663 }
2664 }
2665
2666 Eigen::MatrixXd fun;
2667 polyfem::io::Evaluator::interpolate_function(
2668 mesh, problem.is_scalar(), state.bases, state.disc_orders, state.disc_ordersq,
2669 state.polys, state.polys_3d, ref_element_sampler,
2670 pts_index, sol, fun, /*use_sampler*/ true, false);
2671
2672 Eigen::MatrixXd exact_fun, err;
2673
2674 if (problem.has_exact_sol())
2675 {
2676 problem.exact(points, t, exact_fun);
2677 err = (fun - exact_fun).eval().rowwise().norm();
2678 }
2679
2680 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2681 if (opts.use_hdf5)
2682 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2683 else
2684 tmpw = std::make_shared<paraviewo::VTUWriter>();
2685 paraviewo::ParaviewWriter &writer = *tmpw;
2686
2687 if (problem.has_exact_sol())
2688 {
2689 if (opts.export_field("exact"))
2690 writer.add_field("exact", exact_fun);
2691 if (opts.export_field("error"))
2692 writer.add_field("error", err);
2693 }
2694
2695 if (fun.cols() != 1)
2696 {
2697 std::vector<assembler::Assembler::NamedMatrix> scalar_val;
2698 polyfem::io::Evaluator::compute_scalar_value(
2699 mesh, problem.is_scalar(), state.bases, gbases,
2700 state.disc_orders, state.disc_ordersq, state.polys, state.polys_3d,
2701 *state.assembler,
2702 ref_element_sampler, pts_index, sol, t, scalar_val, /*use_sampler*/ true, false);
2703 for (const auto &v : scalar_val)
2704 {
2705 if (opts.export_field(v.first))
2706 writer.add_field(v.first, v.second);
2707 }
2708 }
2709 // Write the solution last so it is the default for warp-by-vector
2710 writer.add_field("solution", fun);
2711
2712 writer.write_mesh(name, points, edges, CellType::Line);
2713 }
2714
2715 void OutGeometryData::save_points(
2716 const std::string &path,
2717 const State &state,
2718 const Eigen::MatrixXd &sol,
2719 const ExportOptions &opts) const
2720 {
2721 const auto &dirichlet_nodes = state.dirichlet_nodes;
2722 const auto &dirichlet_nodes_position = state.dirichlet_nodes_position;
2723 const mesh::Mesh &mesh = *state.mesh;
2724 const assembler::Problem &problem = *state.problem;
2725
2726 int actual_dim = 1;
2727 if (!problem.is_scalar())
2728 actual_dim = mesh.dimension();
2729
2730 Eigen::MatrixXd fun(dirichlet_nodes_position.size(), actual_dim);
2731 Eigen::MatrixXd b_sidesets(dirichlet_nodes_position.size(), 1);
2732 b_sidesets.setZero();
2733 Eigen::MatrixXd points(dirichlet_nodes_position.size(), mesh.dimension());
2734 std::vector<CellElement> cells(dirichlet_nodes_position.size());
2735
2736 for (int i = 0; i < dirichlet_nodes_position.size(); ++i)
2737 {
2738 const int n_id = dirichlet_nodes[i];
2739 const auto s_id = mesh.get_node_id(n_id);
2740 if (s_id > 0)
2741 {
2742 b_sidesets(i) = s_id;
2743 }
2744
2745 for (int j = 0; j < actual_dim; ++j)
2746 {
2747 fun(i, j) = sol(n_id * actual_dim + j);
2748 }
2749
2750 points.row(i) = dirichlet_nodes_position[i];
2751 cells[i].vertices.push_back(i);
2752 cells[i].ctype = CellType::Vertex;
2753 }
2754
2755 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2756 if (opts.use_hdf5)
2757 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2758 else
2759 tmpw = std::make_shared<paraviewo::VTUWriter>();
2760 paraviewo::ParaviewWriter &writer = *tmpw;
2761
2762 if (opts.export_field("sidesets"))
2763 writer.add_field("sidesets", b_sidesets);
2764 // Write the solution last so it is the default for warp-by-vector
2765 writer.add_field("solution", fun);
2766 writer.write_mesh(path, points, cells);
2767 }
2768
2769 void OutGeometryData::save_pvd(
2770 const std::string &name,
2771 const std::function<std::string(int)> &vtu_names,
2772 int time_steps, double t0, double dt, int skip_frame) const
2773 {
2774 paraviewo::PVDWriter::save_pvd(name, vtu_names, time_steps, t0, dt, skip_frame);
2775 }
2776
2777 void OutGeometryData::init_sampler(const polyfem::mesh::Mesh &mesh, const double vismesh_rel_area)
2778 {
2779 ref_element_sampler.init(mesh.is_volume(), mesh.n_elements(), vismesh_rel_area);
2780 }
2781
2782 void OutGeometryData::build_grid(const polyfem::mesh::Mesh &mesh, const double spacing)
2783 {
2784 if (spacing <= 0)
2785 return;
2786
2787 RowVectorNd min, max;
2788 mesh.bounding_box(min, max);
2789 const RowVectorNd delta = max - min;
2790 const int nx = delta[0] / spacing + 1;
2791 const int ny = delta[1] / spacing + 1;
2792 const int nz = delta.cols() >= 3 ? (delta[2] / spacing + 1) : 1;
2793 const int n = nx * ny * nz;
2794
2795 grid_points.resize(n, delta.cols());
2796 int index = 0;
2797 for (int i = 0; i < nx; ++i)
2798 {
2799 const double x = (delta[0] / (nx - 1)) * i + min[0];
2800
2801 for (int j = 0; j < ny; ++j)
2802 {
2803 const double y = (delta[1] / (ny - 1)) * j + min[1];
2804
2805 if (delta.cols() <= 2)
2806 {
2807 grid_points.row(index++) << x, y;
2808 }
2809 else
2810 {
2811 for (int k = 0; k < nz; ++k)
2812 {
2813 const double z = (delta[2] / (nz - 1)) * k + min[2];
2814 grid_points.row(index++) << x, y, z;
2815 }
2816 }
2817 }
2818 }
2819
2820 assert(index == n);
2821
2822 std::vector<std::array<Eigen::Vector3d, 2>> boxes;
2823 mesh.elements_boxes(boxes);
2824
2825 SimpleBVH::BVH bvh;
2826 bvh.init(boxes);
2827
2828 const double eps = 1e-6;
2829
2830 grid_points_to_elements.resize(grid_points.rows(), 1);
2831 grid_points_to_elements.setConstant(-1);
2832
2833 grid_points_bc.resize(grid_points.rows(), mesh.is_volume() ? 4 : 3);
2834
2835 for (int i = 0; i < grid_points.rows(); ++i)
2836 {
2837 const Eigen::Vector3d min(
2838 grid_points(i, 0) - eps,
2839 grid_points(i, 1) - eps,
2840 (mesh.is_volume() ? grid_points(i, 2) : 0) - eps);
2841
2842 const Eigen::Vector3d max(
2843 grid_points(i, 0) + eps,
2844 grid_points(i, 1) + eps,
2845 (mesh.is_volume() ? grid_points(i, 2) : 0) + eps);
2846
2847 std::vector<unsigned int> candidates;
2848
2849 bvh.intersect_box(min, max, candidates);
2850
2851 for (const auto cand : candidates)
2852 {
2853 if (!mesh.is_simplex(cand))
2854 {
2855 logger().warn("Element {} is not simplex, skipping", cand);
2856 continue;
2857 }
2858
2859 Eigen::MatrixXd coords;
2860 mesh.barycentric_coords(grid_points.row(i), cand, coords);
2861
2862 for (int d = 0; d < coords.size(); ++d)
2863 {
2864 if (fabs(coords(d)) < 1e-8)
2865 coords(d) = 0;
2866 else if (fabs(coords(d) - 1) < 1e-8)
2867 coords(d) = 1;
2868 }
2869
2870 if (coords.array().minCoeff() >= 0 && coords.array().maxCoeff() <= 1)
2871 {
2872 grid_points_to_elements(i) = cand;
2873 grid_points_bc.row(i) = coords;
2874 break;
2875 }
2876 }
2877 }
2878 }
2879
2880 EnergyCSVWriter::EnergyCSVWriter(const std::string &path, const solver::SolveData &solve_data)
2881 : file(path), solve_data(solve_data)
2882 {
2883 file << "i,";
2884 for (const auto &[name, _] : solve_data.named_forms())
2885 {
2886 file << name << ",";
2887 }
2888 file << "total_energy" << std::endl;
2889 }
2890
2891 EnergyCSVWriter::~EnergyCSVWriter()
2892 {
2893 file.close();
2894 }
2895
2896 void EnergyCSVWriter::write(const int i, const Eigen::MatrixXd &sol)
2897 {
2898 const double s = solve_data.time_integrator
2899 ? solve_data.time_integrator->acceleration_scaling()
2900 : 1;
2901 file << i << ",";
2902 for (const auto &[_, form] : solve_data.named_forms())
2903 {
2904 // Divide by acceleration scaling to get the energy (units of J)
2905 file << ((form && form->enabled()) ? form->value(sol) : 0) / s << ",";
2906 }
2907 file << solve_data.nl_problem->value(sol) / s << "\n";
2908 file.flush();
2909 }
2910
2911 RuntimeStatsCSVWriter::RuntimeStatsCSVWriter(const std::string &path, const State &state, const double t0, const double dt)
2912 : file(path), state(state), t0(t0), dt(dt)
2913 {
2914 file << "step,time,forward,remeshing,global_relaxation,peak_mem,#V,#T" << std::endl;
2915 }
2916
2921
2922 void RuntimeStatsCSVWriter::write(const int t, const double forward, const double remeshing, const double global_relaxation, const Eigen::MatrixXd &sol)
2923 {
2924 total_forward_solve_time += forward;
2925 total_remeshing_time += remeshing;
2926 total_global_relaxation_time += global_relaxation;
2927
2928 // logger().debug(
2929 // "Forward (cur, avg, total): {} s, {} s, {} s",
2930 // forward, total_forward_solve_time / t, total_forward_solve_time);
2931 // logger().debug(
2932 // "Remeshing (cur, avg, total): {} s, {} s, {} s",
2933 // remeshing, total_remeshing_time / t, total_remeshing_time);
2934 // logger().debug(
2935 // "Global relaxation (cur, avg, total): {} s, {} s, {} s",
2936 // global_relaxation, total_global_relaxation_time / t, total_global_relaxation_time);
2937
2938 const double peak_mem = getPeakRSS() / double(1 << 30);
2939 // logger().debug("Peak mem: {} GiB", peak_mem);
2940
2941 file << fmt::format(
2942 "{},{},{},{},{},{},{},{}\n",
2943 t, t0 + dt * t, forward, remeshing, global_relaxation, peak_mem,
2944 state.n_bases, state.mesh->n_elements());
2945 file.flush();
2946 }
2947
2948} // namespace polyfem::legacy::io
vals
ElementAssemblyValues vals
Definition Assembler.cpp:22
BodyForm.hpp
faces
std::vector< Eigen::VectorXi > faces
Definition PressureAssembler.cpp:52
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::Assembler::parameters
virtual std::map< std::string, ParamFunc > parameters() const =0
polyfem::assembler::Assembler::compute_tensor_value
virtual void compute_tensor_value(const OutputData &data, std::vector< NamedMatrix > &result) const
Definition Assembler.hpp:133
polyfem::assembler::AssemblyValues
stores per local bases evaluations
Definition AssemblyValues.hpp:13
polyfem::assembler::AssemblyValues::global
std::vector< basis::Local2Global > global
Definition AssemblyValues.hpp:19
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::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:51
polyfem::assembler::Problem::is_time_dependent
virtual bool is_time_dependent() const
Definition Problem.hpp:56
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::Evaluator::interpolate_at_local_vals
static void interpolate_at_local_vals(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const int el_index, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &fun, Eigen::MatrixXd &result, Eigen::MatrixXd &result_grad)
interpolate solution and gradient at element (calls interpolate_at_local_vals with sol)
Definition Evaluator.cpp:660
polyfem::io::Evaluator::average_grad_based_function
static void average_grad_based_function(const mesh::Mesh &mesh, const bool is_problem_scalar, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXi &disc_orders, const Eigen::VectorXi &disc_ordersq, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const assembler::Assembler &assembler, const utils::RefElementSampler &sampler, const double t, const int n_points, const Eigen::MatrixXd &fun, std::vector< assembler::Assembler::NamedMatrix > &result_scalar, std::vector< assembler::Assembler::NamedMatrix > &result_tensor, const bool use_sampler, const bool boundary_only)
calls compute_scalar_value (i.e von mises for elasticity and norm of velocity for fluid) and compute_...
Definition Evaluator.cpp:161
polyfem::io::Evaluator::compute_scalar_value
static void compute_scalar_value(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXi &disc_orders, const Eigen::VectorXi &disc_ordersq, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const assembler::Assembler &assembler, const utils::RefElementSampler &sampler, const int n_points, const Eigen::MatrixXd &fun, const double t, std::vector< assembler::Assembler::NamedMatrix > &result, const bool use_sampler, const bool boundary_only)
computes scalar quantity of funtion (ie von mises for elasticity and norm of velocity for fluid)
Definition Evaluator.cpp:856
polyfem::io::Evaluator::compute_tensor_value
static void compute_tensor_value(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXi &disc_orders, const Eigen::VectorXi &disc_ordersq, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const assembler::Assembler &assembler, const utils::RefElementSampler &sampler, const int n_points, const Eigen::MatrixXd &fun, const double t, std::vector< assembler::Assembler::NamedMatrix > &result, const bool use_sampler, const bool boundary_only)
compute tensor quantity (ie stress tensor or velocity)
Definition Evaluator.cpp:968
polyfem::io::Evaluator::mark_flipped_cells
static void mark_flipped_cells(const mesh::Mesh &mesh, const std::vector< basis::ElementBases > &gbasis, const std::vector< basis::ElementBases > &basis, const Eigen::VectorXi &disc_orders, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const utils::RefElementSampler &sampler, const int n_points, const Eigen::MatrixXd &fun, Eigen::Vector< bool, -1 > &result, const bool use_sampler, const bool boundary_only)
Definition Evaluator.cpp:477
polyfem::io::Evaluator::interpolate_function
static void interpolate_function(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const Eigen::VectorXi &disc_orders, const Eigen::VectorXi &disc_ordersq, const std::map< int, Eigen::MatrixXd > &polys, const std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &polys_3d, const utils::RefElementSampler &sampler, const int n_points, const Eigen::MatrixXd &fun, Eigen::MatrixXd &result, const bool use_sampler, const bool boundary_only)
interpolate the function fun.
Definition Evaluator.cpp:454
polyfem::io::Evaluator::compute_stress_at_quadrature_points
static void compute_stress_at_quadrature_points(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXi &disc_orders, const Eigen::VectorXi &disc_ordersq, const assembler::Assembler &assembler, const Eigen::MatrixXd &fun, const double t, Eigen::MatrixXd &result, Eigen::VectorXd &von_mises)
compute von mises stress at quadrature points for the function fun, also compute the interpolated fun...
Definition Evaluator.cpp:342
polyfem::legacy::State
main class that contains the polyfem solver and all its state
Definition State.hpp:114
polyfem::legacy::State::problem
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
Definition State.hpp:204
polyfem::legacy::State::geom_bases
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:264
polyfem::legacy::State::collision_mesh
ipc::CollisionMesh collision_mesh
IPC collision mesh.
Definition State.hpp:672
polyfem::legacy::State::dirichlet_nodes_position
std::vector< RowVectorNd > dirichlet_nodes_position
Definition State.hpp:572
polyfem::legacy::State::dirichlet_nodes
std::vector< int > dirichlet_nodes
per node dirichlet
Definition State.hpp:571
polyfem::legacy::State::mass_matrix_assembler
std::shared_ptr< assembler::Mass > mass_matrix_assembler
Definition State.hpp:192
polyfem::legacy::State::mesh
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:594
polyfem::legacy::State::rhs
Eigen::MatrixXd rhs
System right-hand side.
Definition State.hpp:248
polyfem::legacy::State::args
json args
main input arguments containing all defaults
Definition State.hpp:136
polyfem::legacy::State::disc_ordersq
Eigen::VectorXi disc_ordersq
Definition State.hpp:226
polyfem::legacy::State::polys
std::map< int, Eigen::MatrixXd > polys
polygons, used since poly have no geom mapping
Definition State.hpp:221
polyfem::legacy::State::n_bases
int n_bases
number of bases
Definition State.hpp:214
polyfem::legacy::State::pressure_bases
std::vector< basis::ElementBases > pressure_bases
FE pressure bases for mixed elements, the size is #elements.
Definition State.hpp:209
polyfem::legacy::State::obstacle
mesh::Obstacle obstacle
Obstacles used in collisions.
Definition State.hpp:596
polyfem::legacy::State::disc_orders
Eigen::VectorXi disc_orders
vector of discretization orders, used when not all elements have the same degree, one per element
Definition State.hpp:226
polyfem::legacy::State::n_boundary_samples
QuadratureOrders n_boundary_samples() const
quadrature used for projecting boundary conditions
Definition State.hpp:304
polyfem::legacy::State::assembler
std::shared_ptr< assembler::Assembler > assembler
assemblers
Definition State.hpp:190
polyfem::legacy::State::bases
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:207
polyfem::legacy::State::polys_3d
std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > polys_3d
polyhedra, used since poly have no geom mapping
Definition State.hpp:223
polyfem::legacy::State::is_adhesion_enabled
bool is_adhesion_enabled() const
does the simulation have adhesion
Definition State.hpp:716
polyfem::legacy::State::total_local_boundary
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
Definition State.hpp:559
polyfem::legacy::State::solve_data
solver::SolveData solve_data
timedependent stuff cached
Definition State.hpp:382
polyfem::legacy::State::in_node_to_node
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
Definition State.hpp:578
polyfem::legacy::State::mixed_assembler
std::shared_ptr< assembler::MixedAssembler > mixed_assembler
Definition State.hpp:195
polyfem::legacy::io::EnergyCSVWriter::EnergyCSVWriter
EnergyCSVWriter(const std::string &path, const solver::SolveData &solve_data)
Definition OutData.cpp:2880
polyfem::legacy::io::EnergyCSVWriter::write
void write(const int i, const Eigen::MatrixXd &sol)
Definition OutData.cpp:2896
polyfem::legacy::io::EnergyCSVWriter::solve_data
const solver::SolveData & solve_data
Definition OutData.hpp:360
polyfem::legacy::io::OutGeometryData::ref_element_sampler
utils::RefElementSampler ref_element_sampler
used to sample the solution
Definition OutData.hpp:256
polyfem::legacy::io::OutGeometryData::build_high_order_vis_mesh
void build_high_order_vis_mesh(const mesh::Mesh &mesh, const Eigen::VectorXi &disc_orders, const Eigen::VectorXi &disc_ordersq, const std::vector< basis::ElementBases > &bases, Eigen::MatrixXd &points, std::vector< paraviewo::CellElement > &elements, Eigen::MatrixXi &el_id, Eigen::MatrixXd &discr) const
builds high-der visualzation mesh per element all disconnected it also retuns the mapping to element ...
polyfem::legacy::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:2769
polyfem::legacy::io::OutGeometryData::grid_points_bc
Eigen::MatrixXd grid_points_bc
grid mesh boundaries
Definition OutData.hpp:263
polyfem::legacy::io::OutGeometryData::export_data
void export_data(const State &state, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure, const bool is_time_dependent, const double tend_in, const double dt, const ExportOptions &opts, const std::string &vis_mesh_path, const std::string &nodes_path, const std::string &solution_path, const std::string &stress_path, const std::string &mises_path, const bool is_contact_enabled) const
exports everytihng, txt, vtu, etc
Definition OutData.cpp:1205
polyfem::legacy::io::OutGeometryData::save_vtu
void save_vtu(const std::string &path, const State &state, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure, const double t, const double dt, const ExportOptions &opts, const bool is_contact_enabled) const
saves the vtu file for time t
Definition OutData.cpp:1390
polyfem::legacy::io::OutGeometryData::save_wire
void save_wire(const std::string &name, const State &state, const Eigen::MatrixXd &sol, const double t, const ExportOptions &opts) const
saves the wireframe
Definition OutData.cpp:2489
polyfem::legacy::io::OutGeometryData::save_points
void save_points(const std::string &path, const State &state, const Eigen::MatrixXd &sol, const ExportOptions &opts) const
saves the nodal values
Definition OutData.cpp:2715
polyfem::legacy::io::OutGeometryData::save_volume_vector_field
void save_volume_vector_field(const State &state, const Eigen::MatrixXd &points, const ExportOptions &opts, const std::string &name, const Eigen::VectorXd &field, paraviewo::ParaviewWriter &writer) const
Definition OutData.cpp:2077
polyfem::legacy::io::OutGeometryData::save_volume
void save_volume(const std::string &path, const State &state, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure, const double t, const double dt, const ExportOptions &opts) const
saves the volume vtu file
Definition OutData.cpp:1474
polyfem::legacy::io::OutGeometryData::build_vis_boundary_mesh
void build_vis_boundary_mesh(const mesh::Mesh &mesh, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const std::vector< mesh::LocalBoundary > &total_local_boundary, const Eigen::MatrixXd &solution, const int problem_dim, 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, Eigen::MatrixXd &displaced_boundary_vis_normals) const
builds the boundary mesh for visualization
Definition OutData.cpp:551
polyfem::legacy::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) const
builds visualzation mesh, upsampled mesh used for visualization the visualization mesh is a dense mes...
Definition OutData.cpp:813
polyfem::legacy::io::OutGeometryData::save_contact_surface
void save_contact_surface(const std::string &export_surface, const State &state, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure, const double t, const double dt_in, const ExportOptions &opts, const bool is_contact_enabled) const
saves the surface vtu file for for constact quantites, eg contact or friction forces
Definition OutData.cpp:2298
polyfem::legacy::io::OutGeometryData::grid_points_to_elements
Eigen::MatrixXi grid_points_to_elements
grid mesh mapping to fe elements
Definition OutData.hpp:261
polyfem::legacy::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:151
polyfem::legacy::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:2782
polyfem::legacy::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:2777
polyfem::legacy::io::OutGeometryData::grid_points
Eigen::MatrixXd grid_points
grid mesh points to export solution sampled on a grid
Definition OutData.hpp:259
polyfem::legacy::io::OutGeometryData::save_surface
void save_surface(const std::string &export_surface, const State &state, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure, const double t, const double dt_in, const ExportOptions &opts, const bool is_contact_enabled) const
saves the surface vtu file for for surface quantites, eg traction forces
Definition OutData.cpp:2102
polyfem::legacy::io::RuntimeStatsCSVWriter::RuntimeStatsCSVWriter
RuntimeStatsCSVWriter(const std::string &path, const State &state, const double t0, const double dt)
Definition OutData.cpp:2911
polyfem::legacy::io::RuntimeStatsCSVWriter::write
void write(const int t, const double forward, const double remeshing, const double global_relaxation, const Eigen::MatrixXd &sol)
Definition OutData.cpp:2922
polyfem::mesh::LocalBoundary
Boundary primitive IDs for a single element.
Definition LocalBoundary.hpp:27
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:41
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:163
polyfem::mesh::Mesh::get_body_id
virtual int get_body_id(const int primitive) const
Get the volume selection of an element (cell in 3d, face in 2d)
Definition Mesh.hpp:514
polyfem::mesh::Mesh::edge_length
virtual double edge_length(const int gid) const
edge length
Definition Mesh.hpp:308
polyfem::mesh::Mesh::is_polytope
bool is_polytope(const int el_id) const
checks if element is polygon compatible
Definition Mesh.cpp:365
polyfem::mesh::Mesh::tri_area
virtual double tri_area(const int gid) const
area of a tri face of a tet mesh
Definition Mesh.hpp:326
polyfem::mesh::Mesh::is_simplicial
bool is_simplicial() const
checks if the mesh is simplicial
Definition Mesh.hpp:618
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:352
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:488
polyfem::mesh::Mesh::is_simplex
bool is_simplex(const int el_id) const
checks if element is simplex
Definition Mesh.cpp:422
polyfem::mesh::Mesh::quad_area
virtual double quad_area(const int gid) const
area of a quad face of an hex mesh
Definition Mesh.hpp:317
polyfem::mesh::Mesh::is_prism
bool is_prism(const int el_id) const
checks if element is a prism
Definition Mesh.cpp:427
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:576
polyfem::mesh::Mesh::dimension
int dimension() const
utily for dimension
Definition Mesh.hpp:153
polyfem::mesh::Mesh::n_faces
virtual int n_faces() const =0
number of faces
polyfem::mesh::Mesh::is_pyramid
bool is_pyramid(const int el_id) const
checks if element is a pyramid
Definition Mesh.cpp:432
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:497
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::solver::SolveData
class to store time stepping data
Definition SolveData.hpp:59
polyfem::solver::SolveData::friction_form
std::shared_ptr< solver::FrictionForm > friction_form
Definition SolveData.hpp:178
polyfem::solver::SolveData::nl_problem
std::shared_ptr< solver::NLProblem > nl_problem
Definition SolveData.hpp:170
polyfem::solver::SolveData::normal_adhesion_form
std::shared_ptr< solver::NormalAdhesionForm > normal_adhesion_form
Definition SolveData.hpp:181
polyfem::solver::SolveData::contact_form
std::shared_ptr< solver::ContactForm > contact_form
Definition SolveData.hpp:175
polyfem::solver::SolveData::named_forms
std::vector< std::pair< std::string, std::shared_ptr< solver::Form > > > named_forms() const
Definition SolveData.cpp:495
polyfem::solver::SolveData::elastic_form
std::shared_ptr< solver::ElasticForm > elastic_form
Definition SolveData.hpp:177
polyfem::solver::SolveData::time_integrator
std::shared_ptr< time_integrator::ImplicitTimeIntegrator > time_integrator
Definition SolveData.hpp:186
polyfem::solver::SolveData::tangential_adhesion_form
std::shared_ptr< solver::TangentialAdhesionForm > tangential_adhesion_form
Definition SolveData.hpp:182
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::boundary_quadrature
static bool boundary_quadrature(const mesh::LocalBoundary &local_boundary, const QuadratureOrders &order, const mesh::Mesh &mesh, const bool skip_computation, Eigen::MatrixXd &uv, Eigen::MatrixXd &points, Eigen::MatrixXd &normals, Eigen::VectorXd &weights, Eigen::VectorXi &global_primitive_ids)
Definition BoundarySampler.cpp:808
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::RefElementSampler::init
void init(const bool is_volume, const int n_elements, const double target_rel_area)
Definition RefElementSampler.cpp:217
polyfem::utils::RefElementSampler::simplex_points
const Eigen::MatrixXd & simplex_points() const
Definition RefElementSampler.hpp:51
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::legacy::io::CellType
paraviewo::CellType CellType
Definition OutData.cpp:63
polyfem::legacy::io::CellElement
paraviewo::CellElement CellElement
Definition OutData.cpp:64
polyfem::utils::unflatten
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
Definition MatrixUtils.cpp:53
polyfem::utils::append_rows_of_zeros
void append_rows_of_zeros(DstMat &dst, const size_t n_zero_rows)
Definition MatrixUtils.hpp:161
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
pyramid_duffy.max_order
max_order
Definition pyramid_duffy.py:116
par_for.hpp
prism_bases.hpp
polyfem::legacy::io::OutGeometryData::ExportOptions::export_field
bool export_field(const std::string &field) const
Definition OutData.cpp:1347
polyfem::legacy::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:1352
polyfem::legacy::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:84