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