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