44#include <paraviewo/VTMWriter.hpp>
45#include <paraviewo/PVDWriter.hpp>
47#include <SimpleBVH/BVH.hpp>
49#include <igl/write_triangle_mesh.h>
51#include <igl/facet_adjacency_matrix.h>
52#include <igl/connected_components.h>
62 void compute_traction_forces(
const State &state,
const Eigen::MatrixXd &solution,
const double t, Eigen::MatrixXd &traction_forces,
bool skip_dirichlet =
true)
65 if (!state.problem->is_scalar())
66 actual_dim = state.mesh->dimension();
70 const std::vector<basis::ElementBases> &bases = state.bases;
71 const std::vector<basis::ElementBases> &gbases = state.geom_bases();
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;
79 traction_forces.setZero(state.n_bases * actual_dim, 1);
81 for (
const auto &lb : state.total_local_boundary)
83 const int e = lb.element_id();
89 const basis::ElementBases &gbs = gbases[
e];
90 const basis::ElementBases &bs = bases[
e];
92 vals.
compute(e, state.mesh->is_volume(), points, bs, gbs);
94 for (
int n = 0; n < normals.rows(); ++n)
98 if (solution.size() > 0)
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)
105 for (
const auto &g : b.global)
107 for (
int d = 0; d < actual_dim; ++d)
109 deform_mat.row(d) += solution(
g.index * actual_dim + d) * b.grad.row(n);
117 normals.row(n) = normals.row(n) * trafo.inverse();
118 normals.row(n).normalize();
121 std::vector<assembler::Assembler::NamedMatrix> tensor_flat;
122 state.assembler->compute_tensor_value(assembler::OutputData(t, e, bs, gbs, points, solution), tensor_flat);
128 const int g_index = v.
global[0].index * actual_dim;
130 for (
int q = 0; q <
points.rows(); ++q)
133 assert(tensor_flat[0].first ==
"cauchy_stess");
134 assert(tensor_flat[0].second.row(q).size() == actual_dim * actual_dim);
136 Eigen::MatrixXd stress_tensor =
utils::unflatten(tensor_flat[0].second.row(q), actual_dim);
138 traction_forces.block(g_index, 0, actual_dim, 1) += stress_tensor * normals.row(q).transpose() * v.
val(q) * weights(q);
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)
157 displacement_map_entries.clear();
163 logger().warn(
"Skipping as the mesh has polygons");
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);
173 std::vector<std::tuple<int, int, int>> tris;
175 std::vector<bool> visited_node(n_bases,
false);
177 std::stringstream print_warning;
183 for (
int j = 0; j < lb.size(); ++j)
185 const int eid = lb.global_primitive_id(j);
186 const int lid = lb[j];
189 if (mesh.
is_cube(lb.element_id()))
191 assert(!is_simplicial);
193 std::vector<int> loc_nodes;
196 for (
long n = 0; n < nodes.size(); ++n)
198 auto &bs = b.
bases[nodes(n)];
199 const auto &glob = bs.global();
200 if (glob.size() != 1)
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);
209 if (loc_nodes.size() != 4)
211 logger().trace(
"skipping element {} since it is not Q1", eid);
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);
224 for (
int q = 0; q < 4; ++q)
226 if (!visited_node[loc_nodes[q]])
227 displacement_map_entries.emplace_back(loc_nodes[q], loc_nodes[q], 1);
229 visited_node[loc_nodes[q]] =
true;
230 displacement_map_entries.emplace_back(new_node, loc_nodes[q], 0.25);
238 logger().trace(
"skipping element {} since it is not a simplex or hex", eid);
244 std::vector<int> loc_nodes;
246 bool is_follower =
false;
249 for (
long n = 0; n < nodes.size(); ++n)
251 auto &bs = b.
bases[nodes(n)];
252 const auto &glob = bs.global();
253 if (glob.size() != 1)
264 for (
long n = 0; n < nodes.size(); ++n)
267 const std::vector<basis::Local2Global> &glob = bs.
global();
268 if (glob.size() != 1)
271 int gindex = glob.front().index;
272 node_positions.row(gindex) = glob.front().node;
273 loc_nodes.push_back(gindex);
276 if (loc_nodes.size() == 3)
278 tris.emplace_back(loc_nodes[0], loc_nodes[1], loc_nodes[2]);
280 else if (loc_nodes.size() == 6)
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]);
287 else if (loc_nodes.size() == 10)
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]);
299 else if (loc_nodes.size() == 15)
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]);
320 print_warning << loc_nodes.size() <<
" ";
326 for (
int k = 0; k < loc_nodes.size(); ++k)
328 if (!visited_node[loc_nodes[k]])
329 displacement_map_entries.emplace_back(loc_nodes[k], loc_nodes[k], 1);
331 visited_node[loc_nodes[k]] =
true;
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());
340 boundary_triangles.resize(tris.size(), 3);
341 for (
int i = 0; i < tris.size(); ++i)
343 boundary_triangles.row(i) << std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]);
346 if (boundary_triangles.rows() > 0)
348 igl::edges(boundary_triangles, boundary_edges);
353 node_positions.resize(n_bases, 2);
354 node_positions.setZero();
355 const Mesh2D &mesh2d =
dynamic_cast<const Mesh2D &
>(mesh);
357 std::vector<std::pair<int, int>> edges;
363 for (
int j = 0; j < lb.size(); ++j)
365 const int eid = lb.global_primitive_id(j);
366 const int lid = lb[j];
371 for (
long n = 0; n < nodes.size(); ++n)
374 const std::vector<basis::Local2Global> &glob = bs.
global();
375 if (glob.size() != 1)
378 int gindex = glob.front().index;
379 node_positions.row(gindex) = glob.front().node.head<2>();
382 edges.emplace_back(prev_node, gindex);
389 boundary_triangles.resize(0, 0);
390 boundary_edges.resize(edges.size(), 2);
391 for (
int i = 0; i < edges.size(); ++i)
393 boundary_edges.row(i) << edges[i].first, edges[i].second;
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
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;
423 std::vector<std::pair<int, int>> edges;
424 std::vector<std::tuple<int, int, int>> tris;
426 for (
auto it = total_local_boundary.begin(); it != total_local_boundary.end(); ++it)
428 const auto &lb = *it;
429 const auto &gbs = gbases[lb.element_id()];
430 const auto &bs = bases[lb.element_id()];
432 for (
int k = 0; k < lb.size(); ++k)
436 case BoundaryType::TRI_LINE:
440 case BoundaryType::QUAD_LINE:
444 case BoundaryType::QUAD:
448 case BoundaryType::TRI:
452 case BoundaryType::POLYGON:
456 case BoundaryType::POLYHEDRON:
459 case BoundaryType::INVALID:
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();
476 if (lb.type() == BoundaryType::QUAD)
478 const auto map = [n_samples, size](
int i,
int j) {
return j * n_samples + i + size; };
480 for (
int j = 0; j < n_samples - 1; ++j)
482 for (
int i = 0; i < n_samples - 1; ++i)
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));
489 else if (lb.type() == BoundaryType::TRI)
492 std::vector<int> mapp(n_samples * n_samples, -1);
493 for (
int j = 0; j < n_samples; ++j)
495 for (
int i = 0; i < n_samples - j; ++i)
497 mapp[j * n_samples + i] = index;
501 const auto map = [mapp, n_samples](
int i,
int j) {
502 if (j * n_samples + i >= mapp.size())
504 return mapp[j * n_samples + i];
507 for (
int j = 0; j < n_samples - 1; ++j)
509 for (
int i = 0; i < n_samples - j; ++i)
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);
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);
526 for (
int i = 0; i < vertices.back().rows() - 1; ++i)
527 edges.emplace_back(i + size, i + size + 1);
530 normals.resize(
vals.jac_it.size(), tmp_n.cols());
531 displaced_normals.resize(
vals.jac_it.size(), tmp_n.cols());
533 for (
int n = 0; n <
vals.jac_it.size(); ++n)
535 trafo =
vals.jac_it[n].inverse();
537 if (problem_dim == 2 || problem_dim == 3)
540 if (solution.size() > 0)
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);
553 normals.row(n) = tmp_n *
vals.jac_it[n];
554 normals.row(n).normalize();
556 displaced_normals.row(n) = tmp_n * trafo.inverse();
557 displaced_normals.row(n).normalize();
560 allnormals.push_back(normals);
561 displaced_allnormals.push_back(displaced_normals);
564 for (
int n = 0; n <
vals.jac_it.size(); ++n)
566 tmp_n += normals.row(n);
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);
574 Eigen::Vector3d n = e1.cross(e2);
575 Eigen::Vector3d nn = tmp_n.transpose();
579 for (
int i = tris_start; i < tris.size(); ++i)
581 tris[i] = std::tuple<int, int, int>(std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]));
586 size += vertices.back().rows();
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());
598 boundary_vis_elements.resize(tris.size(), 3);
600 boundary_vis_elements.resize(edges.size(), 2);
604 for (
const auto &v : vertices)
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++]);
614 for (
const auto &n : allnormals)
616 boundary_vis_normals.block(index, 0, n.rows(), n.cols()) = n;
621 for (
const auto &n : displaced_allnormals)
623 displaced_boundary_vis_normals.block(index, 0, n.rows(), n.cols()) = n;
630 for (
const auto &t : tris)
632 boundary_vis_elements.row(index) << std::get<0>(t), std::get<1>(t), std::get<2>(t);
638 for (
const auto &e : edges)
640 boundary_vis_elements.row(index) << e.first, e.second;
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
671 const auto ¤t_bases = gbases;
672 int tet_total_size = 0;
673 int pts_total_size = 0;
675 Eigen::MatrixXd vis_pts_poly;
676 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
678 for (
size_t i = 0; i < current_bases.size(); ++i)
680 const auto &bs = current_bases[i];
688 pts_total_size += sampler.simplex_points().rows();
692 tet_total_size += sampler.cube_volume().rows();
693 pts_total_size += sampler.cube_points().rows();
699 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
701 tet_total_size += vis_faces_poly.rows();
702 pts_total_size += vis_pts_poly.rows();
706 sampler.sample_polygon(polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
708 tet_total_size += vis_faces_poly.rows();
709 pts_total_size += vis_pts_poly.rows();
714 points.resize(pts_total_size, mesh.
dimension());
715 tets.resize(tet_total_size, mesh.
is_volume() ? 4 : 3);
717 el_id.resize(pts_total_size, 1);
718 discr.resize(pts_total_size, 1);
720 Eigen::MatrixXd mapped, tmp;
721 int tet_index = 0, pts_index = 0;
723 for (
size_t i = 0; i < current_bases.size(); ++i)
725 const auto &bs = current_bases[i];
732 bs.eval_geom_mapping(sampler.simplex_points(), mapped);
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();
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();
744 bs.eval_geom_mapping(sampler.cube_points(), mapped);
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();
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();
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);
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();
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();
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);
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();
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();
785 assert(pts_index == points.rows());
786 assert(tet_index == tets.rows());
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
810 std::vector<RowVectorNd> nodes;
811 int pts_total_size = 0;
812 elements.resize(bases.size());
813 Eigen::MatrixXd ref_pts;
815 for (
size_t i = 0; i < bases.size(); ++i)
817 const auto &bs = bases[i];
835 const int n_v =
static_cast<const mesh::Mesh2D &
>(mesh).n_face_vertices(i);
836 ref_pts.resize(n_v, 2);
840 pts_total_size += ref_pts.rows();
843 points.resize(pts_total_size, mesh.
dimension());
845 el_id.resize(pts_total_size, 1);
846 discr.resize(pts_total_size, 1);
848 Eigen::MatrixXd mapped;
851 std::string error_msg =
"";
853 for (
size_t i = 0; i < bases.size(); ++i)
855 const auto &bs = bases[i];
875 bs.eval_geom_mapping(ref_pts, mapped);
877 for (
int j = 0; j < mapped.rows(); ++j)
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);
891 const int n_nodes = elements[i].size();
892 if (disc_orders(i) >= 3)
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]);
898 if (disc_orders(i) > 4)
899 error_msg =
"Saving high-order meshes not implemented for P5+ elements!";
903 if (disc_orders(i) == 4)
905 const int n_nodes = elements[i].size();
906 std::swap(elements[i][n_nodes - 1], elements[i][n_nodes - 2]);
908 if (disc_orders(i) > 4)
909 error_msg =
"Saving high-order meshes not implemented for P5+ elements!";
912 else if (disc_orders(i) > 1)
913 error_msg =
"Saving high-order meshes not implemented for Q2+ elements!";
916 if (!error_msg.empty())
919 for (
size_t i = 0; i < bases.size(); ++i)
924 const auto &mesh2d =
static_cast<const mesh::Mesh2D &
>(mesh);
927 for (
int j = 0; j < n_v; ++j)
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);
938 assert(pts_index == points.rows());
943 const Eigen::MatrixXd &sol,
944 const Eigen::MatrixXd &pressure,
945 const bool is_time_dependent,
946 const double tend_in,
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
958 logger().error(
"Load the mesh first!");
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();
966 const Eigen::MatrixXd &rhs = state.
rhs;
971 logger().error(
"Build the bases first!");
981 logger().error(
"Solve the problem first!");
985 if (!solution_path.empty())
987 std::ofstream out(solution_path);
989 out << std::scientific;
993 Eigen::VectorXi reordering(n_bases);
994 reordering.setConstant(-1);
996 for (
int i = 0; i < in_node_to_node.size(); ++i)
998 reordering[in_node_to_node[i]] = i;
1001 Eigen::MatrixXd tmp(tmp_sol.rows(), tmp_sol.cols());
1003 for (
int i = 0; i < reordering.size(); ++i)
1005 if (reordering[i] < 0)
1008 tmp.row(reordering[i]) = tmp_sol.row(i);
1011 for (
int i = 0; i < tmp.rows(); ++i)
1013 for (
int j = 0; j < tmp.cols(); ++j)
1014 out << tmp(i, j) <<
" ";
1020 out << sol << std::endl;
1024 double tend = tend_in;
1028 if (!vis_mesh_path.empty() && !is_time_dependent)
1031 vis_mesh_path, state, sol, pressure,
1033 is_contact_enabled);
1035 if (!nodes_path.empty())
1037 Eigen::MatrixXd nodes(n_bases, mesh.
dimension());
1043 for (
size_t ii = 0; ii < b.global().size(); ++ii)
1045 const auto &lg = b.global()[ii];
1046 nodes.row(lg.index) = lg.node;
1050 std::ofstream out(nodes_path);
1055 if (!stress_path.empty())
1057 Eigen::MatrixXd result;
1058 Eigen::VectorXd mises;
1062 sol, tend, result, mises);
1063 std::ofstream out(stress_path);
1067 if (!mises_path.empty())
1069 Eigen::MatrixXd result;
1070 Eigen::VectorXd mises;
1074 sol, tend, result, mises);
1075 std::ofstream out(mises_path);
1088 fields = args[
"output"][
"paraview"][
"fields"];
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;
1099 if (args[
"output"][
"paraview"][
"options"][
"force_high_order"])
1100 use_sampler =
false;
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;
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;
1117 use_spline = args[
"space"][
"basis_type"] ==
"Spline";
1119 reorder_output = args[
"output"][
"data"][
"advanced"][
"reorder_nodes"];
1121 use_hdf5 = args[
"output"][
"paraview"][
"options"][
"use_hdf5"];
1125 const std::string &path,
1127 const Eigen::MatrixXd &sol,
1128 const Eigen::MatrixXd &pressure,
1132 const bool is_contact_enabled)
const
1136 logger().error(
"Load the mesh first!");
1140 const Eigen::MatrixXd &rhs = state.
rhs;
1144 logger().error(
"Build the bases first!");
1152 if (sol.size() <= 0)
1154 logger().error(
"Solve the problem first!");
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();
1170 is_contact_enabled);
1176 is_contact_enabled);
1189 paraviewo::VTMWriter vtm(t);
1191 vtm.add_dataset(
"Volume",
"data", path_stem + opts.
file_extension());
1193 vtm.add_dataset(
"Surface",
"data", path_stem +
"_surf" + opts.
file_extension());
1195 vtm.add_dataset(
"Contact",
"data", path_stem +
"_surf_contact" + opts.
file_extension());
1197 vtm.add_dataset(
"Wireframe",
"data", path_stem +
"_wire" + opts.
file_extension());
1199 vtm.add_dataset(
"Points",
"data", path_stem +
"_points" + opts.
file_extension());
1200 vtm.save(base_path +
".vtm");
1204 const std::string &path,
1206 const Eigen::MatrixXd &sol,
1207 const Eigen::MatrixXd &pressure,
1212 const Eigen::VectorXi &disc_orders = state.
disc_orders;
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;
1225 Eigen::MatrixXd points;
1226 Eigen::MatrixXi tets;
1227 Eigen::MatrixXi el_id;
1228 Eigen::MatrixXd discr;
1229 std::vector<std::vector<int>> elements;
1234 points, tets, el_id, discr);
1237 points, elements, el_id, discr);
1239 Eigen::MatrixXd fun, exact_fun, err, node_fun;
1244 Eigen::MatrixXd tmp, tmp_grad;
1245 Eigen::MatrixXd tmp_p, tmp_grad_p;
1247 res.setConstant(std::numeric_limits<double>::quiet_NaN());
1249 res_grad.setConstant(std::numeric_limits<double>::quiet_NaN());
1252 res_p.setConstant(std::numeric_limits<double>::quiet_NaN());
1254 res_grad_p.setConstant(std::numeric_limits<double>::quiet_NaN());
1263 Eigen::MatrixXd pt(1, bc.cols() - 1);
1264 for (
int d = 1; d < bc.cols(); ++d)
1267 mesh, problem.
is_scalar(), bases, gbases,
1268 el_id, pt, sol, tmp, tmp_grad);
1271 res_grad.row(i) = tmp_grad;
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;
1283 std::ofstream os(path +
"_sol.txt");
1286 std::ofstream osg(path +
"_grad.txt");
1289 std::ofstream osgg(path +
"_grid.txt");
1294 std::ofstream osp(path +
"_p_sol.txt");
1297 std::ofstream osgp(path +
"_p_grad.txt");
1302 Eigen::Vector<bool, -1> validity;
1315 Eigen::MatrixXd tmp = Eigen::VectorXd::LinSpaced(sol.size(), 0, sol.size() - 1);
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();
1335 problem.
exact(points, t, exact_fun);
1336 err = (fun - exact_fun).eval().rowwise().norm();
1340 exact_fun.conservativeResize(exact_fun.rows() + obstacle.
n_vertices(), exact_fun.cols());
1344 err.conservativeResize(err.rows() + obstacle.
n_vertices(), 1);
1345 err.bottomRows(obstacle.
n_vertices()).setZero();
1349 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
1351 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
1353 tmpw = std::make_shared<paraviewo::VTUWriter>();
1354 paraviewo::ParaviewWriter &writer = *tmpw;
1357 writer.add_field(
"validity", validity.cast<
double>());
1360 writer.add_field(
"nodes", node_fun);
1364 bool is_time_integrator_valid = time_integrator !=
nullptr;
1368 const Eigen::VectorXd velocity =
1369 is_time_integrator_valid ? (time_integrator->v_prev()) : Eigen::VectorXd::Zero(sol.size());
1375 const Eigen::VectorXd acceleration =
1376 is_time_integrator_valid ? (time_integrator->a_prev()) : Eigen::VectorXd::Zero(sol.size());
1390 if (form ==
nullptr)
1393 Eigen::VectorXd force;
1394 if (form->enabled())
1396 form->first_derivative(sol, force);
1401 force.setZero(sol.size());
1411 Eigen::MatrixXd interp_p;
1419 interp_p.conservativeResize(interp_p.size() + obstacle.
n_vertices(), 1);
1420 interp_p.bottomRows(obstacle.
n_vertices()).setZero();
1423 writer.add_field(
"pressure", interp_p);
1428 discr.conservativeResize(discr.size() + obstacle.
n_vertices(), 1);
1429 discr.bottomRows(obstacle.
n_vertices()).setZero();
1433 writer.add_field(
"discr", discr);
1438 writer.add_field(
"exact", exact_fun);
1440 writer.add_field(
"error", err);
1443 if (fun.cols() != 1)
1445 std::vector<assembler::Assembler::NamedMatrix>
vals, tvals;
1447 mesh, problem.
is_scalar(), bases, gbases,
1452 for (
auto &[_, v] :
vals)
1457 for (
const auto &[name, v] :
vals)
1460 writer.add_field(name, v);
1471 for (
auto &[_, v] : tvals)
1474 for (
const auto &[name, v] : tvals)
1477 assert(v.cols() % stride == 0);
1482 for (
int i = 0; i < v.cols(); i += stride)
1484 const Eigen::MatrixXd tmp = v.middleCols(i, stride);
1485 assert(tmp.cols() == stride);
1487 const int ii = (i / stride) + 1;
1488 writer.add_field(fmt::format(
"{:s}_{:d}", name, ii), tmp);
1503 for (
auto &v :
vals)
1505 v.second.conservativeResize(v.second.size() + obstacle.
n_vertices(), 1);
1506 v.second.bottomRows(obstacle.
n_vertices()).setZero();
1512 for (
const auto &v :
vals)
1514 if (opts.
export_field(fmt::format(
"{:s}_avg", v.first)))
1515 writer.add_field(fmt::format(
"{:s}_avg", v.first), v.second);
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);
1530 Eigen::MatrixXd local_pts;
1531 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
1535 for (
int e = 0; e < int(bases.size()); ++e)
1543 local_pts = sampler.simplex_points();
1545 local_pts = sampler.cube_points();
1549 sampler.sample_polyhedron(polys_3d.at(e).first, polys_3d.at(e).second, local_pts, vis_faces_poly, vis_edges_poly);
1551 sampler.sample_polygon(polys.at(e), local_pts, vis_faces_poly, vis_edges_poly);
1573 const auto &mesh2d =
static_cast<const mesh::Mesh2D &
>(mesh);
1575 local_pts.resize(n_v, 2);
1577 for (
int j = 0; j < n_v; ++j)
1579 local_pts.row(j) = mesh2d.point(mesh2d.face_vertex(e, j));
1588 for (
int j = 0; j <
vals.val.rows(); ++j)
1590 for (
const auto &[p, func] : params)
1591 param_val.at(p)(index) = func(local_pts.row(j),
vals.val.row(j), t, e);
1593 rhos(index) = density(local_pts.row(j),
vals.val.row(j), t, e);
1599 assert(index == points.rows());
1603 for (
auto &[_, tmp] : param_val)
1605 tmp.conservativeResize(tmp.size() + obstacle.
n_vertices(), 1);
1606 tmp.bottomRows(obstacle.
n_vertices()).setZero();
1609 rhos.conservativeResize(rhos.size() + obstacle.
n_vertices(), 1);
1610 rhos.bottomRows(obstacle.
n_vertices()).setZero();
1612 for (
const auto &[p, tmp] : param_val)
1615 writer.add_field(p, tmp);
1618 writer.add_field(
"rho", rhos);
1624 Eigen::MatrixXd ids(points.rows(), 1);
1626 for (
int i = 0; i < points.rows(); ++i)
1633 ids.conservativeResize(ids.size() + obstacle.
n_vertices(), 1);
1634 ids.bottomRows(obstacle.
n_vertices()).setZero();
1637 writer.add_field(
"body_ids", ids);
1648 Eigen::MatrixXd traction_forces, traction_forces_fun;
1649 compute_traction_forces(state, sol, t, traction_forces,
false);
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();
1662 writer.add_field(
"traction_force", traction_forces_fun);
1669 Eigen::VectorXd potential_grad;
1670 Eigen::MatrixXd potential_grad_fun;
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();
1685 writer.add_field(
"gradient_of_elastic_potential", potential_grad_fun);
1687 catch (std::exception &)
1696 Eigen::VectorXd potential_grad;
1697 Eigen::MatrixXd potential_grad_fun;
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();
1714 writer.add_field(
"gradient_of_contact_potential", potential_grad_fun);
1717 catch (std::exception &)
1723 writer.add_field(
"solution", fun);
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();
1731 if (elements.empty())
1733 for (
int i = 0; i < tets.rows(); ++i)
1735 elements.emplace_back();
1736 for (
int j = 0; j < tets.cols(); ++j)
1737 elements.back().push_back(tets(i, j));
1743 elements.emplace_back();
1750 elements.emplace_back();
1757 elements.emplace_back();
1762 if (elements.empty())
1763 writer.write_mesh(path, points, tets);
1765 writer.write_mesh(path, points, elements,
true, disc_orders.maxCoeff() == 1);
1770 const Eigen::MatrixXd &points,
1772 const std::string &name,
1773 const Eigen::VectorXd &field,
1774 paraviewo::ParaviewWriter &writer)
const
1776 Eigen::MatrixXd inerpolated_field;
1784 inerpolated_field.conservativeResize(
1790 writer.add_field(name, inerpolated_field);
1794 const std::string &export_surface,
1796 const Eigen::MatrixXd &sol,
1797 const Eigen::MatrixXd &pressure,
1801 const bool is_contact_enabled)
const
1804 const Eigen::VectorXi &disc_orders = state.
disc_orders;
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();
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;
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);
1827 Eigen::MatrixXd fun, interp_p, discr, vect, b_sidesets;
1829 Eigen::MatrixXd lsol, lp, lgrad, lpgrad;
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());
1840 b_sidesets.resize(boundary_vis_vertices.rows(), 1);
1841 b_sidesets.setZero();
1843 for (
int i = 0; i < boundary_vis_vertices.rows(); ++i)
1845 const auto s_id = mesh.
get_boundary_id(boundary_vis_primitive_ids(i));
1848 b_sidesets(i) = s_id;
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);
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);
1865 discr(i) = disc_orders(el_index);
1866 for (
int j = 0; j < actual_dim; ++j)
1868 fun(i, j) = lsol(j);
1871 if (actual_dim == 1)
1873 assert(lgrad.size() == mesh.
dimension());
1874 for (
int j = 0; j < mesh.
dimension(); ++j)
1876 vect(i, j) = lgrad(j);
1881 assert(lgrad.size() == actual_dim * actual_dim);
1882 std::vector<assembler::Assembler::NamedMatrix> tensor_flat;
1887 assert(tensor_flat[0].first ==
"cauchy_stess");
1888 assert(tensor_flat[0].second.size() == actual_dim * actual_dim);
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;
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));
1902 area = mesh.
edge_length(boundary_vis_primitive_ids(i));
1904 vect.row(i) *= area;
1908 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
1910 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
1912 tmpw = std::make_shared<paraviewo::VTUWriter>();
1913 paraviewo::ParaviewWriter &writer = *tmpw;
1916 writer.add_field(
"normals", boundary_vis_normals);
1918 writer.add_field(
"displaced_normals", displaced_boundary_vis_normals);
1920 writer.add_field(
"pressure", interp_p);
1922 writer.add_field(
"discr", discr);
1924 writer.add_field(
"sidesets", b_sidesets);
1926 if (actual_dim == 1 && opts.
export_field(
"solution_grad"))
1927 writer.add_field(
"solution_grad", vect);
1930 writer.add_field(
"traction_force", vect);
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);
1942 for (
int i = 0; i < boundary_vis_vertices.rows(); ++i)
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));
1949 rhos(i) = density(boundary_vis_local_vertices.row(i), boundary_vis_vertices.row(i), t, boundary_vis_elements_ids(i));
1952 for (
const auto &[p, tmp] : param_val)
1955 writer.add_field(p, tmp);
1958 writer.add_field(
"rho", rhos);
1964 Eigen::MatrixXd ids(boundary_vis_vertices.rows(), 1);
1966 for (
int i = 0; i < boundary_vis_vertices.rows(); ++i)
1968 ids(i) = mesh.
get_body_id(boundary_vis_elements_ids(i));
1971 writer.add_field(
"body_ids", ids);
1975 writer.add_field(
"solution", fun);
1976 writer.write_mesh(export_surface, boundary_vis_vertices, boundary_vis_elements);
1980 const std::string &export_surface,
1982 const Eigen::MatrixXd &sol,
1983 const Eigen::MatrixXd &pressure,
1987 const bool is_contact_enabled)
const
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"];
2004 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2006 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2008 tmpw = std::make_shared<paraviewo::VTUWriter>();
2009 paraviewo::ParaviewWriter &writer = *tmpw;
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);
2015 const Eigen::MatrixXd displaced_surface = collision_mesh.displace_vertices(full_displacements);
2017 ipc::NormalCollisions collision_set;
2019 if (state.
args[
"contact"][
"use_convergent_formulation"])
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"]);
2025 collision_set.build(
2026 collision_mesh, displaced_surface, dhat,
2027 0, ipc::build_broad_phase(state.
args[
"solver"][
"contact"][
"CCD"][
"broad_phase"]));
2029 ipc::BarrierPotential barrier_potential(dhat);
2030 if (state.
args[
"contact"][
"use_convergent_formulation"])
2032 barrier_potential.set_use_physical_barrier(state.
args[
"contact"][
"use_physical_barrier"]);
2035 const double barrier_stiffness = contact_form !=
nullptr ? contact_form->barrier_stiffness() : 1;
2039 Eigen::MatrixXd forces = -barrier_stiffness * barrier_potential.gradient(collision_set, collision_mesh, displaced_surface);
2043 assert(forces_reshaped.rows() == surface_displacements.rows());
2044 assert(forces_reshaped.cols() == surface_displacements.cols());
2045 writer.add_field(
"contact_forces", forces_reshaped);
2048 if (contact_form && state.
args[
"contact"][
"use_gcp_formulation"] && state.
args[
"contact"][
"use_adaptive_dhat"] && opts.
export_field(
"adaptive_dhat"))
2050 const auto form = std::dynamic_pointer_cast<solver::SmoothContactForm>(contact_form);
2052 const auto &set = form->collision_set();
2054 if (problem_dim == 2)
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);
2061 writer.add_cell_field(
"dhat", dhats);
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);
2070 writer.add_cell_field(
"dhat_face", fdhats);
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);
2077 writer.add_field(
"dhat_vert", vdhats);
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);
2088 ipc::FrictionPotential friction_potential(epsv);
2090 Eigen::MatrixXd velocities;
2095 velocities = collision_mesh.map_displacements(
utils::unflatten(velocities, collision_mesh.dim()));
2097 Eigen::MatrixXd forces = -friction_potential.gradient(
2098 friction_collision_set, collision_mesh, velocities);
2102 assert(forces_reshaped.rows() == surface_displacements.rows());
2103 assert(forces_reshaped.cols() == surface_displacements.cols());
2104 writer.add_field(
"friction_forces", forces_reshaped);
2107 ipc::NormalCollisions adhesion_collision_set;
2108 adhesion_collision_set.build(
2109 collision_mesh, displaced_surface, dhat_a,
2110 0, ipc::build_broad_phase(state.
args[
"solver"][
"contact"][
"CCD"][
"broad_phase"]));
2112 ipc::NormalAdhesionPotential normal_adhesion_potential(dhat_p, dhat_a, Y, 1);
2116 Eigen::MatrixXd forces = -1 * normal_adhesion_potential.gradient(adhesion_collision_set, collision_mesh, displaced_surface);
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);
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);
2132 ipc::TangentialAdhesionPotential tangential_adhesion_potential(epsa);
2134 Eigen::MatrixXd velocities;
2139 velocities = collision_mesh.map_displacements(
utils::unflatten(velocities, collision_mesh.dim()));
2141 Eigen::MatrixXd forces = -tangential_adhesion_potential.gradient(
2142 tangential_collision_set, collision_mesh, velocities);
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);
2151 assert(collision_mesh.rest_positions().rows() == surface_displacements.rows());
2152 assert(collision_mesh.rest_positions().cols() == surface_displacements.cols());
2155 writer.add_field(
"solution", surface_displacements);
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());
2164 const std::string &name,
2166 const Eigen::MatrixXd &sol,
2170 const std::vector<basis::ElementBases> &gbases = state.
geom_bases();
2176 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
2177 Eigen::MatrixXd vis_pts_poly;
2179 const auto ¤t_bases = gbases;
2180 int seg_total_size = 0;
2181 int pts_total_size = 0;
2182 int faces_total_size = 0;
2184 for (
size_t i = 0; i < current_bases.size(); ++i)
2186 const auto &bs = current_bases[i];
2191 seg_total_size += sampler.simplex_edges().rows();
2192 faces_total_size += sampler.simplex_faces().rows();
2196 pts_total_size += sampler.cube_points().rows();
2197 seg_total_size += sampler.cube_edges().rows();
2198 faces_total_size += sampler.cube_faces().rows();
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);
2205 sampler.sample_polygon(state.
polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
2207 pts_total_size += vis_pts_poly.rows();
2208 seg_total_size += vis_edges_poly.rows();
2209 faces_total_size += vis_faces_poly.rows();
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);
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)
2222 const auto &bs = current_bases[i];
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();
2230 faces.block(face_index, 0, sampler.simplex_faces().rows(), 3) = sampler.simplex_faces().array() + pts_index;
2231 face_index += sampler.simplex_faces().rows();
2233 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2234 pts_index += mapped.rows();
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();
2242 faces.block(face_index, 0, sampler.cube_faces().rows(), 3) = sampler.cube_faces().array() + pts_index;
2243 face_index += sampler.cube_faces().rows();
2245 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2246 pts_index += mapped.rows();
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);
2253 sampler.sample_polygon(state.
polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
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();
2258 faces.block(face_index, 0, vis_faces_poly.rows(), 3) = vis_faces_poly.array() + pts_index;
2259 face_index += vis_faces_poly.rows();
2261 points.block(pts_index, 0, vis_pts_poly.rows(), points.cols()) = vis_pts_poly;
2262 pts_index += vis_pts_poly.rows();
2266 assert(pts_index == points.rows());
2267 assert(face_index ==
faces.rows());
2272 for (
long i = 0; i <
faces.rows(); ++i)
2274 const int v0 =
faces(i, 0);
2275 const int v1 =
faces(i, 1);
2276 const int v2 =
faces(i, 2);
2278 int tmpc =
faces(i, 2);
2285 Eigen::Matrix2d mmat;
2286 for (
long i = 0; i <
faces.rows(); ++i)
2288 const int v0 =
faces(i, 0);
2289 const int v1 =
faces(i, 1);
2290 const int v2 =
faces(i, 2);
2292 mmat.row(0) = points.row(v2) - points.row(v0);
2293 mmat.row(1) = points.row(v1) - points.row(v0);
2295 if (mmat.determinant() > 0)
2297 int tmpc =
faces(i, 2);
2304 Eigen::MatrixXd fun;
2308 pts_index, sol, fun,
true,
false);
2310 Eigen::MatrixXd exact_fun, err;
2314 problem.
exact(points, t, exact_fun);
2315 err = (fun - exact_fun).eval().rowwise().norm();
2318 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2320 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2322 tmpw = std::make_shared<paraviewo::VTUWriter>();
2323 paraviewo::ParaviewWriter &writer = *tmpw;
2328 writer.add_field(
"exact", exact_fun);
2330 writer.add_field(
"error", err);
2333 if (fun.cols() != 1)
2335 std::vector<assembler::Assembler::NamedMatrix> scalar_val;
2341 for (
const auto &v : scalar_val)
2344 writer.add_field(v.first, v.second);
2348 writer.add_field(
"solution", fun);
2350 writer.write_mesh(name, points, edges);
2354 const std::string &path,
2356 const Eigen::MatrixXd &sol,
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());
2374 for (
int i = 0; i < dirichlet_nodes_position.size(); ++i)
2376 const int n_id = dirichlet_nodes[i];
2380 b_sidesets(i) = s_id;
2383 for (
int j = 0; j < actual_dim; ++j)
2385 fun(i, j) = sol(n_id * actual_dim + j);
2388 points.row(i) = dirichlet_nodes_position[i];
2389 cells[i].push_back(i);
2392 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2394 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2396 tmpw = std::make_shared<paraviewo::VTUWriter>();
2397 paraviewo::ParaviewWriter &writer = *tmpw;
2400 writer.add_field(
"sidesets", b_sidesets);
2402 writer.add_field(
"solution", fun);
2403 writer.write_mesh(path, points, cells,
false,
false);
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
2411 paraviewo::PVDWriter::save_pvd(name, vtu_names, time_steps, t0, dt, skip_frame);
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;
2434 for (
int i = 0; i < nx; ++i)
2436 const double x = (delta[0] / (nx - 1)) * i + min[0];
2438 for (
int j = 0; j < ny; ++j)
2440 const double y = (delta[1] / (ny - 1)) * j + min[1];
2442 if (delta.cols() <= 2)
2448 for (
int k = 0; k < nz; ++k)
2450 const double z = (delta[2] / (nz - 1)) * k + min[2];
2459 std::vector<std::array<Eigen::Vector3d, 2>> boxes;
2465 const double eps = 1e-6;
2474 const Eigen::Vector3d min(
2479 const Eigen::Vector3d max(
2484 std::vector<unsigned int> candidates;
2486 bvh.intersect_box(min, max, candidates);
2488 for (
const auto cand : candidates)
2492 logger().warn(
"Element {} is not simplex, skipping", cand);
2496 Eigen::MatrixXd coords;
2499 for (
int d = 0; d < coords.size(); ++d)
2501 if (fabs(coords(d)) < 1e-8)
2503 else if (fabs(coords(d) - 1) < 1e-8)
2507 if (coords.array().minCoeff() >= 0 && coords.array().maxCoeff() <= 1)
2519 Eigen::MatrixXd samples_simplex, samples_cube, mapped, p0, p1, p;
2522 average_edge_length = 0;
2523 min_edge_length = std::numeric_limits<double>::max();
2525 if (!use_curved_mesh_size)
2529 min_edge_length = p.rowwise().norm().minCoeff();
2530 average_edge_length = p.rowwise().norm().mean();
2531 mesh_size = p.rowwise().norm().maxCoeff();
2533 logger().info(
"hmin: {}", min_edge_length);
2534 logger().info(
"hmax: {}", mesh_size);
2535 logger().info(
"havg: {}", average_edge_length);
2552 for (
size_t i = 0; i < bases_in.size(); ++i)
2561 bases_in[i].eval_geom_mapping(samples_simplex, mapped);
2566 bases_in[i].eval_geom_mapping(samples_cube, mapped);
2569 for (
int j = 0; j < n_edges; ++j)
2571 double current_edge = 0;
2572 for (
int k = 0; k < n_samples - 1; ++k)
2574 p0 = mapped.row(j * n_samples + k);
2575 p1 = mapped.row(j * n_samples + k + 1);
2578 current_edge += p.norm();
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;
2588 average_edge_length /= n;
2590 logger().info(
"hmin: {}", min_edge_length);
2591 logger().info(
"hmax: {}", mesh_size);
2592 logger().info(
"havg: {}", average_edge_length);
2606 using namespace mesh;
2608 logger().info(
"Counting flipped elements...");
2612 for (
size_t i = 0; i < gbases.size(); ++i)
2618 if (!
vals.is_geom_mapping_positive(mesh.
is_volume(), gbases[i]))
2622 static const std::vector<std::string> element_type_names{{
2624 "RegularInteriorCube",
2625 "RegularBoundaryCube",
2626 "SimpleSingularInteriorCube",
2627 "MultiSingularInteriorCube",
2628 "SimpleSingularBoundaryCube",
2630 "MultiSingularBoundaryCube",
2636 log_and_throw_error(
"element {} is flipped, type {}", i, element_type_names[
static_cast<int>(els_tag[i])]);
2651 const std::vector<polyfem::basis::ElementBases> &bases,
2652 const std::vector<polyfem::basis::ElementBases> &gbases,
2656 const Eigen::MatrixXd &sol)
2660 logger().error(
"Build the bases first!");
2663 if (sol.size() <= 0)
2665 logger().error(
"Solve the problem first!");
2675 logger().info(
"Computing errors...");
2678 const int n_el = int(bases.size());
2680 Eigen::MatrixXd v_exact, v_approx;
2681 Eigen::MatrixXd v_exact_grad(0, 0), v_approx_grad;
2691 static const int p = 8;
2696 for (
int e = 0; e < n_el; ++e)
2706 v_approx.resize(
vals.val.rows(), actual_dim);
2709 v_approx_grad.resize(
vals.val.rows(), mesh.
dimension() * actual_dim);
2710 v_approx_grad.setZero();
2712 const int n_loc_bases = int(
vals.basis_values.size());
2714 for (
int i = 0; i < n_loc_bases; ++i)
2716 const auto &
val =
vals.basis_values[i];
2718 for (
size_t ii = 0; ii <
val.global.size(); ++ii)
2720 for (
int d = 0; d < actual_dim; ++d)
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;
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();
2734 linf_err = std::max(linf_err, err.maxCoeff());
2735 grad_max_err = std::max(linf_err, err_grad.maxCoeff());
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();
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));
2786 lp_err = pow(fabs(lp_err), 1. / p);
2791 const double computing_errors_time = timer.getElapsedTime();
2792 logger().info(
" took {}s", computing_errors_time);
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);
2800 logger().info(
"-- Linf error: {}", linf_err);
2801 logger().info(
"-- grad max error: {}", grad_max_err);
2816 regular_boundary_count = 0;
2817 simple_singular_count = 0;
2818 multi_singular_count = 0;
2820 non_regular_boundary_count = 0;
2821 non_regular_count = 0;
2822 undefined_count = 0;
2823 multi_singular_boundary_count = 0;
2827 for (
size_t i = 0; i < els_tag.size(); ++i)
2833 case ElementType::SIMPLEX:
2836 case ElementType::REGULAR_INTERIOR_CUBE:
2839 case ElementType::REGULAR_BOUNDARY_CUBE:
2840 regular_boundary_count++;
2842 case ElementType::SIMPLE_SINGULAR_INTERIOR_CUBE:
2843 simple_singular_count++;
2845 case ElementType::MULTI_SINGULAR_INTERIOR_CUBE:
2846 multi_singular_count++;
2848 case ElementType::SIMPLE_SINGULAR_BOUNDARY_CUBE:
2851 case ElementType::INTERFACE_CUBE:
2852 case ElementType::MULTI_SINGULAR_BOUNDARY_CUBE:
2853 multi_singular_boundary_count++;
2855 case ElementType::BOUNDARY_POLYTOPE:
2856 non_regular_boundary_count++;
2858 case ElementType::INTERIOR_POLYTOPE:
2859 non_regular_count++;
2861 case ElementType::UNDEFINED:
2865 throw std::runtime_error(
"Unknown element type");
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);
2883 const nlohmann::json &args,
2884 const int n_bases,
const int n_pressure_bases,
2885 const Eigen::MatrixXd &sol,
2887 const Eigen::VectorXi &disc_orders,
2890 const std::string &formulation,
2891 const bool isoparametric,
2892 const int sol_at_node_id,
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;
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();
2919 j[
"mesh_size"] = mesh_size;
2920 j[
"max_angle"] = max_angle;
2922 j[
"sigma_max"] = sigma_max;
2923 j[
"sigma_min"] = sigma_min;
2924 j[
"sigma_avg"] = sigma_avg;
2926 j[
"min_edge_length"] = min_edge_length;
2927 j[
"average_edge_length"] = average_edge_length;
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;
2936 j[
"spectrum"] = {spectrum(0), spectrum(1), spectrum(2), spectrum(3)};
2937 j[
"spectrum_condest"] = std::abs(spectrum(3)) / std::abs(spectrum(0));
2950 j[
"solver_info"] = solver_info;
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;
2963 j[
"is_simplicial"] = mesh.
n_elements() == simplex_count;
2965 j[
"peak_memory"] =
getPeakRSS() / (1024 * 1024);
2969 std::vector<double> mmin(actual_dim);
2970 std::vector<double> mmax(actual_dim);
2972 for (
int d = 0; d < actual_dim; ++d)
2974 mmin[d] = std::numeric_limits<double>::max();
2975 mmax[d] = -std::numeric_limits<double>::max();
2978 for (
int i = 0; i < sol.size(); i += actual_dim)
2980 for (
int d = 0; d < actual_dim; ++d)
2982 mmin[d] = std::min(mmin[d], sol(i + d));
2983 mmax[d] = std::max(mmax[d], sol(i + d));
2987 std::vector<double> sol_at_node(actual_dim);
2989 if (sol_at_node_id >= 0)
2991 const int node_id = sol_at_node_id;
2993 for (
int d = 0; d < actual_dim; ++d)
2995 sol_at_node[d] = sol(node_id * actual_dim + d);
2999 j[
"sol_at_node"] = sol_at_node;
3000 j[
"sol_min"] = mmin;
3001 j[
"sol_max"] = mmax;
3003#if defined(POLYFEM_WITH_CPP_THREADS)
3005#elif defined(POLYFEM_WITH_TBB)
3008 j[
"num_threads"] = 1;
3011 j[
"formulation"] = formulation;
3017 : file(path), solve_data(solve_data)
3022 file << name <<
",";
3024 file <<
"total_energy" << std::endl;
3041 file << ((form && form->enabled()) ? form->value(sol) : 0) / s <<
",";
3048 : file(path), state(state), t0(t0), dt(dt)
3050 file <<
"step,time,forward,remeshing,global_relaxation,peak_mem,#V,#T" << std::endl;
3074 const double peak_mem =
getPeakRSS() / double(1 << 30);
3077 file << fmt::format(
3078 "{},{},{},{},{},{},{},{}\n",
3079 t,
t0 +
dt * t, forward, remeshing, global_relaxation, peak_mem,
ElementAssemblyValues vals
std::vector< Eigen::VectorXi > faces
main class that contains the polyfem solver and all its state
int n_bases
number of bases
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
mesh::Obstacle obstacle
Obstacles used in collisions.
std::shared_ptr< assembler::Assembler > assembler
assemblers
ipc::CollisionMesh collision_mesh
IPC collision mesh.
std::vector< basis::ElementBases > pressure_bases
FE pressure bases for mixed elements, the size is #elements.
std::shared_ptr< assembler::Mass > mass_matrix_assembler
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
std::vector< int > dirichlet_nodes
per node dirichlet
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
std::vector< RowVectorNd > dirichlet_nodes_position
std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > polys_3d
polyhedra, used since poly have no geom mapping
json args
main input arguments containing all defaults
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
solver::SolveData solve_data
timedependent stuff cached
Eigen::VectorXi disc_orders
vector of discretization orders, used when not all elements have the same degree, one per element
std::shared_ptr< assembler::MixedAssembler > mixed_assembler
std::map< int, Eigen::MatrixXd > polys
polygons, used since poly have no geom mapping
Eigen::MatrixXd rhs
System right-hand side.
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
std::vector< AssemblyValues > basis_values
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
virtual void exact_grad(const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const
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
virtual bool is_time_dependent() const
Represents one basis function and its gradient.
const std::vector< Local2Global > & global() const
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)
const solver::SolveData & solve_data
void write(const int i, const Eigen::MatrixXd &sol)
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...
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 ...
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
Eigen::MatrixXd grid_points_bc
grid mesh boundaries
Eigen::MatrixXd grid_points
grid mesh points to export solution sampled on a grid
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
void save_points(const std::string &path, const State &state, const Eigen::MatrixXd &sol, const ExportOptions &opts) const
saves the nodal values
void build_grid(const polyfem::mesh::Mesh &mesh, const double spacing)
builds the grid to export the solution
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
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
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
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
void save_wire(const std::string &name, const State &state, const Eigen::MatrixXd &sol, const double t, const ExportOptions &opts) const
saves the wireframe
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
void init_sampler(const polyfem::mesh::Mesh &mesh, const double vismesh_rel_area)
unitalize the ref element sampler
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
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
Eigen::MatrixXi grid_points_to_elements
grid mesh mapping to fe elements
utils::RefElementSampler ref_element_sampler
used to sample the solution
double loading_mesh_time
time to load the mesh
double assembling_stiffness_mat_time
time to assembly
double assigning_rhs_time
time to computing the rhs
double assembling_mass_mat_time
time to assembly mass
double building_basis_time
time to construct the basis
double solving_time
time to solve
double computing_poly_basis_time
time to build the polygonal/polyhedral bases
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
void count_flipped_elements(const polyfem::mesh::Mesh &mesh, const std::vector< polyfem::basis::ElementBases > &gbases)
counts the number of flipped elements
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
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...
void reset()
clears all stats
void compute_mesh_stats(const polyfem::mesh::Mesh &mesh)
compute stats (counts els type, mesh lenght, etc), step 1 of solve
double total_forward_solve_time
void write(const int t, const double forward, const double remeshing, const double global_relaxation, const Eigen::MatrixXd &sol)
double total_remeshing_time
RuntimeStatsCSVWriter(const std::string &path, const State &state, const double t0, const double dt)
double total_global_relaxation_time
Boundary primitive IDs for a single element.
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
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)
virtual double edge_length(const int gid) const
edge length
bool is_polytope(const int el_id) const
checks if element is polygon compatible
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
bool is_simplicial() const
checks if the mesh is simplicial
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
const Eigen::MatrixXi & orders() const
order of each element
virtual int get_boundary_id(const int primitive) const
Get the boundary selection of an element (face in 3d, edge in 2d)
bool is_simplex(const int el_id) const
checks if element is simples compatible
virtual double quad_area(const int gid) const
area of a quad face of an hex mesh
virtual bool is_volume() const =0
checks if mesh is volume
bool has_poly() const
checks if the mesh has polytopes
int dimension() const
utily for dimension
virtual int n_faces() const =0
number of faces
const std::vector< ElementType > & elements_tag() const
Returns the elements types.
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.
const Eigen::MatrixXi & get_edge_connectivity() const
const Eigen::MatrixXi & get_face_connectivity() const
const Eigen::MatrixXd & v() const
const Eigen::VectorXi & get_vertex_connectivity() const
class to store time stepping data
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,...
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.
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.
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
void log_and_throw_error(const std::string &msg)
bool tangential_adhesion_forces
std::string file_extension() const
return the extension of the output paraview files depending on use_hdf5
std::vector< std::string > fields
bool discretization_order
bool normal_adhesion_forces
ExportOptions(const json &args, const bool is_mesh_linear, const bool is_problem_scalar)
initialize the flags based on the input args
bool export_field(const std::string &field) const