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);
388 boundary_triangles.resize(0, 0);
389 boundary_edges.resize(edges.size(), 2);
390 for (
int i = 0; i < edges.size(); ++i)
392 boundary_edges.row(i) << edges[i].first, edges[i].second;
399 const std::vector<basis::ElementBases> &bases,
400 const std::vector<basis::ElementBases> &gbases,
401 const std::vector<mesh::LocalBoundary> &total_local_boundary,
402 const Eigen::MatrixXd &solution,
403 const int problem_dim,
404 Eigen::MatrixXd &boundary_vis_vertices,
405 Eigen::MatrixXd &boundary_vis_local_vertices,
406 Eigen::MatrixXi &boundary_vis_elements,
407 Eigen::MatrixXi &boundary_vis_elements_ids,
408 Eigen::MatrixXi &boundary_vis_primitive_ids,
409 Eigen::MatrixXd &boundary_vis_normals,
410 Eigen::MatrixXd &displaced_boundary_vis_normals)
const
414 std::vector<Eigen::MatrixXd> lv, vertices, allnormals, displaced_allnormals;
415 std::vector<int> el_ids, global_primitive_ids;
416 Eigen::MatrixXd uv, local_pts, tmp_n, normals, displaced_normals, trafo, deform_mat;
422 std::vector<std::pair<int, int>> edges;
423 std::vector<std::tuple<int, int, int>> tris;
425 for (
auto it = total_local_boundary.begin(); it != total_local_boundary.end(); ++it)
427 const auto &lb = *it;
428 const auto &gbs = gbases[lb.element_id()];
429 const auto &bs = bases[lb.element_id()];
431 for (
int k = 0; k < lb.size(); ++k)
435 case BoundaryType::TRI_LINE:
439 case BoundaryType::QUAD_LINE:
443 case BoundaryType::QUAD:
447 case BoundaryType::TRI:
451 case BoundaryType::POLYGON:
455 case BoundaryType::POLYHEDRON:
458 case BoundaryType::INVALID:
465 vertices.emplace_back();
466 lv.emplace_back(local_pts);
467 el_ids.push_back(lb.element_id());
468 global_primitive_ids.push_back(lb.global_primitive_id(k));
469 gbs.eval_geom_mapping(local_pts, vertices.back());
470 vals.compute(lb.element_id(), mesh.
is_volume(), local_pts, bs, gbs);
471 const int tris_start = tris.size();
475 if (lb.type() == BoundaryType::QUAD)
477 const auto map = [n_samples, size](
int i,
int j) {
return j * n_samples + i + size; };
479 for (
int j = 0; j < n_samples - 1; ++j)
481 for (
int i = 0; i < n_samples - 1; ++i)
483 tris.emplace_back(map(i, j), map(i + 1, j), map(i, j + 1));
484 tris.emplace_back(map(i + 1, j + 1), map(i, j + 1), map(i + 1, j));
488 else if (lb.type() == BoundaryType::TRI)
491 std::vector<int> mapp(n_samples * n_samples, -1);
492 for (
int j = 0; j < n_samples; ++j)
494 for (
int i = 0; i < n_samples - j; ++i)
496 mapp[j * n_samples + i] = index;
500 const auto map = [mapp, n_samples](
int i,
int j) {
501 if (j * n_samples + i >= mapp.size())
503 return mapp[j * n_samples + i];
506 for (
int j = 0; j < n_samples - 1; ++j)
508 for (
int i = 0; i < n_samples - j; ++i)
510 if (map(i, j) >= 0 && map(i + 1, j) >= 0 && map(i, j + 1) >= 0)
511 tris.emplace_back(map(i, j) + size, map(i + 1, j) + size, map(i, j + 1) + size);
513 if (map(i + 1, j + 1) >= 0 && map(i, j + 1) >= 0 && map(i + 1, j) >= 0)
514 tris.emplace_back(map(i + 1, j + 1) + size, map(i, j + 1) + size, map(i + 1, j) + size);
525 for (
int i = 0; i < vertices.back().rows() - 1; ++i)
526 edges.emplace_back(i + size, i + size + 1);
529 normals.resize(
vals.jac_it.size(), tmp_n.cols());
530 displaced_normals.resize(
vals.jac_it.size(), tmp_n.cols());
532 for (
int n = 0; n <
vals.jac_it.size(); ++n)
534 trafo =
vals.jac_it[n].inverse();
536 if (problem_dim == 2 || problem_dim == 3)
539 if (solution.size() > 0)
541 deform_mat.resize(problem_dim, problem_dim);
542 deform_mat.setZero();
543 for (
const auto &b :
vals.basis_values)
544 for (
const auto &g : b.global)
545 for (
int d = 0; d < problem_dim; ++d)
546 deform_mat.row(d) += solution(g.index * problem_dim + d) * b.grad.row(n);
552 normals.row(n) = tmp_n *
vals.jac_it[n];
553 normals.row(n).normalize();
555 displaced_normals.row(n) = tmp_n * trafo.inverse();
556 displaced_normals.row(n).normalize();
559 allnormals.push_back(normals);
560 displaced_allnormals.push_back(displaced_normals);
563 for (
int n = 0; n <
vals.jac_it.size(); ++n)
565 tmp_n += normals.row(n);
570 Eigen::Vector3d e1 = vertices.back().row(std::get<1>(tris.back()) - size) - vertices.back().row(std::get<0>(tris.back()) - size);
571 Eigen::Vector3d e2 = vertices.back().row(std::get<2>(tris.back()) - size) - vertices.back().row(std::get<0>(tris.back()) - size);
573 Eigen::Vector3d n = e1.cross(e2);
574 Eigen::Vector3d nn = tmp_n.transpose();
578 for (
int i = tris_start; i < tris.size(); ++i)
580 tris[i] = std::tuple<int, int, int>(std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]));
585 size += vertices.back().rows();
589 boundary_vis_vertices.resize(size, vertices.front().cols());
590 boundary_vis_local_vertices.resize(size, vertices.front().cols());
591 boundary_vis_elements_ids.resize(size, 1);
592 boundary_vis_primitive_ids.resize(size, 1);
593 boundary_vis_normals.resize(size, vertices.front().cols());
594 displaced_boundary_vis_normals.resize(size, vertices.front().cols());
597 boundary_vis_elements.resize(tris.size(), 3);
599 boundary_vis_elements.resize(edges.size(), 2);
603 for (
const auto &v : vertices)
605 boundary_vis_vertices.block(index, 0, v.rows(), v.cols()) = v;
606 boundary_vis_local_vertices.block(index, 0, v.rows(), v.cols()) = lv[ii];
607 boundary_vis_elements_ids.block(index, 0, v.rows(), 1).setConstant(el_ids[ii]);
608 boundary_vis_primitive_ids.block(index, 0, v.rows(), 1).setConstant(global_primitive_ids[ii++]);
613 for (
const auto &n : allnormals)
615 boundary_vis_normals.block(index, 0, n.rows(), n.cols()) = n;
620 for (
const auto &n : displaced_allnormals)
622 displaced_boundary_vis_normals.block(index, 0, n.rows(), n.cols()) = n;
629 for (
const auto &t : tris)
631 boundary_vis_elements.row(index) << std::get<0>(t), std::get<1>(t), std::get<2>(t);
637 for (
const auto &e : edges)
639 boundary_vis_elements.row(index) << e.first, e.second;
647 const Eigen::VectorXi &disc_orders,
648 const std::vector<basis::ElementBases> &gbases,
649 const std::map<int, Eigen::MatrixXd> &polys,
650 const std::map<
int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
651 const bool boundary_only,
652 Eigen::MatrixXd &points,
653 Eigen::MatrixXi &tets,
654 Eigen::MatrixXi &el_id,
655 Eigen::MatrixXd &discr)
const
670 const auto ¤t_bases = gbases;
671 int tet_total_size = 0;
672 int pts_total_size = 0;
674 Eigen::MatrixXd vis_pts_poly;
675 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
677 for (
size_t i = 0; i < current_bases.size(); ++i)
679 const auto &bs = current_bases[i];
687 pts_total_size += sampler.simplex_points().rows();
691 tet_total_size += sampler.cube_volume().rows();
692 pts_total_size += sampler.cube_points().rows();
698 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
700 tet_total_size += vis_faces_poly.rows();
701 pts_total_size += vis_pts_poly.rows();
705 sampler.sample_polygon(polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
707 tet_total_size += vis_faces_poly.rows();
708 pts_total_size += vis_pts_poly.rows();
713 points.resize(pts_total_size, mesh.
dimension());
714 tets.resize(tet_total_size, mesh.
is_volume() ? 4 : 3);
716 el_id.resize(pts_total_size, 1);
717 discr.resize(pts_total_size, 1);
719 Eigen::MatrixXd mapped, tmp;
720 int tet_index = 0, pts_index = 0;
722 for (
size_t i = 0; i < current_bases.size(); ++i)
724 const auto &bs = current_bases[i];
731 bs.eval_geom_mapping(sampler.simplex_points(), mapped);
733 tets.block(tet_index, 0, sampler.simplex_volume().rows(), tets.cols()) = sampler.simplex_volume().array() + pts_index;
734 tet_index += sampler.simplex_volume().rows();
736 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
737 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
738 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
739 pts_index += mapped.rows();
743 bs.eval_geom_mapping(sampler.cube_points(), mapped);
745 tets.block(tet_index, 0, sampler.cube_volume().rows(), tets.cols()) = sampler.cube_volume().array() + pts_index;
746 tet_index += sampler.cube_volume().rows();
748 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
749 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
750 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
751 pts_index += mapped.rows();
757 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
758 bs.eval_geom_mapping(vis_pts_poly, mapped);
760 tets.block(tet_index, 0, vis_faces_poly.rows(), tets.cols()) = vis_faces_poly.array() + pts_index;
761 tet_index += vis_faces_poly.rows();
763 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
764 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(-1);
765 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
766 pts_index += mapped.rows();
770 sampler.sample_polygon(polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
771 bs.eval_geom_mapping(vis_pts_poly, mapped);
773 tets.block(tet_index, 0, vis_faces_poly.rows(), tets.cols()) = vis_faces_poly.array() + pts_index;
774 tet_index += vis_faces_poly.rows();
776 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
777 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(-1);
778 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
779 pts_index += mapped.rows();
784 assert(pts_index == points.rows());
785 assert(tet_index == tets.rows());
790 const Eigen::VectorXi &disc_orders,
791 const std::vector<basis::ElementBases> &bases,
792 Eigen::MatrixXd &points,
793 std::vector<std::vector<int>> &elements,
794 Eigen::MatrixXi &el_id,
795 Eigen::MatrixXd &discr)
const
809 std::vector<RowVectorNd> nodes;
810 int pts_total_size = 0;
811 elements.resize(bases.size());
812 Eigen::MatrixXd ref_pts;
814 for (
size_t i = 0; i < bases.size(); ++i)
816 const auto &bs = bases[i];
834 const int n_v =
static_cast<const mesh::Mesh2D &
>(mesh).n_face_vertices(i);
835 ref_pts.resize(n_v, 2);
839 pts_total_size += ref_pts.rows();
842 points.resize(pts_total_size, mesh.
dimension());
844 el_id.resize(pts_total_size, 1);
845 discr.resize(pts_total_size, 1);
847 Eigen::MatrixXd mapped;
850 std::string error_msg =
"";
852 for (
size_t i = 0; i < bases.size(); ++i)
854 const auto &bs = bases[i];
874 bs.eval_geom_mapping(ref_pts, mapped);
876 for (
int j = 0; j < mapped.rows(); ++j)
878 points.row(pts_index) = mapped.row(j);
879 el_id(pts_index) = i;
880 discr(pts_index) = disc_orders(i);
881 elements[i].push_back(pts_index);
890 const int n_nodes = elements[i].size();
891 if (disc_orders(i) >= 3)
893 std::swap(elements[i][16], elements[i][17]);
894 std::swap(elements[i][17], elements[i][18]);
895 std::swap(elements[i][18], elements[i][19]);
897 if (disc_orders(i) > 4)
898 error_msg =
"Saving high-order meshes not implemented for P5+ elements!";
902 if (disc_orders(i) == 4)
904 const int n_nodes = elements[i].size();
905 std::swap(elements[i][n_nodes - 1], elements[i][n_nodes - 2]);
907 if (disc_orders(i) > 4)
908 error_msg =
"Saving high-order meshes not implemented for P5+ elements!";
911 else if (disc_orders(i) > 1)
912 error_msg =
"Saving high-order meshes not implemented for Q2+ elements!";
915 if (!error_msg.empty())
918 for (
size_t i = 0; i < bases.size(); ++i)
923 const auto &mesh2d =
static_cast<const mesh::Mesh2D &
>(mesh);
926 for (
int j = 0; j < n_v; ++j)
928 points.row(pts_index) = mesh2d.point(mesh2d.face_vertex(i, j));
929 el_id(pts_index) = i;
930 discr(pts_index) = disc_orders(i);
931 elements[i].push_back(pts_index);
937 assert(pts_index == points.rows());
942 const Eigen::MatrixXd &sol,
943 const Eigen::MatrixXd &pressure,
944 const bool is_time_dependent,
945 const double tend_in,
948 const std::string &vis_mesh_path,
949 const std::string &nodes_path,
950 const std::string &solution_path,
951 const std::string &stress_path,
952 const std::string &mises_path,
953 const bool is_contact_enabled)
const
957 logger().error(
"Load the mesh first!");
960 const int n_bases = state.
n_bases;
961 const std::vector<basis::ElementBases> &bases = state.
bases;
962 const std::vector<basis::ElementBases> &gbases = state.
geom_bases();
965 const Eigen::MatrixXd &rhs = state.
rhs;
970 logger().error(
"Build the bases first!");
980 logger().error(
"Solve the problem first!");
984 if (!solution_path.empty())
986 std::ofstream out(solution_path);
988 out << std::scientific;
992 Eigen::VectorXi reordering(n_bases);
993 reordering.setConstant(-1);
995 for (
int i = 0; i < in_node_to_node.size(); ++i)
997 reordering[in_node_to_node[i]] = i;
1000 Eigen::MatrixXd tmp(tmp_sol.rows(), tmp_sol.cols());
1002 for (
int i = 0; i < reordering.size(); ++i)
1004 if (reordering[i] < 0)
1007 tmp.row(reordering[i]) = tmp_sol.row(i);
1010 for (
int i = 0; i < tmp.rows(); ++i)
1012 for (
int j = 0; j < tmp.cols(); ++j)
1013 out << tmp(i, j) <<
" ";
1019 out << sol << std::endl;
1023 double tend = tend_in;
1027 if (!vis_mesh_path.empty() && !is_time_dependent)
1030 vis_mesh_path, state, sol, pressure,
1032 is_contact_enabled);
1034 if (!nodes_path.empty())
1036 Eigen::MatrixXd nodes(n_bases, mesh.
dimension());
1042 for (
size_t ii = 0; ii < b.global().size(); ++ii)
1044 const auto &lg = b.global()[ii];
1045 nodes.row(lg.index) = lg.node;
1049 std::ofstream out(nodes_path);
1054 if (!stress_path.empty())
1056 Eigen::MatrixXd result;
1057 Eigen::VectorXd mises;
1061 sol, tend, result, mises);
1062 std::ofstream out(stress_path);
1066 if (!mises_path.empty())
1068 Eigen::MatrixXd result;
1069 Eigen::VectorXd mises;
1073 sol, tend, result, mises);
1074 std::ofstream out(mises_path);
1087 fields = args[
"output"][
"paraview"][
"fields"];
1089 volume = args[
"output"][
"paraview"][
"volume"];
1090 surface = args[
"output"][
"paraview"][
"surface"];
1091 wire = args[
"output"][
"paraview"][
"wireframe"];
1092 points = args[
"output"][
"paraview"][
"points"];
1093 contact_forces = args[
"output"][
"paraview"][
"options"][
"contact_forces"] && !is_problem_scalar;
1094 friction_forces = args[
"output"][
"paraview"][
"options"][
"friction_forces"] && !is_problem_scalar;
1095 normal_adhesion_forces = args[
"output"][
"paraview"][
"options"][
"normal_adhesion_forces"] && !is_problem_scalar;
1096 tangential_adhesion_forces = args[
"output"][
"paraview"][
"options"][
"tangential_adhesion_forces"] && !is_problem_scalar;
1098 if (args[
"output"][
"paraview"][
"options"][
"force_high_order"])
1099 use_sampler =
false;
1101 use_sampler = !(is_mesh_linear && args[
"output"][
"paraview"][
"high_order_mesh"]);
1102 boundary_only = use_sampler && args[
"output"][
"advanced"][
"vis_boundary_only"];
1103 material_params = args[
"output"][
"paraview"][
"options"][
"material"];
1104 body_ids = args[
"output"][
"paraview"][
"options"][
"body_ids"];
1105 sol_on_grid = args[
"output"][
"advanced"][
"sol_on_grid"] > 0;
1106 velocity = args[
"output"][
"paraview"][
"options"][
"velocity"];
1107 acceleration = args[
"output"][
"paraview"][
"options"][
"acceleration"];
1108 forces = args[
"output"][
"paraview"][
"options"][
"forces"] && !is_problem_scalar;
1109 jacobian_validity = args[
"output"][
"paraview"][
"options"][
"jacobian_validity"] && !is_problem_scalar;
1111 scalar_values = args[
"output"][
"paraview"][
"options"][
"scalar_values"];
1112 tensor_values = args[
"output"][
"paraview"][
"options"][
"tensor_values"] && !is_problem_scalar;
1113 discretization_order = args[
"output"][
"paraview"][
"options"][
"discretization_order"];
1114 nodes = args[
"output"][
"paraview"][
"options"][
"nodes"] && !is_problem_scalar;
1116 use_spline = args[
"space"][
"basis_type"] ==
"Spline";
1118 reorder_output = args[
"output"][
"data"][
"advanced"][
"reorder_nodes"];
1120 use_hdf5 = args[
"output"][
"paraview"][
"options"][
"use_hdf5"];
1124 const std::string &path,
1126 const Eigen::MatrixXd &sol,
1127 const Eigen::MatrixXd &pressure,
1131 const bool is_contact_enabled)
const
1135 logger().error(
"Load the mesh first!");
1139 const Eigen::MatrixXd &rhs = state.
rhs;
1143 logger().error(
"Build the bases first!");
1151 if (sol.size() <= 0)
1153 logger().error(
"Solve the problem first!");
1157 const std::filesystem::path fs_path(path);
1158 const std::string path_stem = fs_path.stem().string();
1159 const std::string base_path = (fs_path.parent_path() / path_stem).
string();
1169 is_contact_enabled);
1175 is_contact_enabled);
1188 paraviewo::VTMWriter vtm(t);
1190 vtm.add_dataset(
"Volume",
"data", path_stem + opts.
file_extension());
1192 vtm.add_dataset(
"Surface",
"data", path_stem +
"_surf" + opts.
file_extension());
1194 vtm.add_dataset(
"Contact",
"data", path_stem +
"_surf_contact" + opts.
file_extension());
1196 vtm.add_dataset(
"Wireframe",
"data", path_stem +
"_wire" + opts.
file_extension());
1198 vtm.add_dataset(
"Points",
"data", path_stem +
"_points" + opts.
file_extension());
1199 vtm.save(base_path +
".vtm");
1203 const std::string &path,
1205 const Eigen::MatrixXd &sol,
1206 const Eigen::MatrixXd &pressure,
1211 const Eigen::VectorXi &disc_orders = state.
disc_orders;
1213 const std::vector<basis::ElementBases> &bases = state.
bases;
1214 const std::vector<basis::ElementBases> &pressure_bases = state.
pressure_bases;
1215 const std::vector<basis::ElementBases> &gbases = state.
geom_bases();
1216 const std::map<int, Eigen::MatrixXd> &polys = state.
polys;
1217 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d = state.
polys_3d;
1224 Eigen::MatrixXd points;
1225 Eigen::MatrixXi tets;
1226 Eigen::MatrixXi el_id;
1227 Eigen::MatrixXd discr;
1228 std::vector<std::vector<int>> elements;
1233 points, tets, el_id, discr);
1236 points, elements, el_id, discr);
1238 Eigen::MatrixXd fun, exact_fun, err, node_fun;
1243 Eigen::MatrixXd tmp, tmp_grad;
1244 Eigen::MatrixXd tmp_p, tmp_grad_p;
1246 res.setConstant(std::numeric_limits<double>::quiet_NaN());
1248 res_grad.setConstant(std::numeric_limits<double>::quiet_NaN());
1251 res_p.setConstant(std::numeric_limits<double>::quiet_NaN());
1253 res_grad_p.setConstant(std::numeric_limits<double>::quiet_NaN());
1262 Eigen::MatrixXd pt(1, bc.cols() - 1);
1263 for (
int d = 1; d < bc.cols(); ++d)
1266 mesh, problem.
is_scalar(), bases, gbases,
1267 el_id, pt, sol, tmp, tmp_grad);
1270 res_grad.row(i) = tmp_grad;
1275 mesh, 1, pressure_bases, gbases,
1276 el_id, pt, pressure, tmp_p, tmp_grad_p);
1277 res_p.row(i) = tmp_p;
1278 res_grad_p.row(i) = tmp_grad_p;
1282 std::ofstream os(path +
"_sol.txt");
1285 std::ofstream osg(path +
"_grad.txt");
1288 std::ofstream osgg(path +
"_grid.txt");
1293 std::ofstream osp(path +
"_p_sol.txt");
1296 std::ofstream osgp(path +
"_p_grad.txt");
1301 Eigen::Vector<bool, -1> validity;
1314 Eigen::MatrixXd tmp = Eigen::VectorXd::LinSpaced(sol.size(), 0, sol.size() - 1);
1324 fun.conservativeResize(fun.rows() + obstacle.
n_vertices(), fun.cols());
1325 node_fun.conservativeResize(node_fun.rows() + obstacle.
n_vertices(), node_fun.cols());
1326 node_fun.bottomRows(obstacle.
n_vertices()).setZero();
1334 problem.
exact(points, t, exact_fun);
1335 err = (fun - exact_fun).eval().rowwise().norm();
1339 exact_fun.conservativeResize(exact_fun.rows() + obstacle.
n_vertices(), exact_fun.cols());
1343 err.conservativeResize(err.rows() + obstacle.
n_vertices(), 1);
1344 err.bottomRows(obstacle.
n_vertices()).setZero();
1348 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
1350 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
1352 tmpw = std::make_shared<paraviewo::VTUWriter>();
1353 paraviewo::ParaviewWriter &writer = *tmpw;
1356 writer.add_field(
"validity", validity.cast<
double>());
1359 writer.add_field(
"nodes", node_fun);
1363 bool is_time_integrator_valid = time_integrator !=
nullptr;
1367 const Eigen::VectorXd velocity =
1368 is_time_integrator_valid ? (time_integrator->v_prev()) : Eigen::VectorXd::Zero(sol.size());
1374 const Eigen::VectorXd acceleration =
1375 is_time_integrator_valid ? (time_integrator->a_prev()) : Eigen::VectorXd::Zero(sol.size());
1389 if (form ==
nullptr)
1392 Eigen::VectorXd force;
1393 if (form->enabled())
1395 form->first_derivative(sol, force);
1400 force.setZero(sol.size());
1410 Eigen::MatrixXd interp_p;
1418 interp_p.conservativeResize(interp_p.size() + obstacle.
n_vertices(), 1);
1419 interp_p.bottomRows(obstacle.
n_vertices()).setZero();
1422 writer.add_field(
"pressure", interp_p);
1427 discr.conservativeResize(discr.size() + obstacle.
n_vertices(), 1);
1428 discr.bottomRows(obstacle.
n_vertices()).setZero();
1432 writer.add_field(
"discr", discr);
1437 writer.add_field(
"exact", exact_fun);
1439 writer.add_field(
"error", err);
1442 if (fun.cols() != 1)
1444 std::vector<assembler::Assembler::NamedMatrix>
vals, tvals;
1446 mesh, problem.
is_scalar(), bases, gbases,
1451 for (
auto &[_, v] :
vals)
1456 for (
const auto &[name, v] :
vals)
1459 writer.add_field(name, v);
1470 for (
auto &[_, v] : tvals)
1473 for (
const auto &[name, v] : tvals)
1476 assert(v.cols() % stride == 0);
1481 for (
int i = 0; i < v.cols(); i += stride)
1483 const Eigen::MatrixXd tmp = v.middleCols(i, stride);
1484 assert(tmp.cols() == stride);
1486 const int ii = (i / stride) + 1;
1487 writer.add_field(fmt::format(
"{:s}_{:d}", name, ii), tmp);
1502 for (
auto &v :
vals)
1504 v.second.conservativeResize(v.second.size() + obstacle.
n_vertices(), 1);
1505 v.second.bottomRows(obstacle.
n_vertices()).setZero();
1511 for (
const auto &v :
vals)
1513 if (opts.
export_field(fmt::format(
"{:s}_avg", v.first)))
1514 writer.add_field(fmt::format(
"{:s}_avg", v.first), v.second);
1524 std::map<std::string, Eigen::MatrixXd> param_val;
1525 for (
const auto &[p, _] : params)
1526 param_val[p] = Eigen::MatrixXd(points.rows(), 1);
1527 Eigen::MatrixXd rhos(points.rows(), 1);
1529 Eigen::MatrixXd local_pts;
1530 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
1534 for (
int e = 0; e < int(bases.size()); ++e)
1542 local_pts = sampler.simplex_points();
1544 local_pts = sampler.cube_points();
1548 sampler.sample_polyhedron(polys_3d.at(e).first, polys_3d.at(e).second, local_pts, vis_faces_poly, vis_edges_poly);
1550 sampler.sample_polygon(polys.at(e), local_pts, vis_faces_poly, vis_edges_poly);
1572 const auto &mesh2d =
static_cast<const mesh::Mesh2D &
>(mesh);
1574 local_pts.resize(n_v, 2);
1576 for (
int j = 0; j < n_v; ++j)
1578 local_pts.row(j) = mesh2d.point(mesh2d.face_vertex(e, j));
1587 for (
int j = 0; j <
vals.val.rows(); ++j)
1589 for (
const auto &[p, func] : params)
1590 param_val.at(p)(index) = func(local_pts.row(j),
vals.val.row(j), t, e);
1592 rhos(index) = density(local_pts.row(j),
vals.val.row(j), t, e);
1598 assert(index == points.rows());
1602 for (
auto &[_, tmp] : param_val)
1604 tmp.conservativeResize(tmp.size() + obstacle.
n_vertices(), 1);
1605 tmp.bottomRows(obstacle.
n_vertices()).setZero();
1608 rhos.conservativeResize(rhos.size() + obstacle.
n_vertices(), 1);
1609 rhos.bottomRows(obstacle.
n_vertices()).setZero();
1611 for (
const auto &[p, tmp] : param_val)
1614 writer.add_field(p, tmp);
1617 writer.add_field(
"rho", rhos);
1623 Eigen::MatrixXd ids(points.rows(), 1);
1625 for (
int i = 0; i < points.rows(); ++i)
1632 ids.conservativeResize(ids.size() + obstacle.
n_vertices(), 1);
1633 ids.bottomRows(obstacle.
n_vertices()).setZero();
1636 writer.add_field(
"body_ids", ids);
1647 Eigen::MatrixXd traction_forces, traction_forces_fun;
1648 compute_traction_forces(state, sol, t, traction_forces,
false);
1657 traction_forces_fun.conservativeResize(traction_forces_fun.rows() + obstacle.
n_vertices(), traction_forces_fun.cols());
1658 traction_forces_fun.bottomRows(obstacle.
n_vertices()).setZero();
1661 writer.add_field(
"traction_force", traction_forces_fun);
1668 Eigen::MatrixXd potential_grad, potential_grad_fun;
1678 potential_grad_fun.conservativeResize(potential_grad_fun.rows() + obstacle.
n_vertices(), potential_grad_fun.cols());
1679 potential_grad_fun.bottomRows(obstacle.
n_vertices()).setZero();
1681 writer.add_field(
"gradient_of_potential", potential_grad_fun);
1683 catch (std::exception &)
1689 writer.add_field(
"solution", fun);
1693 const int orig_p = points.rows();
1694 points.conservativeResize(points.rows() + obstacle.
n_vertices(), points.cols());
1695 points.bottomRows(obstacle.
n_vertices()) = obstacle.
v();
1697 if (elements.empty())
1699 for (
int i = 0; i < tets.rows(); ++i)
1701 elements.emplace_back();
1702 for (
int j = 0; j < tets.cols(); ++j)
1703 elements.back().push_back(tets(i, j));
1709 elements.emplace_back();
1716 elements.emplace_back();
1723 elements.emplace_back();
1728 if (elements.empty())
1729 writer.write_mesh(path, points, tets);
1731 writer.write_mesh(path, points, elements,
true, disc_orders.maxCoeff() == 1);
1736 const Eigen::MatrixXd &points,
1738 const std::string &name,
1739 const Eigen::VectorXd &field,
1740 paraviewo::ParaviewWriter &writer)
const
1742 Eigen::MatrixXd inerpolated_field;
1750 inerpolated_field.conservativeResize(
1756 writer.add_field(name, inerpolated_field);
1760 const std::string &export_surface,
1762 const Eigen::MatrixXd &sol,
1763 const Eigen::MatrixXd &pressure,
1767 const bool is_contact_enabled)
const
1770 const Eigen::VectorXi &disc_orders = state.
disc_orders;
1772 const std::vector<basis::ElementBases> &bases = state.
bases;
1773 const std::vector<basis::ElementBases> &pressure_bases = state.
pressure_bases;
1774 const std::vector<basis::ElementBases> &gbases = state.
geom_bases();
1780 Eigen::MatrixXd boundary_vis_vertices;
1781 Eigen::MatrixXd boundary_vis_local_vertices;
1782 Eigen::MatrixXi boundary_vis_elements;
1783 Eigen::MatrixXi boundary_vis_elements_ids;
1784 Eigen::MatrixXi boundary_vis_primitive_ids;
1785 Eigen::MatrixXd boundary_vis_normals;
1786 Eigen::MatrixXd displaced_boundary_vis_normals;
1789 boundary_vis_vertices, boundary_vis_local_vertices, boundary_vis_elements,
1790 boundary_vis_elements_ids, boundary_vis_primitive_ids, boundary_vis_normals,
1791 displaced_boundary_vis_normals);
1793 Eigen::MatrixXd fun, interp_p, discr, vect, b_sidesets;
1795 Eigen::MatrixXd lsol, lp, lgrad, lpgrad;
1801 discr.resize(boundary_vis_vertices.rows(), 1);
1802 fun.resize(boundary_vis_vertices.rows(), actual_dim);
1803 interp_p.resize(boundary_vis_vertices.rows(), 1);
1804 vect.resize(boundary_vis_vertices.rows(), mesh.
dimension());
1806 b_sidesets.resize(boundary_vis_vertices.rows(), 1);
1807 b_sidesets.setZero();
1809 for (
int i = 0; i < boundary_vis_vertices.rows(); ++i)
1811 const auto s_id = mesh.
get_boundary_id(boundary_vis_primitive_ids(i));
1814 b_sidesets(i) = s_id;
1817 const int el_index = boundary_vis_elements_ids(i);
1819 mesh, problem.
is_scalar(), bases, gbases,
1820 el_index, boundary_vis_local_vertices.row(i), sol, lsol, lgrad);
1821 assert(lsol.size() == actual_dim);
1825 mesh, 1, pressure_bases, gbases,
1826 el_index, boundary_vis_local_vertices.row(i), pressure, lp, lpgrad);
1827 assert(lp.size() == 1);
1828 interp_p(i) = lp(0);
1831 discr(i) = disc_orders(el_index);
1832 for (
int j = 0; j < actual_dim; ++j)
1834 fun(i, j) = lsol(j);
1837 if (actual_dim == 1)
1839 assert(lgrad.size() == mesh.
dimension());
1840 for (
int j = 0; j < mesh.
dimension(); ++j)
1842 vect(i, j) = lgrad(j);
1847 assert(lgrad.size() == actual_dim * actual_dim);
1848 std::vector<assembler::Assembler::NamedMatrix> tensor_flat;
1853 assert(tensor_flat[0].first ==
"cauchy_stess");
1854 assert(tensor_flat[0].second.size() == actual_dim * actual_dim);
1856 Eigen::Map<Eigen::MatrixXd> tensor(tensor_flat[0].second.data(), actual_dim, actual_dim);
1857 vect.row(i) = displaced_boundary_vis_normals.row(i) * tensor;
1863 area = mesh.
tri_area(boundary_vis_primitive_ids(i));
1864 else if (mesh.
is_cube(el_index))
1865 area = mesh.
quad_area(boundary_vis_primitive_ids(i));
1868 area = mesh.
edge_length(boundary_vis_primitive_ids(i));
1870 vect.row(i) *= area;
1874 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
1876 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
1878 tmpw = std::make_shared<paraviewo::VTUWriter>();
1879 paraviewo::ParaviewWriter &writer = *tmpw;
1882 writer.add_field(
"normals", boundary_vis_normals);
1884 writer.add_field(
"displaced_normals", displaced_boundary_vis_normals);
1886 writer.add_field(
"pressure", interp_p);
1888 writer.add_field(
"discr", discr);
1890 writer.add_field(
"sidesets", b_sidesets);
1892 if (actual_dim == 1 && opts.
export_field(
"solution_grad"))
1893 writer.add_field(
"solution_grad", vect);
1896 writer.add_field(
"traction_force", vect);
1903 std::map<std::string, Eigen::MatrixXd> param_val;
1904 for (
const auto &[p, _] : params)
1905 param_val[p] = Eigen::MatrixXd(boundary_vis_vertices.rows(), 1);
1906 Eigen::MatrixXd rhos(boundary_vis_vertices.rows(), 1);
1908 for (
int i = 0; i < boundary_vis_vertices.rows(); ++i)
1912 for (
const auto &[p, func] : params)
1913 param_val.at(p)(i) = func(boundary_vis_local_vertices.row(i), boundary_vis_vertices.row(i), t, boundary_vis_elements_ids(i));
1915 rhos(i) = density(boundary_vis_local_vertices.row(i), boundary_vis_vertices.row(i), t, boundary_vis_elements_ids(i));
1918 for (
const auto &[p, tmp] : param_val)
1921 writer.add_field(p, tmp);
1924 writer.add_field(
"rho", rhos);
1930 Eigen::MatrixXd ids(boundary_vis_vertices.rows(), 1);
1932 for (
int i = 0; i < boundary_vis_vertices.rows(); ++i)
1934 ids(i) = mesh.
get_body_id(boundary_vis_elements_ids(i));
1937 writer.add_field(
"body_ids", ids);
1941 writer.add_field(
"solution", fun);
1942 writer.write_mesh(export_surface, boundary_vis_vertices, boundary_vis_elements);
1946 const std::string &export_surface,
1948 const Eigen::MatrixXd &sol,
1949 const Eigen::MatrixXd &pressure,
1953 const bool is_contact_enabled)
const
1957 const double dhat = state.
args[
"contact"][
"dhat"];
1958 const double friction_coefficient = state.
args[
"contact"][
"friction_coefficient"];
1959 const double epsv = state.
args[
"contact"][
"epsv"];
1960 const double dhat_a = state.
args[
"contact"][
"adhesion"][
"dhat_a"];
1961 const double dhat_p = state.
args[
"contact"][
"adhesion"][
"dhat_p"];
1962 const double Y = state.
args[
"contact"][
"adhesion"][
"adhesion_strength"];
1963 const double epsa = state.
args[
"contact"][
"adhesion"][
"epsa"];
1964 const double tangential_adhesion_coefficient = state.
args[
"contact"][
"adhesion"][
"tangential_adhesion_coefficient"];
1970 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
1972 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
1974 tmpw = std::make_shared<paraviewo::VTUWriter>();
1975 paraviewo::ParaviewWriter &writer = *tmpw;
1977 const int problem_dim = mesh.
dimension();
1978 const Eigen::MatrixXd full_displacements =
utils::unflatten(sol, problem_dim);
1979 const Eigen::MatrixXd surface_displacements = collision_mesh.map_displacements(full_displacements);
1981 const Eigen::MatrixXd displaced_surface = collision_mesh.displace_vertices(full_displacements);
1983 ipc::NormalCollisions collision_set;
1985 if (state.
args[
"contact"][
"use_convergent_formulation"])
1987 collision_set.set_use_improved_max_approximator(state.
args[
"contact"][
"use_improved_max_operator"]);
1988 collision_set.set_use_area_weighting(state.
args[
"contact"][
"use_area_weighting"]);
1991 collision_set.build(
1992 collision_mesh, displaced_surface, dhat,
1993 0, state.
args[
"solver"][
"contact"][
"CCD"][
"broad_phase"]);
1995 ipc::BarrierPotential barrier_potential(dhat);
1996 if (state.
args[
"contact"][
"use_convergent_formulation"])
1998 barrier_potential.set_use_physical_barrier(state.
args[
"contact"][
"use_physical_barrier"]);
2001 const double barrier_stiffness = contact_form !=
nullptr ? contact_form->barrier_stiffness() : 1;
2005 Eigen::MatrixXd forces = -barrier_stiffness * barrier_potential.gradient(collision_set, collision_mesh, displaced_surface);
2009 assert(forces_reshaped.rows() == surface_displacements.rows());
2010 assert(forces_reshaped.cols() == surface_displacements.cols());
2011 writer.add_field(
"contact_forces", forces_reshaped);
2016 ipc::TangentialCollisions friction_collision_set;
2017 friction_collision_set.build(
2018 collision_mesh, displaced_surface, collision_set,
2019 barrier_potential, barrier_stiffness, friction_coefficient);
2021 ipc::FrictionPotential friction_potential(epsv);
2023 Eigen::MatrixXd velocities;
2028 velocities = collision_mesh.map_displacements(
utils::unflatten(velocities, collision_mesh.dim()));
2030 Eigen::MatrixXd forces = -friction_potential.gradient(
2031 friction_collision_set, collision_mesh, velocities);
2035 assert(forces_reshaped.rows() == surface_displacements.rows());
2036 assert(forces_reshaped.cols() == surface_displacements.cols());
2037 writer.add_field(
"friction_forces", forces_reshaped);
2040 ipc::NormalCollisions adhesion_collision_set;
2041 adhesion_collision_set.build(
2042 collision_mesh, displaced_surface, dhat_a,
2043 0, state.
args[
"solver"][
"contact"][
"CCD"][
"broad_phase"]);
2045 ipc::NormalAdhesionPotential normal_adhesion_potential(dhat_p, dhat_a, Y, 1);
2049 Eigen::MatrixXd forces = -1 * normal_adhesion_potential.gradient(adhesion_collision_set, collision_mesh, displaced_surface);
2053 assert(forces_reshaped.rows() == surface_displacements.rows());
2054 assert(forces_reshaped.cols() == surface_displacements.cols());
2055 writer.add_field(
"normal_adhesion_forces", forces_reshaped);
2060 ipc::TangentialCollisions tangential_collision_set;
2061 tangential_collision_set.build(
2062 collision_mesh, displaced_surface, adhesion_collision_set,
2063 normal_adhesion_potential, 1, tangential_adhesion_coefficient);
2065 ipc::TangentialAdhesionPotential tangential_adhesion_potential(epsa);
2067 Eigen::MatrixXd velocities;
2072 velocities = collision_mesh.map_displacements(
utils::unflatten(velocities, collision_mesh.dim()));
2074 Eigen::MatrixXd forces = -tangential_adhesion_potential.gradient(
2075 tangential_collision_set, collision_mesh, velocities);
2079 assert(forces_reshaped.rows() == surface_displacements.rows());
2080 assert(forces_reshaped.cols() == surface_displacements.cols());
2081 writer.add_field(
"tangential_adhesion_forces", forces_reshaped);
2084 assert(collision_mesh.rest_positions().rows() == surface_displacements.rows());
2085 assert(collision_mesh.rest_positions().cols() == surface_displacements.cols());
2088 writer.add_field(
"solution", surface_displacements);
2091 export_surface.substr(0, export_surface.length() - 4) +
"_contact.vtu",
2092 collision_mesh.rest_positions(),
2093 problem_dim == 3 ? collision_mesh.faces() : collision_mesh.edges());
2097 const std::string &name,
2099 const Eigen::MatrixXd &sol,
2103 const std::vector<basis::ElementBases> &gbases = state.
geom_bases();
2109 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
2110 Eigen::MatrixXd vis_pts_poly;
2112 const auto ¤t_bases = gbases;
2113 int seg_total_size = 0;
2114 int pts_total_size = 0;
2115 int faces_total_size = 0;
2117 for (
size_t i = 0; i < current_bases.size(); ++i)
2119 const auto &bs = current_bases[i];
2124 seg_total_size += sampler.simplex_edges().rows();
2125 faces_total_size += sampler.simplex_faces().rows();
2129 pts_total_size += sampler.cube_points().rows();
2130 seg_total_size += sampler.cube_edges().rows();
2131 faces_total_size += sampler.cube_faces().rows();
2136 sampler.sample_polyhedron(state.
polys_3d.at(i).first, state.
polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
2138 sampler.sample_polygon(state.
polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
2140 pts_total_size += vis_pts_poly.rows();
2141 seg_total_size += vis_edges_poly.rows();
2142 faces_total_size += vis_faces_poly.rows();
2146 Eigen::MatrixXd points(pts_total_size, mesh.
dimension());
2147 Eigen::MatrixXi edges(seg_total_size, 2);
2148 Eigen::MatrixXi
faces(faces_total_size, 3);
2151 Eigen::MatrixXd mapped, tmp;
2152 int seg_index = 0, pts_index = 0, face_index = 0;
2153 for (
size_t i = 0; i < current_bases.size(); ++i)
2155 const auto &bs = current_bases[i];
2159 bs.eval_geom_mapping(sampler.simplex_points(), mapped);
2160 edges.block(seg_index, 0, sampler.simplex_edges().rows(), edges.cols()) = sampler.simplex_edges().array() + pts_index;
2161 seg_index += sampler.simplex_edges().rows();
2163 faces.block(face_index, 0, sampler.simplex_faces().rows(), 3) = sampler.simplex_faces().array() + pts_index;
2164 face_index += sampler.simplex_faces().rows();
2166 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2167 pts_index += mapped.rows();
2171 bs.eval_geom_mapping(sampler.cube_points(), mapped);
2172 edges.block(seg_index, 0, sampler.cube_edges().rows(), edges.cols()) = sampler.cube_edges().array() + pts_index;
2173 seg_index += sampler.cube_edges().rows();
2175 faces.block(face_index, 0, sampler.cube_faces().rows(), 3) = sampler.cube_faces().array() + pts_index;
2176 face_index += sampler.cube_faces().rows();
2178 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2179 pts_index += mapped.rows();
2184 sampler.sample_polyhedron(state.
polys_3d.at(i).first, state.
polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
2186 sampler.sample_polygon(state.
polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
2188 edges.block(seg_index, 0, vis_edges_poly.rows(), edges.cols()) = vis_edges_poly.array() + pts_index;
2189 seg_index += vis_edges_poly.rows();
2191 faces.block(face_index, 0, vis_faces_poly.rows(), 3) = vis_faces_poly.array() + pts_index;
2192 face_index += vis_faces_poly.rows();
2194 points.block(pts_index, 0, vis_pts_poly.rows(), points.cols()) = vis_pts_poly;
2195 pts_index += vis_pts_poly.rows();
2199 assert(pts_index == points.rows());
2200 assert(face_index ==
faces.rows());
2205 for (
long i = 0; i <
faces.rows(); ++i)
2207 const int v0 =
faces(i, 0);
2208 const int v1 =
faces(i, 1);
2209 const int v2 =
faces(i, 2);
2211 int tmpc =
faces(i, 2);
2218 Eigen::Matrix2d mmat;
2219 for (
long i = 0; i <
faces.rows(); ++i)
2221 const int v0 =
faces(i, 0);
2222 const int v1 =
faces(i, 1);
2223 const int v2 =
faces(i, 2);
2225 mmat.row(0) = points.row(v2) - points.row(v0);
2226 mmat.row(1) = points.row(v1) - points.row(v0);
2228 if (mmat.determinant() > 0)
2230 int tmpc =
faces(i, 2);
2237 Eigen::MatrixXd fun;
2241 pts_index, sol, fun,
true,
false);
2243 Eigen::MatrixXd exact_fun, err;
2247 problem.
exact(points, t, exact_fun);
2248 err = (fun - exact_fun).eval().rowwise().norm();
2251 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2253 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2255 tmpw = std::make_shared<paraviewo::VTUWriter>();
2256 paraviewo::ParaviewWriter &writer = *tmpw;
2261 writer.add_field(
"exact", exact_fun);
2263 writer.add_field(
"error", err);
2266 if (fun.cols() != 1)
2268 std::vector<assembler::Assembler::NamedMatrix> scalar_val;
2274 for (
const auto &v : scalar_val)
2277 writer.add_field(v.first, v.second);
2281 writer.add_field(
"solution", fun);
2283 writer.write_mesh(name, points, edges);
2287 const std::string &path,
2289 const Eigen::MatrixXd &sol,
2301 Eigen::MatrixXd fun(dirichlet_nodes_position.size(), actual_dim);
2302 Eigen::MatrixXd b_sidesets(dirichlet_nodes_position.size(), 1);
2303 b_sidesets.setZero();
2304 Eigen::MatrixXd points(dirichlet_nodes_position.size(), mesh.
dimension());
2305 std::vector<std::vector<int>> cells(dirichlet_nodes_position.size());
2307 for (
int i = 0; i < dirichlet_nodes_position.size(); ++i)
2309 const int n_id = dirichlet_nodes[i];
2313 b_sidesets(i) = s_id;
2316 for (
int j = 0; j < actual_dim; ++j)
2318 fun(i, j) = sol(n_id * actual_dim + j);
2321 points.row(i) = dirichlet_nodes_position[i];
2322 cells[i].push_back(i);
2325 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2327 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2329 tmpw = std::make_shared<paraviewo::VTUWriter>();
2330 paraviewo::ParaviewWriter &writer = *tmpw;
2333 writer.add_field(
"sidesets", b_sidesets);
2335 writer.add_field(
"solution", fun);
2336 writer.write_mesh(path, points, cells,
false,
false);
2340 const std::string &name,
2341 const std::function<std::string(
int)> &vtu_names,
2342 int time_steps,
double t0,
double dt,
int skip_frame)
const
2344 paraviewo::PVDWriter::save_pvd(name, vtu_names, time_steps, t0, dt, skip_frame);
2360 const int nx = delta[0] / spacing + 1;
2361 const int ny = delta[1] / spacing + 1;
2362 const int nz = delta.cols() >= 3 ? (delta[2] / spacing + 1) : 1;
2363 const int n = nx * ny * nz;
2367 for (
int i = 0; i < nx; ++i)
2369 const double x = (delta[0] / (nx - 1)) * i + min[0];
2371 for (
int j = 0; j < ny; ++j)
2373 const double y = (delta[1] / (ny - 1)) * j + min[1];
2375 if (delta.cols() <= 2)
2381 for (
int k = 0; k < nz; ++k)
2383 const double z = (delta[2] / (nz - 1)) * k + min[2];
2392 std::vector<std::array<Eigen::Vector3d, 2>> boxes;
2398 const double eps = 1e-6;
2407 const Eigen::Vector3d min(
2412 const Eigen::Vector3d max(
2417 std::vector<unsigned int> candidates;
2419 bvh.intersect_box(min, max, candidates);
2421 for (
const auto cand : candidates)
2425 logger().warn(
"Element {} is not simplex, skipping", cand);
2429 Eigen::MatrixXd coords;
2432 for (
int d = 0; d < coords.size(); ++d)
2434 if (fabs(coords(d)) < 1e-8)
2436 else if (fabs(coords(d) - 1) < 1e-8)
2440 if (coords.array().minCoeff() >= 0 && coords.array().maxCoeff() <= 1)
2452 Eigen::MatrixXd samples_simplex, samples_cube, mapped, p0, p1, p;
2455 average_edge_length = 0;
2456 min_edge_length = std::numeric_limits<double>::max();
2458 if (!use_curved_mesh_size)
2462 min_edge_length = p.rowwise().norm().minCoeff();
2463 average_edge_length = p.rowwise().norm().mean();
2464 mesh_size = p.rowwise().norm().maxCoeff();
2466 logger().info(
"hmin: {}", min_edge_length);
2467 logger().info(
"hmax: {}", mesh_size);
2468 logger().info(
"havg: {}", average_edge_length);
2485 for (
size_t i = 0; i < bases_in.size(); ++i)
2494 bases_in[i].eval_geom_mapping(samples_simplex, mapped);
2499 bases_in[i].eval_geom_mapping(samples_cube, mapped);
2502 for (
int j = 0; j < n_edges; ++j)
2504 double current_edge = 0;
2505 for (
int k = 0; k < n_samples - 1; ++k)
2507 p0 = mapped.row(j * n_samples + k);
2508 p1 = mapped.row(j * n_samples + k + 1);
2511 current_edge += p.norm();
2514 mesh_size = std::max(current_edge, mesh_size);
2515 min_edge_length = std::min(current_edge, min_edge_length);
2516 average_edge_length += current_edge;
2521 average_edge_length /= n;
2523 logger().info(
"hmin: {}", min_edge_length);
2524 logger().info(
"hmax: {}", mesh_size);
2525 logger().info(
"havg: {}", average_edge_length);
2539 using namespace mesh;
2541 logger().info(
"Counting flipped elements...");
2545 for (
size_t i = 0; i < gbases.size(); ++i)
2551 if (!
vals.is_geom_mapping_positive(mesh.
is_volume(), gbases[i]))
2555 static const std::vector<std::string> element_type_names{{
2557 "RegularInteriorCube",
2558 "RegularBoundaryCube",
2559 "SimpleSingularInteriorCube",
2560 "MultiSingularInteriorCube",
2561 "SimpleSingularBoundaryCube",
2563 "MultiSingularBoundaryCube",
2569 log_and_throw_error(
"element {} is flipped, type {}", i, element_type_names[
static_cast<int>(els_tag[i])]);
2584 const std::vector<polyfem::basis::ElementBases> &bases,
2585 const std::vector<polyfem::basis::ElementBases> &gbases,
2589 const Eigen::MatrixXd &sol)
2593 logger().error(
"Build the bases first!");
2596 if (sol.size() <= 0)
2598 logger().error(
"Solve the problem first!");
2608 logger().info(
"Computing errors...");
2611 const int n_el = int(bases.size());
2613 Eigen::MatrixXd v_exact, v_approx;
2614 Eigen::MatrixXd v_exact_grad(0, 0), v_approx_grad;
2624 static const int p = 8;
2629 for (
int e = 0; e < n_el; ++e)
2639 v_approx.resize(
vals.val.rows(), actual_dim);
2642 v_approx_grad.resize(
vals.val.rows(), mesh.
dimension() * actual_dim);
2643 v_approx_grad.setZero();
2645 const int n_loc_bases = int(
vals.basis_values.size());
2647 for (
int i = 0; i < n_loc_bases; ++i)
2649 const auto &
val =
vals.basis_values[i];
2651 for (
size_t ii = 0; ii <
val.global.size(); ++ii)
2653 for (
int d = 0; d < actual_dim; ++d)
2655 v_approx.col(d) +=
val.global[ii].val * sol(
val.global[ii].index * actual_dim + d) *
val.val;
2656 v_approx_grad.block(0, d *
val.grad_t_m.cols(), v_approx_grad.rows(),
val.grad_t_m.cols()) +=
val.global[ii].val * sol(
val.global[ii].index * actual_dim + d) *
val.grad_t_m;
2661 const auto err = problem.
has_exact_sol() ? (v_exact - v_approx).eval().rowwise().norm().eval() : (v_approx).eval().rowwise().norm().eval();
2662 const auto err_grad = problem.
has_exact_sol() ? (v_exact_grad - v_approx_grad).eval().rowwise().norm().eval() : (v_approx_grad).eval().rowwise().norm().eval();
2667 linf_err = std::max(linf_err, err.maxCoeff());
2668 grad_max_err = std::max(linf_err, err_grad.maxCoeff());
2710 l2_err += (err.array() * err.array() *
vals.det.array() *
vals.quadrature.weights.array()).sum();
2711 h1_err += (err_grad.array() * err_grad.array() *
vals.det.array() *
vals.quadrature.weights.array()).sum();
2712 lp_err += (err.array().pow(p) *
vals.det.array() *
vals.quadrature.weights.array()).sum();
2715 h1_semi_err = sqrt(fabs(h1_err));
2716 h1_err = sqrt(fabs(l2_err) + fabs(h1_err));
2717 l2_err = sqrt(fabs(l2_err));
2719 lp_err = pow(fabs(lp_err), 1. / p);
2724 const double computing_errors_time = timer.getElapsedTime();
2725 logger().info(
" took {}s", computing_errors_time);
2727 logger().info(
"-- L2 error: {}", l2_err);
2728 logger().info(
"-- Lp error: {}", lp_err);
2729 logger().info(
"-- H1 error: {}", h1_err);
2730 logger().info(
"-- H1 semi error: {}", h1_semi_err);
2733 logger().info(
"-- Linf error: {}", linf_err);
2734 logger().info(
"-- grad max error: {}", grad_max_err);
2749 regular_boundary_count = 0;
2750 simple_singular_count = 0;
2751 multi_singular_count = 0;
2753 non_regular_boundary_count = 0;
2754 non_regular_count = 0;
2755 undefined_count = 0;
2756 multi_singular_boundary_count = 0;
2760 for (
size_t i = 0; i < els_tag.size(); ++i)
2766 case ElementType::SIMPLEX:
2769 case ElementType::REGULAR_INTERIOR_CUBE:
2772 case ElementType::REGULAR_BOUNDARY_CUBE:
2773 regular_boundary_count++;
2775 case ElementType::SIMPLE_SINGULAR_INTERIOR_CUBE:
2776 simple_singular_count++;
2778 case ElementType::MULTI_SINGULAR_INTERIOR_CUBE:
2779 multi_singular_count++;
2781 case ElementType::SIMPLE_SINGULAR_BOUNDARY_CUBE:
2784 case ElementType::INTERFACE_CUBE:
2785 case ElementType::MULTI_SINGULAR_BOUNDARY_CUBE:
2786 multi_singular_boundary_count++;
2788 case ElementType::BOUNDARY_POLYTOPE:
2789 non_regular_boundary_count++;
2791 case ElementType::INTERIOR_POLYTOPE:
2792 non_regular_count++;
2794 case ElementType::UNDEFINED:
2798 throw std::runtime_error(
"Unknown element type");
2802 logger().info(
"simplex_count: \t{}", simplex_count);
2803 logger().info(
"regular_count: \t{}", regular_count);
2804 logger().info(
"regular_boundary_count: \t{}", regular_boundary_count);
2805 logger().info(
"simple_singular_count: \t{}", simple_singular_count);
2806 logger().info(
"multi_singular_count: \t{}", multi_singular_count);
2807 logger().info(
"boundary_count: \t{}", boundary_count);
2808 logger().info(
"multi_singular_boundary_count: \t{}", multi_singular_boundary_count);
2809 logger().info(
"non_regular_count: \t{}", non_regular_count);
2810 logger().info(
"non_regular_boundary_count: \t{}", non_regular_boundary_count);
2811 logger().info(
"undefined_count: \t{}", undefined_count);
2816 const nlohmann::json &args,
2817 const int n_bases,
const int n_pressure_bases,
2818 const Eigen::MatrixXd &sol,
2820 const Eigen::VectorXi &disc_orders,
2823 const std::string &formulation,
2824 const bool isoparametric,
2825 const int sol_at_node_id,
2831 j[
"geom_order"] = mesh.
orders().size() > 0 ? mesh.
orders().maxCoeff() : 1;
2832 j[
"geom_order_min"] = mesh.
orders().size() > 0 ? mesh.
orders().minCoeff() : 1;
2833 j[
"discr_order_min"] = disc_orders.minCoeff();
2834 j[
"discr_order_max"] = disc_orders.maxCoeff();
2835 j[
"iso_parametric"] = isoparametric;
2836 j[
"problem"] = problem.
name();
2837 j[
"mat_size"] = mat_size;
2838 j[
"num_bases"] = n_bases;
2839 j[
"num_pressure_bases"] = n_pressure_bases;
2840 j[
"num_non_zero"] = nn_zero;
2841 j[
"num_flipped"] = n_flipped;
2842 j[
"num_dofs"] = num_dofs;
2846 j[
"num_p1"] = (disc_orders.array() == 1).count();
2847 j[
"num_p2"] = (disc_orders.array() == 2).count();
2848 j[
"num_p3"] = (disc_orders.array() == 3).count();
2849 j[
"num_p4"] = (disc_orders.array() == 4).count();
2850 j[
"num_p5"] = (disc_orders.array() == 5).count();
2852 j[
"mesh_size"] = mesh_size;
2853 j[
"max_angle"] = max_angle;
2855 j[
"sigma_max"] = sigma_max;
2856 j[
"sigma_min"] = sigma_min;
2857 j[
"sigma_avg"] = sigma_avg;
2859 j[
"min_edge_length"] = min_edge_length;
2860 j[
"average_edge_length"] = average_edge_length;
2862 j[
"err_l2"] = l2_err;
2863 j[
"err_h1"] = h1_err;
2864 j[
"err_h1_semi"] = h1_semi_err;
2865 j[
"err_linf"] = linf_err;
2866 j[
"err_linf_grad"] = grad_max_err;
2867 j[
"err_lp"] = lp_err;
2869 j[
"spectrum"] = {spectrum(0), spectrum(1), spectrum(2), spectrum(3)};
2870 j[
"spectrum_condest"] = std::abs(spectrum(3)) / std::abs(spectrum(0));
2883 j[
"solver_info"] = solver_info;
2885 j[
"count_simplex"] = simplex_count;
2886 j[
"count_regular"] = regular_count;
2887 j[
"count_regular_boundary"] = regular_boundary_count;
2888 j[
"count_simple_singular"] = simple_singular_count;
2889 j[
"count_multi_singular"] = multi_singular_count;
2890 j[
"count_boundary"] = boundary_count;
2891 j[
"count_non_regular_boundary"] = non_regular_boundary_count;
2892 j[
"count_non_regular"] = non_regular_count;
2893 j[
"count_undefined"] = undefined_count;
2894 j[
"count_multi_singular_boundary"] = multi_singular_boundary_count;
2896 j[
"is_simplicial"] = mesh.
n_elements() == simplex_count;
2898 j[
"peak_memory"] =
getPeakRSS() / (1024 * 1024);
2902 std::vector<double> mmin(actual_dim);
2903 std::vector<double> mmax(actual_dim);
2905 for (
int d = 0; d < actual_dim; ++d)
2907 mmin[d] = std::numeric_limits<double>::max();
2908 mmax[d] = -std::numeric_limits<double>::max();
2911 for (
int i = 0; i < sol.size(); i += actual_dim)
2913 for (
int d = 0; d < actual_dim; ++d)
2915 mmin[d] = std::min(mmin[d], sol(i + d));
2916 mmax[d] = std::max(mmax[d], sol(i + d));
2920 std::vector<double> sol_at_node(actual_dim);
2922 if (sol_at_node_id >= 0)
2924 const int node_id = sol_at_node_id;
2926 for (
int d = 0; d < actual_dim; ++d)
2928 sol_at_node[d] = sol(node_id * actual_dim + d);
2932 j[
"sol_at_node"] = sol_at_node;
2933 j[
"sol_min"] = mmin;
2934 j[
"sol_max"] = mmax;
2936#if defined(POLYFEM_WITH_CPP_THREADS)
2938#elif defined(POLYFEM_WITH_TBB)
2941 j[
"num_threads"] = 1;
2944 j[
"formulation"] = formulation;
2950 : file(path), solve_data(solve_data)
2955 file << name <<
",";
2957 file <<
"total_energy" << std::endl;
2974 file << ((form && form->enabled()) ? form->value(sol) : 0) / s <<
",";
2981 : file(path), state(state), t0(t0), dt(dt)
2983 file <<
"step,time,forward,remeshing,global_relaxation,peak_mem,#V,#T" << std::endl;
3007 const double peak_mem =
getPeakRSS() / double(1 << 30);
3010 file << fmt::format(
3011 "{},{},{},{},{},{},{},{}\n",
3012 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
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
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< 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