45#include <paraviewo/VTMWriter.hpp>
46#include <paraviewo/PVDWriter.hpp>
48#include <SimpleBVH/BVH.hpp>
50#include <igl/write_triangle_mesh.h>
52#include <igl/facet_adjacency_matrix.h>
53#include <igl/connected_components.h>
63 void compute_traction_forces(
const State &state,
const Eigen::MatrixXd &solution,
const double t, Eigen::MatrixXd &traction_forces,
bool skip_dirichlet =
true)
66 if (!state.problem->is_scalar())
67 actual_dim = state.mesh->dimension();
71 const std::vector<basis::ElementBases> &bases = state.bases;
72 const std::vector<basis::ElementBases> &gbases = state.geom_bases();
74 Eigen::MatrixXd uv, samples, gtmp, rhs_fun, deform_mat, trafo;
75 Eigen::VectorXi global_primitive_ids;
76 Eigen::MatrixXd
points, normals;
77 Eigen::VectorXd weights;
80 traction_forces.setZero(state.n_bases * actual_dim, 1);
82 for (
const auto &lb : state.total_local_boundary)
90 const basis::ElementBases &gbs = gbases[e];
91 const basis::ElementBases &bs = bases[e];
93 vals.
compute(e, state.mesh->is_volume(), points, bs, gbs);
95 for (
int n = 0; n < normals.rows(); ++n)
99 if (solution.size() > 0)
101 assert(actual_dim == 2 || actual_dim == 3);
102 deform_mat.resize(actual_dim, actual_dim);
103 deform_mat.setZero();
104 for (
const auto &b :
vals.basis_values)
106 for (
const auto &g : b.global)
108 for (
int d = 0; d < actual_dim; ++d)
110 deform_mat.row(d) += solution(g.index * actual_dim + d) * b.grad.row(n);
118 normals.row(n) = normals.row(n) * trafo.inverse();
119 normals.row(n).normalize();
122 std::vector<assembler::Assembler::NamedMatrix> tensor_flat;
123 state.assembler->compute_tensor_value(assembler::OutputData(t, e, bs, gbs, points, solution), tensor_flat);
129 const int g_index = v.
global[0].index * actual_dim;
131 for (
int q = 0; q <
points.rows(); ++q)
134 assert(tensor_flat[0].first ==
"cauchy_stess");
135 assert(tensor_flat[0].second.row(q).size() == actual_dim * actual_dim);
137 Eigen::MatrixXd stress_tensor =
utils::unflatten(tensor_flat[0].second.row(q), actual_dim);
139 traction_forces.block(g_index, 0, actual_dim, 1) += stress_tensor * normals.row(q).transpose() * v.
val(q) * weights(q);
149 const std::vector<basis::ElementBases> &bases,
150 const std::vector<mesh::LocalBoundary> &total_local_boundary,
151 Eigen::MatrixXd &node_positions,
152 Eigen::MatrixXi &boundary_edges,
153 Eigen::MatrixXi &boundary_triangles,
154 std::vector<Eigen::Triplet<double>> &displacement_map_entries)
158 displacement_map_entries.clear();
164 node_positions.resize(n_bases + (is_simplicial ? 0 : mesh.
n_faces()), 3);
165 node_positions.setZero();
166 const Mesh3D &mesh3d =
dynamic_cast<const Mesh3D &
>(mesh);
168 std::vector<std::tuple<int, int, int>> tris;
170 std::vector<bool> visited_node(n_bases,
false);
172 std::stringstream print_warning;
178 for (
int j = 0; j < lb.size(); ++j)
180 const int eid = lb.global_primitive_id(j);
181 const int lid = lb[j];
184 if (mesh.
is_cube(lb.element_id()))
186 assert(!is_simplicial);
188 std::vector<int> loc_nodes;
191 for (
long n = 0; n < nodes.size(); ++n)
193 auto &bs = b.
bases[nodes(n)];
194 const auto &glob = bs.global();
195 if (glob.size() != 1)
198 int gindex = glob.front().index;
199 node_positions.row(gindex) = glob.front().node;
200 bary += glob.front().node;
201 loc_nodes.push_back(gindex);
204 if (loc_nodes.size() != 4)
206 logger().trace(
"skipping element {} since it is not Q1", eid);
212 const int new_node = n_bases + eid;
213 node_positions.row(new_node) = bary;
214 tris.emplace_back(loc_nodes[1], loc_nodes[0], new_node);
215 tris.emplace_back(loc_nodes[2], loc_nodes[1], new_node);
216 tris.emplace_back(loc_nodes[3], loc_nodes[2], new_node);
217 tris.emplace_back(loc_nodes[0], loc_nodes[3], new_node);
219 for (
int q = 0; q < 4; ++q)
221 if (!visited_node[loc_nodes[q]])
222 displacement_map_entries.emplace_back(loc_nodes[q], loc_nodes[q], 1);
224 visited_node[loc_nodes[q]] =
true;
225 displacement_map_entries.emplace_back(new_node, loc_nodes[q], 0.25);
233 logger().trace(
"skipping element {} since it is not a simplex or hex", eid);
239 std::vector<int> loc_nodes;
241 bool is_follower =
false;
244 for (
long n = 0; n < nodes.size(); ++n)
246 auto &bs = b.
bases[nodes(n)];
247 const auto &glob = bs.global();
248 if (glob.size() != 1)
259 for (
long n = 0; n < nodes.size(); ++n)
262 const std::vector<basis::Local2Global> &glob = bs.
global();
263 if (glob.size() != 1)
266 int gindex = glob.front().index;
267 node_positions.row(gindex) = glob.front().node;
268 loc_nodes.push_back(gindex);
271 if (loc_nodes.size() == 3)
273 tris.emplace_back(loc_nodes[0], loc_nodes[1], loc_nodes[2]);
275 else if (loc_nodes.size() == 6)
277 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[5]);
278 tris.emplace_back(loc_nodes[3], loc_nodes[1], loc_nodes[4]);
279 tris.emplace_back(loc_nodes[4], loc_nodes[2], loc_nodes[5]);
280 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[5]);
282 else if (loc_nodes.size() == 10)
284 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[8]);
285 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[9]);
286 tris.emplace_back(loc_nodes[4], loc_nodes[1], loc_nodes[5]);
287 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[9]);
288 tris.emplace_back(loc_nodes[6], loc_nodes[2], loc_nodes[7]);
289 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[9]);
290 tris.emplace_back(loc_nodes[8], loc_nodes[3], loc_nodes[9]);
291 tris.emplace_back(loc_nodes[9], loc_nodes[4], loc_nodes[5]);
292 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[9]);
294 else if (loc_nodes.size() == 15)
296 tris.emplace_back(loc_nodes[0], loc_nodes[3], loc_nodes[11]);
297 tris.emplace_back(loc_nodes[3], loc_nodes[4], loc_nodes[12]);
298 tris.emplace_back(loc_nodes[3], loc_nodes[12], loc_nodes[11]);
299 tris.emplace_back(loc_nodes[12], loc_nodes[10], loc_nodes[11]);
300 tris.emplace_back(loc_nodes[4], loc_nodes[5], loc_nodes[13]);
301 tris.emplace_back(loc_nodes[4], loc_nodes[13], loc_nodes[12]);
302 tris.emplace_back(loc_nodes[12], loc_nodes[13], loc_nodes[14]);
303 tris.emplace_back(loc_nodes[12], loc_nodes[14], loc_nodes[10]);
304 tris.emplace_back(loc_nodes[14], loc_nodes[9], loc_nodes[10]);
305 tris.emplace_back(loc_nodes[5], loc_nodes[1], loc_nodes[6]);
306 tris.emplace_back(loc_nodes[5], loc_nodes[6], loc_nodes[13]);
307 tris.emplace_back(loc_nodes[6], loc_nodes[7], loc_nodes[13]);
308 tris.emplace_back(loc_nodes[13], loc_nodes[7], loc_nodes[14]);
309 tris.emplace_back(loc_nodes[7], loc_nodes[8], loc_nodes[14]);
310 tris.emplace_back(loc_nodes[14], loc_nodes[8], loc_nodes[9]);
311 tris.emplace_back(loc_nodes[8], loc_nodes[2], loc_nodes[9]);
315 print_warning << loc_nodes.size() <<
" ";
321 for (
int k = 0; k < loc_nodes.size(); ++k)
323 if (!visited_node[loc_nodes[k]])
324 displacement_map_entries.emplace_back(loc_nodes[k], loc_nodes[k], 1);
326 visited_node[loc_nodes[k]] =
true;
332 if (print_warning.str().size() > 0)
333 logger().warn(
"Skipping faces as theys have {} nodes, boundary export supported up to p4", print_warning.str());
335 boundary_triangles.resize(tris.size(), 3);
336 for (
int i = 0; i < tris.size(); ++i)
338 boundary_triangles.row(i) << std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]);
341 if (boundary_triangles.rows() > 0)
343 igl::edges(boundary_triangles, boundary_edges);
348 node_positions.resize(n_bases, 2);
349 node_positions.setZero();
350 const Mesh2D &mesh2d =
dynamic_cast<const Mesh2D &
>(mesh);
352 std::vector<std::pair<int, int>> edges;
358 for (
int j = 0; j < lb.size(); ++j)
360 const int eid = lb.global_primitive_id(j);
361 const int lid = lb[j];
366 for (
long n = 0; n < nodes.size(); ++n)
369 const std::vector<basis::Local2Global> &glob = bs.
global();
370 if (glob.size() != 1)
373 int gindex = glob.front().index;
374 node_positions.row(gindex) = glob.front().node.head<2>();
377 edges.emplace_back(prev_node, gindex);
383 boundary_triangles.resize(0, 0);
384 boundary_edges.resize(edges.size(), 2);
385 for (
int i = 0; i < edges.size(); ++i)
387 boundary_edges.row(i) << edges[i].first, edges[i].second;
394 const std::vector<basis::ElementBases> &bases,
395 const std::vector<basis::ElementBases> &gbases,
396 const std::vector<mesh::LocalBoundary> &total_local_boundary,
397 const Eigen::MatrixXd &solution,
398 const int problem_dim,
399 Eigen::MatrixXd &boundary_vis_vertices,
400 Eigen::MatrixXd &boundary_vis_local_vertices,
401 Eigen::MatrixXi &boundary_vis_elements,
402 Eigen::MatrixXi &boundary_vis_elements_ids,
403 Eigen::MatrixXi &boundary_vis_primitive_ids,
404 Eigen::MatrixXd &boundary_vis_normals,
405 Eigen::MatrixXd &displaced_boundary_vis_normals)
const
409 std::vector<Eigen::MatrixXd> lv, vertices, allnormals, displaced_allnormals;
410 std::vector<int> el_ids, global_primitive_ids;
411 Eigen::MatrixXd uv, local_pts, tmp_n, normals, displaced_normals, trafo, deform_mat;
417 std::vector<std::pair<int, int>> edges;
418 std::vector<std::tuple<int, int, int>> tris;
420 for (
auto it = total_local_boundary.begin(); it != total_local_boundary.end(); ++it)
422 const auto &lb = *it;
423 const auto &gbs = gbases[lb.element_id()];
424 const auto &bs = bases[lb.element_id()];
426 for (
int k = 0; k < lb.size(); ++k)
430 case BoundaryType::TRI_LINE:
434 case BoundaryType::QUAD_LINE:
438 case BoundaryType::QUAD:
442 case BoundaryType::TRI:
446 case BoundaryType::POLYGON:
450 case BoundaryType::POLYHEDRON:
453 case BoundaryType::INVALID:
460 vertices.emplace_back();
461 lv.emplace_back(local_pts);
462 el_ids.push_back(lb.element_id());
463 global_primitive_ids.push_back(lb.global_primitive_id(k));
464 gbs.eval_geom_mapping(local_pts, vertices.back());
465 vals.compute(lb.element_id(), mesh.
is_volume(), local_pts, bs, gbs);
466 const int tris_start = tris.size();
470 if (lb.type() == BoundaryType::QUAD)
472 const auto map = [n_samples, size](
int i,
int j) {
return j * n_samples + i + size; };
474 for (
int j = 0; j < n_samples - 1; ++j)
476 for (
int i = 0; i < n_samples - 1; ++i)
478 tris.emplace_back(map(i, j), map(i + 1, j), map(i, j + 1));
479 tris.emplace_back(map(i + 1, j + 1), map(i, j + 1), map(i + 1, j));
483 else if (lb.type() == BoundaryType::TRI)
486 std::vector<int> mapp(n_samples * n_samples, -1);
487 for (
int j = 0; j < n_samples; ++j)
489 for (
int i = 0; i < n_samples - j; ++i)
491 mapp[j * n_samples + i] = index;
495 const auto map = [mapp, n_samples](
int i,
int j) {
496 if (j * n_samples + i >= mapp.size())
498 return mapp[j * n_samples + i];
501 for (
int j = 0; j < n_samples - 1; ++j)
503 for (
int i = 0; i < n_samples - j; ++i)
505 if (map(i, j) >= 0 && map(i + 1, j) >= 0 && map(i, j + 1) >= 0)
506 tris.emplace_back(map(i, j) + size, map(i + 1, j) + size, map(i, j + 1) + size);
508 if (map(i + 1, j + 1) >= 0 && map(i, j + 1) >= 0 && map(i + 1, j) >= 0)
509 tris.emplace_back(map(i + 1, j + 1) + size, map(i, j + 1) + size, map(i + 1, j) + size);
520 for (
int i = 0; i < vertices.back().rows() - 1; ++i)
521 edges.emplace_back(i + size, i + size + 1);
524 normals.resize(
vals.jac_it.size(), tmp_n.cols());
525 displaced_normals.resize(
vals.jac_it.size(), tmp_n.cols());
527 for (
int n = 0; n <
vals.jac_it.size(); ++n)
529 trafo =
vals.jac_it[n].inverse();
531 if (problem_dim == 2 || problem_dim == 3)
534 if (solution.size() > 0)
536 deform_mat.resize(problem_dim, problem_dim);
537 deform_mat.setZero();
538 for (
const auto &b :
vals.basis_values)
539 for (
const auto &g : b.global)
540 for (
int d = 0; d < problem_dim; ++d)
541 deform_mat.row(d) += solution(g.index * problem_dim + d) * b.grad.row(n);
547 normals.row(n) = tmp_n *
vals.jac_it[n];
548 normals.row(n).normalize();
550 displaced_normals.row(n) = tmp_n * trafo.inverse();
551 displaced_normals.row(n).normalize();
554 allnormals.push_back(normals);
555 displaced_allnormals.push_back(displaced_normals);
558 for (
int n = 0; n <
vals.jac_it.size(); ++n)
560 tmp_n += normals.row(n);
565 Eigen::Vector3d e1 = vertices.back().row(std::get<1>(tris.back()) - size) - vertices.back().row(std::get<0>(tris.back()) - size);
566 Eigen::Vector3d e2 = vertices.back().row(std::get<2>(tris.back()) - size) - vertices.back().row(std::get<0>(tris.back()) - size);
568 Eigen::Vector3d n = e1.cross(e2);
569 Eigen::Vector3d nn = tmp_n.transpose();
573 for (
int i = tris_start; i < tris.size(); ++i)
575 tris[i] = std::tuple<int, int, int>(std::get<0>(tris[i]), std::get<2>(tris[i]), std::get<1>(tris[i]));
580 size += vertices.back().rows();
584 boundary_vis_vertices.resize(size, vertices.front().cols());
585 boundary_vis_local_vertices.resize(size, vertices.front().cols());
586 boundary_vis_elements_ids.resize(size, 1);
587 boundary_vis_primitive_ids.resize(size, 1);
588 boundary_vis_normals.resize(size, vertices.front().cols());
589 displaced_boundary_vis_normals.resize(size, vertices.front().cols());
592 boundary_vis_elements.resize(tris.size(), 3);
594 boundary_vis_elements.resize(edges.size(), 2);
598 for (
const auto &v : vertices)
600 boundary_vis_vertices.block(index, 0, v.rows(), v.cols()) = v;
601 boundary_vis_local_vertices.block(index, 0, v.rows(), v.cols()) = lv[ii];
602 boundary_vis_elements_ids.block(index, 0, v.rows(), 1).setConstant(el_ids[ii]);
603 boundary_vis_primitive_ids.block(index, 0, v.rows(), 1).setConstant(global_primitive_ids[ii++]);
608 for (
const auto &n : allnormals)
610 boundary_vis_normals.block(index, 0, n.rows(), n.cols()) = n;
615 for (
const auto &n : displaced_allnormals)
617 displaced_boundary_vis_normals.block(index, 0, n.rows(), n.cols()) = n;
624 for (
const auto &t : tris)
626 boundary_vis_elements.row(index) << std::get<0>(t), std::get<1>(t), std::get<2>(t);
632 for (
const auto &e : edges)
634 boundary_vis_elements.row(index) << e.first, e.second;
642 const Eigen::VectorXi &disc_orders,
643 const std::vector<basis::ElementBases> &gbases,
644 const std::map<int, Eigen::MatrixXd> &polys,
645 const std::map<
int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d,
646 const bool boundary_only,
647 Eigen::MatrixXd &points,
648 Eigen::MatrixXi &tets,
649 Eigen::MatrixXi &el_id,
650 Eigen::MatrixXd &discr)
const
665 const auto ¤t_bases = gbases;
666 int tet_total_size = 0;
667 int pts_total_size = 0;
669 Eigen::MatrixXd vis_pts_poly;
670 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
672 for (
size_t i = 0; i < current_bases.size(); ++i)
674 const auto &bs = current_bases[i];
682 pts_total_size += sampler.simplex_points().rows();
686 tet_total_size += sampler.cube_volume().rows();
687 pts_total_size += sampler.cube_points().rows();
693 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
695 tet_total_size += vis_faces_poly.rows();
696 pts_total_size += vis_pts_poly.rows();
700 sampler.sample_polygon(polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
702 tet_total_size += vis_faces_poly.rows();
703 pts_total_size += vis_pts_poly.rows();
708 points.resize(pts_total_size, mesh.
dimension());
709 tets.resize(tet_total_size, mesh.
is_volume() ? 4 : 3);
711 el_id.resize(pts_total_size, 1);
712 discr.resize(pts_total_size, 1);
714 Eigen::MatrixXd mapped, tmp;
715 int tet_index = 0, pts_index = 0;
717 for (
size_t i = 0; i < current_bases.size(); ++i)
719 const auto &bs = current_bases[i];
726 bs.eval_geom_mapping(sampler.simplex_points(), mapped);
728 tets.block(tet_index, 0, sampler.simplex_volume().rows(), tets.cols()) = sampler.simplex_volume().array() + pts_index;
729 tet_index += sampler.simplex_volume().rows();
731 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
732 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
733 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
734 pts_index += mapped.rows();
738 bs.eval_geom_mapping(sampler.cube_points(), mapped);
740 tets.block(tet_index, 0, sampler.cube_volume().rows(), tets.cols()) = sampler.cube_volume().array() + pts_index;
741 tet_index += sampler.cube_volume().rows();
743 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
744 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(disc_orders(i));
745 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
746 pts_index += mapped.rows();
752 sampler.sample_polyhedron(polys_3d.at(i).first, polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
753 bs.eval_geom_mapping(vis_pts_poly, mapped);
755 tets.block(tet_index, 0, vis_faces_poly.rows(), tets.cols()) = vis_faces_poly.array() + pts_index;
756 tet_index += vis_faces_poly.rows();
758 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
759 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(-1);
760 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
761 pts_index += mapped.rows();
765 sampler.sample_polygon(polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
766 bs.eval_geom_mapping(vis_pts_poly, mapped);
768 tets.block(tet_index, 0, vis_faces_poly.rows(), tets.cols()) = vis_faces_poly.array() + pts_index;
769 tet_index += vis_faces_poly.rows();
771 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
772 discr.block(pts_index, 0, mapped.rows(), 1).setConstant(-1);
773 el_id.block(pts_index, 0, mapped.rows(), 1).setConstant(i);
774 pts_index += mapped.rows();
779 assert(pts_index == points.rows());
780 assert(tet_index == tets.rows());
785 const Eigen::VectorXi &disc_orders,
786 const std::vector<basis::ElementBases> &bases,
787 Eigen::MatrixXd &points,
788 std::vector<std::vector<int>> &elements,
789 Eigen::MatrixXi &el_id,
790 Eigen::MatrixXd &discr)
const
804 std::vector<RowVectorNd> nodes;
805 int pts_total_size = 0;
806 elements.resize(bases.size());
807 Eigen::MatrixXd ref_pts;
809 for (
size_t i = 0; i < bases.size(); ++i)
811 const auto &bs = bases[i];
829 const int n_v =
static_cast<const mesh::Mesh2D &
>(mesh).n_face_vertices(i);
830 ref_pts.resize(n_v, 2);
834 pts_total_size += ref_pts.rows();
837 points.resize(pts_total_size, mesh.
dimension());
839 el_id.resize(pts_total_size, 1);
840 discr.resize(pts_total_size, 1);
842 Eigen::MatrixXd mapped;
845 std::string error_msg =
"";
847 for (
size_t i = 0; i < bases.size(); ++i)
849 const auto &bs = bases[i];
869 bs.eval_geom_mapping(ref_pts, mapped);
871 for (
int j = 0; j < mapped.rows(); ++j)
873 points.row(pts_index) = mapped.row(j);
874 el_id(pts_index) = i;
875 discr(pts_index) = disc_orders(i);
876 elements[i].push_back(pts_index);
885 const int n_nodes = elements[i].size();
886 if (disc_orders(i) >= 3)
888 std::swap(elements[i][16], elements[i][17]);
889 std::swap(elements[i][17], elements[i][18]);
890 std::swap(elements[i][18], elements[i][19]);
892 if (disc_orders(i) > 4)
893 error_msg =
"Saving high-order meshes not implemented for P5+ elements!";
897 if (disc_orders(i) == 4)
899 const int n_nodes = elements[i].size();
900 std::swap(elements[i][n_nodes - 1], elements[i][n_nodes - 2]);
902 if (disc_orders(i) > 4)
903 error_msg =
"Saving high-order meshes not implemented for P5+ elements!";
906 else if (disc_orders(i) > 1)
907 error_msg =
"Saving high-order meshes not implemented for Q2+ elements!";
910 if (!error_msg.empty())
913 for (
size_t i = 0; i < bases.size(); ++i)
918 const auto &mesh2d =
static_cast<const mesh::Mesh2D &
>(mesh);
921 for (
int j = 0; j < n_v; ++j)
923 points.row(pts_index) = mesh2d.point(mesh2d.face_vertex(i, j));
924 el_id(pts_index) = i;
925 discr(pts_index) = disc_orders(i);
926 elements[i].push_back(pts_index);
932 assert(pts_index == points.rows());
937 const Eigen::MatrixXd &sol,
938 const Eigen::MatrixXd &pressure,
939 const bool is_time_dependent,
940 const double tend_in,
943 const std::string &vis_mesh_path,
944 const std::string &nodes_path,
945 const std::string &solution_path,
946 const std::string &stress_path,
947 const std::string &mises_path,
948 const bool is_contact_enabled,
949 std::vector<SolutionFrame> &solution_frames)
const
953 logger().error(
"Load the mesh first!");
956 const int n_bases = state.
n_bases;
957 const std::vector<basis::ElementBases> &bases = state.
bases;
958 const std::vector<basis::ElementBases> &gbases = state.
geom_bases();
961 const Eigen::MatrixXd &rhs = state.
rhs;
966 logger().error(
"Build the bases first!");
976 logger().error(
"Solve the problem first!");
980 if (!solution_path.empty())
982 std::ofstream out(solution_path);
984 out << std::scientific;
988 Eigen::VectorXi reordering(n_bases);
989 reordering.setConstant(-1);
991 for (
int i = 0; i < in_node_to_node.size(); ++i)
993 reordering[in_node_to_node[i]] = i;
996 Eigen::MatrixXd tmp(tmp_sol.rows(), tmp_sol.cols());
998 for (
int i = 0; i < reordering.size(); ++i)
1000 if (reordering[i] < 0)
1003 tmp.row(reordering[i]) = tmp_sol.row(i);
1006 for (
int i = 0; i < tmp.rows(); ++i)
1008 for (
int j = 0; j < tmp.cols(); ++j)
1009 out << tmp(i, j) <<
" ";
1015 out << sol << std::endl;
1019 double tend = tend_in;
1023 if (!vis_mesh_path.empty() && !is_time_dependent)
1026 vis_mesh_path, state, sol, pressure,
1028 is_contact_enabled, solution_frames);
1030 if (!nodes_path.empty())
1032 Eigen::MatrixXd nodes(n_bases, mesh.
dimension());
1038 for (
size_t ii = 0; ii < b.global().size(); ++ii)
1040 const auto &lg = b.global()[ii];
1041 nodes.row(lg.index) = lg.node;
1045 std::ofstream out(nodes_path);
1050 if (!stress_path.empty())
1052 Eigen::MatrixXd result;
1053 Eigen::VectorXd mises;
1057 sol, tend, result, mises);
1058 std::ofstream out(stress_path);
1062 if (!mises_path.empty())
1064 Eigen::MatrixXd result;
1065 Eigen::VectorXd mises;
1069 sol, tend, result, mises);
1070 std::ofstream out(mises_path);
1078 volume = args[
"output"][
"paraview"][
"volume"];
1079 surface = args[
"output"][
"paraview"][
"surface"];
1080 wire = args[
"output"][
"paraview"][
"wireframe"];
1081 points = args[
"output"][
"paraview"][
"points"];
1082 contact_forces = args[
"output"][
"paraview"][
"options"][
"contact_forces"] && !is_problem_scalar;
1083 friction_forces = args[
"output"][
"paraview"][
"options"][
"friction_forces"] && !is_problem_scalar;
1088 body_ids = args[
"output"][
"paraview"][
"options"][
"body_ids"];
1089 sol_on_grid = args[
"output"][
"advanced"][
"sol_on_grid"] > 0;
1090 velocity = args[
"output"][
"paraview"][
"options"][
"velocity"];
1091 acceleration = args[
"output"][
"paraview"][
"options"][
"acceleration"];
1092 forces = args[
"output"][
"paraview"][
"options"][
"forces"] && !is_problem_scalar;
1094 scalar_values = args[
"output"][
"paraview"][
"options"][
"scalar_values"];
1095 tensor_values = args[
"output"][
"paraview"][
"options"][
"tensor_values"] && !is_problem_scalar;
1096 discretization_order = args[
"output"][
"paraview"][
"options"][
"discretization_order"] && !is_problem_scalar;
1097 nodes = args[
"output"][
"paraview"][
"options"][
"nodes"] && !is_problem_scalar;
1099 use_spline = args[
"space"][
"basis_type"] ==
"Spline";
1101 reorder_output = args[
"output"][
"data"][
"advanced"][
"reorder_nodes"];
1103 use_hdf5 = args[
"output"][
"paraview"][
"options"][
"use_hdf5"];
1109 const std::string &path,
1111 const Eigen::MatrixXd &sol,
1112 const Eigen::MatrixXd &pressure,
1116 const bool is_contact_enabled,
1117 std::vector<SolutionFrame> &solution_frames)
const
1121 logger().error(
"Load the mesh first!");
1125 const Eigen::MatrixXd &rhs = state.
rhs;
1129 logger().error(
"Build the bases first!");
1137 if (sol.size() <= 0)
1139 logger().error(
"Solve the problem first!");
1143 const std::filesystem::path fs_path(path);
1144 const std::string path_stem = fs_path.stem().string();
1145 const std::string base_path = (fs_path.parent_path() / path_stem).
string();
1155 is_contact_enabled, solution_frames);
1161 is_contact_enabled, solution_frames);
1177 paraviewo::VTMWriter vtm(t);
1179 vtm.add_dataset(
"Volume",
"data", path_stem + opts.
file_extension());
1181 vtm.add_dataset(
"Surface",
"data", path_stem +
"_surf" + opts.
file_extension());
1183 vtm.add_dataset(
"Contact",
"data", path_stem +
"_surf_contact" + opts.
file_extension());
1185 vtm.add_dataset(
"Wireframe",
"data", path_stem +
"_wire" + opts.
file_extension());
1187 vtm.add_dataset(
"Points",
"data", path_stem +
"_points" + opts.
file_extension());
1188 vtm.save(base_path +
".vtm");
1192 const std::string &path,
1194 const Eigen::MatrixXd &sol,
1195 const Eigen::MatrixXd &pressure,
1199 std::vector<SolutionFrame> &solution_frames)
const
1201 const Eigen::VectorXi &disc_orders = state.
disc_orders;
1203 const std::vector<basis::ElementBases> &bases = state.
bases;
1204 const std::vector<basis::ElementBases> &pressure_bases = state.
pressure_bases;
1205 const std::vector<basis::ElementBases> &gbases = state.
geom_bases();
1206 const std::map<int, Eigen::MatrixXd> &polys = state.
polys;
1207 const std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &polys_3d = state.
polys_3d;
1214 Eigen::MatrixXd points;
1215 Eigen::MatrixXi tets;
1216 Eigen::MatrixXi el_id;
1217 Eigen::MatrixXd discr;
1218 std::vector<std::vector<int>> elements;
1223 points, tets, el_id, discr);
1226 points, elements, el_id, discr);
1228 Eigen::MatrixXd fun, exact_fun, err, node_fun;
1233 Eigen::MatrixXd tmp, tmp_grad;
1234 Eigen::MatrixXd tmp_p, tmp_grad_p;
1236 res.setConstant(std::numeric_limits<double>::quiet_NaN());
1238 res_grad.setConstant(std::numeric_limits<double>::quiet_NaN());
1241 res_p.setConstant(std::numeric_limits<double>::quiet_NaN());
1243 res_grad_p.setConstant(std::numeric_limits<double>::quiet_NaN());
1252 Eigen::MatrixXd pt(1, bc.cols() - 1);
1253 for (
int d = 1; d < bc.cols(); ++d)
1256 mesh, problem.
is_scalar(), bases, gbases,
1257 el_id, pt, sol, tmp, tmp_grad);
1260 res_grad.row(i) = tmp_grad;
1265 mesh, 1, pressure_bases, gbases,
1266 el_id, pt, pressure, tmp_p, tmp_grad_p);
1267 res_p.row(i) = tmp_p;
1268 res_grad_p.row(i) = tmp_grad_p;
1272 std::ofstream os(path +
"_sol.txt");
1275 std::ofstream osg(path +
"_grad.txt");
1278 std::ofstream osgg(path +
"_grid.txt");
1283 std::ofstream osp(path +
"_p_sol.txt");
1286 std::ofstream osgp(path +
"_p_grad.txt");
1297 Eigen::MatrixXd tmp = Eigen::VectorXd::LinSpaced(sol.size(), 0, sol.size() - 1);
1307 fun.conservativeResize(fun.rows() + obstacle.
n_vertices(), fun.cols());
1308 node_fun.conservativeResize(node_fun.rows() + obstacle.
n_vertices(), node_fun.cols());
1309 node_fun.bottomRows(obstacle.
n_vertices()).setZero();
1317 problem.
exact(points, t, exact_fun);
1318 err = (fun - exact_fun).eval().rowwise().norm();
1322 exact_fun.conservativeResize(exact_fun.rows() + obstacle.
n_vertices(), exact_fun.cols());
1326 err.conservativeResize(err.rows() + obstacle.
n_vertices(), 1);
1327 err.bottomRows(obstacle.
n_vertices()).setZero();
1331 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
1333 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
1335 tmpw = std::make_shared<paraviewo::VTUWriter>();
1336 paraviewo::ParaviewWriter &writer = *tmpw;
1339 writer.add_field(
"nodes", node_fun);
1343 bool is_time_integrator_valid = time_integrator !=
nullptr;
1347 const Eigen::VectorXd velocity =
1348 is_time_integrator_valid ? (time_integrator->v_prev()) : Eigen::VectorXd::Zero(sol.size());
1354 const Eigen::VectorXd acceleration =
1355 is_time_integrator_valid ? (time_integrator->a_prev()) : Eigen::VectorXd::Zero(sol.size());
1369 if (form ==
nullptr)
1372 Eigen::VectorXd force;
1373 if (form->enabled())
1375 form->first_derivative(sol, force);
1380 force.setZero(sol.size());
1390 Eigen::MatrixXd interp_p;
1398 interp_p.conservativeResize(interp_p.size() + obstacle.
n_vertices(), 1);
1399 interp_p.bottomRows(obstacle.
n_vertices()).setZero();
1403 writer.add_field(
"pressure", interp_p);
1405 solution_frames.back().pressure = interp_p;
1410 discr.conservativeResize(discr.size() + obstacle.
n_vertices(), 1);
1411 discr.bottomRows(obstacle.
n_vertices()).setZero();
1415 writer.add_field(
"discr", discr);
1421 writer.add_field(
"exact", exact_fun);
1422 writer.add_field(
"error", err);
1426 solution_frames.back().exact = exact_fun;
1427 solution_frames.back().error = err;
1431 if (fun.cols() != 1)
1433 std::vector<assembler::Assembler::NamedMatrix>
vals, tvals;
1435 mesh, problem.
is_scalar(), bases, gbases,
1440 for (
auto &[_, v] :
vals)
1447 for (
const auto &[name, v] :
vals)
1448 writer.add_field(name, v);
1450 else if (
vals.size() > 0)
1451 solution_frames.back().scalar_value =
vals[0].second;
1461 for (
auto &[_, v] : tvals)
1464 for (
const auto &[name, v] : tvals)
1467 assert(v.cols() % stride == 0);
1469 for (
int i = 0; i < v.cols(); i += stride)
1471 const Eigen::MatrixXd tmp = v.middleCols(i, stride);
1472 assert(tmp.cols() == stride);
1474 const int ii = (i / stride) + 1;
1475 writer.add_field(fmt::format(
"{:s}_{:d}", name, ii), tmp);
1490 for (
auto &v :
vals)
1492 v.second.conservativeResize(v.second.size() + obstacle.
n_vertices(), 1);
1493 v.second.bottomRows(obstacle.
n_vertices()).setZero();
1501 for (
const auto &v :
vals)
1502 writer.add_field(fmt::format(
"{:s}_avg", v.first), v.second);
1504 else if (
vals.size() > 0)
1505 solution_frames.back().scalar_value_avg =
vals[0].second;
1519 std::map<std::string, Eigen::MatrixXd> param_val;
1520 for (
const auto &[p, _] : params)
1521 param_val[p] = Eigen::MatrixXd(points.rows(), 1);
1522 Eigen::MatrixXd rhos(points.rows(), 1);
1524 Eigen::MatrixXd local_pts;
1525 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
1529 for (
int e = 0; e < int(bases.size()); ++e)
1537 local_pts = sampler.simplex_points();
1539 local_pts = sampler.cube_points();
1543 sampler.sample_polyhedron(polys_3d.at(e).first, polys_3d.at(e).second, local_pts, vis_faces_poly, vis_edges_poly);
1545 sampler.sample_polygon(polys.at(e), local_pts, vis_faces_poly, vis_edges_poly);
1567 const auto &mesh2d =
static_cast<const mesh::Mesh2D &
>(mesh);
1569 local_pts.resize(n_v, 2);
1571 for (
int j = 0; j < n_v; ++j)
1573 local_pts.row(j) = mesh2d.point(mesh2d.face_vertex(e, j));
1582 for (
int j = 0; j <
vals.val.rows(); ++j)
1584 for (
const auto &[p, func] : params)
1585 param_val.at(p)(index) = func(local_pts.row(j),
vals.val.row(j), t, e);
1587 rhos(index) = density(local_pts.row(j),
vals.val.row(j), t, e);
1593 assert(index == points.rows());
1597 for (
auto &[_, tmp] : param_val)
1599 tmp.conservativeResize(tmp.size() + obstacle.
n_vertices(), 1);
1600 tmp.bottomRows(obstacle.
n_vertices()).setZero();
1603 rhos.conservativeResize(rhos.size() + obstacle.
n_vertices(), 1);
1604 rhos.bottomRows(obstacle.
n_vertices()).setZero();
1606 for (
const auto &[p, tmp] : param_val)
1607 writer.add_field(p, tmp);
1608 writer.add_field(
"rho", rhos);
1614 Eigen::MatrixXd ids(points.rows(), 1);
1616 for (
int i = 0; i < points.rows(); ++i)
1623 ids.conservativeResize(ids.size() + obstacle.
n_vertices(), 1);
1624 ids.bottomRows(obstacle.
n_vertices()).setZero();
1627 writer.add_field(
"body_ids", ids);
1633 if (fun.cols() != 1)
1635 Eigen::MatrixXd traction_forces, traction_forces_fun;
1636 compute_traction_forces(state, sol, t, traction_forces,
false);
1645 traction_forces_fun.conservativeResize(traction_forces_fun.rows() + obstacle.
n_vertices(), traction_forces_fun.cols());
1646 traction_forces_fun.bottomRows(obstacle.
n_vertices()).setZero();
1649 writer.add_field(
"traction_force", traction_forces_fun);
1652 if (fun.cols() != 1)
1656 Eigen::MatrixXd potential_grad, potential_grad_fun;
1666 potential_grad_fun.conservativeResize(potential_grad_fun.rows() + obstacle.
n_vertices(), potential_grad_fun.cols());
1667 potential_grad_fun.bottomRows(obstacle.
n_vertices()).setZero();
1670 writer.add_field(
"gradient_of_potential", potential_grad_fun);
1672 catch (std::exception &)
1679 writer.add_field(
"solution", fun);
1681 solution_frames.back().solution = fun;
1687 const int orig_p = points.rows();
1688 points.conservativeResize(points.rows() + obstacle.
n_vertices(), points.cols());
1689 points.bottomRows(obstacle.
n_vertices()) = obstacle.
v();
1691 if (elements.empty())
1693 for (
int i = 0; i < tets.rows(); ++i)
1695 elements.emplace_back();
1696 for (
int j = 0; j < tets.cols(); ++j)
1697 elements.back().push_back(tets(i, j));
1703 elements.emplace_back();
1710 elements.emplace_back();
1717 elements.emplace_back();
1722 if (elements.empty())
1723 writer.write_mesh(path, points, tets);
1725 writer.write_mesh(path, points, elements,
true, disc_orders.maxCoeff() == 1);
1729 solution_frames.back().name = path;
1730 solution_frames.back().points = points;
1731 solution_frames.back().connectivity = tets;
1737 const Eigen::MatrixXd &points,
1739 const std::string &name,
1740 const Eigen::VectorXd &field,
1741 paraviewo::ParaviewWriter &writer)
const
1743 Eigen::MatrixXd inerpolated_field;
1751 inerpolated_field.conservativeResize(
1759 writer.add_field(name, inerpolated_field);
1765 const std::string &export_surface,
1767 const Eigen::MatrixXd &sol,
1768 const Eigen::MatrixXd &pressure,
1772 const bool is_contact_enabled,
1773 std::vector<SolutionFrame> &solution_frames)
const
1776 const Eigen::VectorXi &disc_orders = state.
disc_orders;
1778 const std::vector<basis::ElementBases> &bases = state.
bases;
1779 const std::vector<basis::ElementBases> &pressure_bases = state.
pressure_bases;
1780 const std::vector<basis::ElementBases> &gbases = state.
geom_bases();
1786 Eigen::MatrixXd boundary_vis_vertices;
1787 Eigen::MatrixXd boundary_vis_local_vertices;
1788 Eigen::MatrixXi boundary_vis_elements;
1789 Eigen::MatrixXi boundary_vis_elements_ids;
1790 Eigen::MatrixXi boundary_vis_primitive_ids;
1791 Eigen::MatrixXd boundary_vis_normals;
1792 Eigen::MatrixXd displaced_boundary_vis_normals;
1795 boundary_vis_vertices, boundary_vis_local_vertices, boundary_vis_elements,
1796 boundary_vis_elements_ids, boundary_vis_primitive_ids, boundary_vis_normals,
1797 displaced_boundary_vis_normals);
1799 Eigen::MatrixXd fun, interp_p, discr, vect, b_sidesets;
1801 Eigen::MatrixXd lsol, lp, lgrad, lpgrad;
1807 discr.resize(boundary_vis_vertices.rows(), 1);
1808 fun.resize(boundary_vis_vertices.rows(), actual_dim);
1809 interp_p.resize(boundary_vis_vertices.rows(), 1);
1810 vect.resize(boundary_vis_vertices.rows(), mesh.
dimension());
1812 b_sidesets.resize(boundary_vis_vertices.rows(), 1);
1813 b_sidesets.setZero();
1815 for (
int i = 0; i < boundary_vis_vertices.rows(); ++i)
1817 const auto s_id = mesh.
get_boundary_id(boundary_vis_primitive_ids(i));
1820 b_sidesets(i) = s_id;
1823 const int el_index = boundary_vis_elements_ids(i);
1825 mesh, problem.
is_scalar(), bases, gbases,
1826 el_index, boundary_vis_local_vertices.row(i), sol, lsol, lgrad);
1827 assert(lsol.size() == actual_dim);
1831 mesh, 1, pressure_bases, gbases,
1832 el_index, boundary_vis_local_vertices.row(i), pressure, lp, lpgrad);
1833 assert(lp.size() == 1);
1834 interp_p(i) = lp(0);
1837 discr(i) = disc_orders(el_index);
1838 for (
int j = 0; j < actual_dim; ++j)
1840 fun(i, j) = lsol(j);
1843 if (actual_dim == 1)
1845 assert(lgrad.size() == mesh.
dimension());
1846 for (
int j = 0; j < mesh.
dimension(); ++j)
1848 vect(i, j) = lgrad(j);
1853 assert(lgrad.size() == actual_dim * actual_dim);
1854 std::vector<assembler::Assembler::NamedMatrix> tensor_flat;
1859 assert(tensor_flat[0].first ==
"cauchy_stess");
1860 assert(tensor_flat[0].second.size() == actual_dim * actual_dim);
1862 Eigen::Map<Eigen::MatrixXd> tensor(tensor_flat[0].second.data(), actual_dim, actual_dim);
1863 vect.row(i) = displaced_boundary_vis_normals.row(i) * tensor;
1869 area = mesh.
tri_area(boundary_vis_primitive_ids(i));
1870 else if (mesh.
is_cube(el_index))
1871 area = mesh.
quad_area(boundary_vis_primitive_ids(i));
1874 area = mesh.
edge_length(boundary_vis_primitive_ids(i));
1876 vect.row(i) *= area;
1880 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
1882 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
1884 tmpw = std::make_shared<paraviewo::VTUWriter>();
1885 paraviewo::ParaviewWriter &writer = *tmpw;
1890 writer.add_field(
"normals", boundary_vis_normals);
1891 writer.add_field(
"displaced_normals", displaced_boundary_vis_normals);
1893 writer.add_field(
"pressure", interp_p);
1894 writer.add_field(
"discr", discr);
1895 writer.add_field(
"sidesets", b_sidesets);
1897 if (actual_dim == 1)
1898 writer.add_field(
"solution_grad", vect);
1901 writer.add_field(
"traction_force", vect);
1907 solution_frames.back().pressure = interp_p;
1914 std::map<std::string, Eigen::MatrixXd> param_val;
1915 for (
const auto &[p, _] : params)
1916 param_val[p] = Eigen::MatrixXd(boundary_vis_vertices.rows(), 1);
1917 Eigen::MatrixXd rhos(boundary_vis_vertices.rows(), 1);
1919 for (
int i = 0; i < boundary_vis_vertices.rows(); ++i)
1923 for (
const auto &[p, func] : params)
1924 param_val.at(p)(i) = func(boundary_vis_local_vertices.row(i), boundary_vis_vertices.row(i), t, boundary_vis_elements_ids(i));
1926 rhos(i) = density(boundary_vis_local_vertices.row(i), boundary_vis_vertices.row(i), t, boundary_vis_elements_ids(i));
1929 for (
const auto &[p, tmp] : param_val)
1930 writer.add_field(p, tmp);
1931 writer.add_field(
"rho", rhos);
1937 Eigen::MatrixXd ids(boundary_vis_vertices.rows(), 1);
1939 for (
int i = 0; i < boundary_vis_vertices.rows(); ++i)
1941 ids(i) = mesh.
get_body_id(boundary_vis_elements_ids(i));
1944 writer.add_field(
"body_ids", ids);
1949 writer.add_field(
"solution", fun);
1951 solution_frames.back().solution = fun;
1954 writer.write_mesh(export_surface, boundary_vis_vertices, boundary_vis_elements);
1957 solution_frames.back().name = export_surface;
1958 solution_frames.back().points = boundary_vis_vertices;
1959 solution_frames.back().connectivity = boundary_vis_elements;
1964 const std::string &export_surface,
1966 const Eigen::MatrixXd &sol,
1967 const Eigen::MatrixXd &pressure,
1971 const bool is_contact_enabled,
1972 std::vector<SolutionFrame> &solution_frames)
const
1976 const double dhat = state.
args[
"contact"][
"dhat"];
1977 const double friction_coefficient = state.
args[
"contact"][
"friction_coefficient"];
1978 const double epsv = state.
args[
"contact"][
"epsv"];
1984 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
1986 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
1988 tmpw = std::make_shared<paraviewo::VTUWriter>();
1989 paraviewo::ParaviewWriter &writer = *tmpw;
1991 const int problem_dim = mesh.
dimension();
1992 const Eigen::MatrixXd full_displacements =
utils::unflatten(sol, problem_dim);
1993 const Eigen::MatrixXd surface_displacements = collision_mesh.map_displacements(full_displacements);
1995 const Eigen::MatrixXd displaced_surface = collision_mesh.displace_vertices(full_displacements);
1997 ipc::Collisions collision_set;
1998 collision_set.set_use_convergent_formulation(state.
args[
"contact"][
"use_convergent_formulation"]);
1999 collision_set.build(
2000 collision_mesh, displaced_surface, dhat,
2001 0, state.
args[
"solver"][
"contact"][
"CCD"][
"broad_phase"]);
2003 ipc::BarrierPotential barrier_potential(dhat);
2005 const double barrier_stiffness = contact_form !=
nullptr ? contact_form->barrier_stiffness() : 1;
2009 Eigen::MatrixXd forces = -barrier_stiffness * barrier_potential.gradient(collision_set, collision_mesh, displaced_surface);
2013 assert(forces_reshaped.rows() == surface_displacements.rows());
2014 assert(forces_reshaped.cols() == surface_displacements.cols());
2015 writer.add_field(
"contact_forces", forces_reshaped);
2020 ipc::FrictionCollisions friction_collision_set;
2021 friction_collision_set.build(
2022 collision_mesh, displaced_surface, collision_set,
2023 barrier_potential, barrier_stiffness, friction_coefficient);
2025 ipc::FrictionPotential friction_potential(epsv);
2027 Eigen::MatrixXd velocities;
2032 velocities = collision_mesh.map_displacements(
utils::unflatten(velocities, collision_mesh.dim()));
2034 Eigen::MatrixXd forces = -friction_potential.gradient(
2035 friction_collision_set, collision_mesh, velocities);
2039 assert(forces_reshaped.rows() == surface_displacements.rows());
2040 assert(forces_reshaped.cols() == surface_displacements.cols());
2041 writer.add_field(
"friction_forces", forces_reshaped);
2044 assert(collision_mesh.rest_positions().rows() == surface_displacements.rows());
2045 assert(collision_mesh.rest_positions().cols() == surface_displacements.cols());
2048 writer.add_field(
"solution", surface_displacements);
2051 export_surface.substr(0, export_surface.length() - 4) +
"_contact.vtu",
2052 collision_mesh.rest_positions(),
2053 problem_dim == 3 ? collision_mesh.faces() : collision_mesh.edges());
2058 const std::string &name,
2060 const Eigen::MatrixXd &sol,
2063 std::vector<SolutionFrame> &solution_frames)
const
2065 const std::vector<basis::ElementBases> &gbases = state.
geom_bases();
2073 Eigen::MatrixXi vis_faces_poly, vis_edges_poly;
2074 Eigen::MatrixXd vis_pts_poly;
2076 const auto ¤t_bases = gbases;
2077 int seg_total_size = 0;
2078 int pts_total_size = 0;
2079 int faces_total_size = 0;
2081 for (
size_t i = 0; i < current_bases.size(); ++i)
2083 const auto &bs = current_bases[i];
2088 seg_total_size += sampler.simplex_edges().rows();
2089 faces_total_size += sampler.simplex_faces().rows();
2093 pts_total_size += sampler.cube_points().rows();
2094 seg_total_size += sampler.cube_edges().rows();
2095 faces_total_size += sampler.cube_faces().rows();
2100 sampler.sample_polyhedron(state.
polys_3d.at(i).first, state.
polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
2102 sampler.sample_polygon(state.
polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
2104 pts_total_size += vis_pts_poly.rows();
2105 seg_total_size += vis_edges_poly.rows();
2106 faces_total_size += vis_faces_poly.rows();
2110 Eigen::MatrixXd points(pts_total_size, mesh.
dimension());
2111 Eigen::MatrixXi edges(seg_total_size, 2);
2112 Eigen::MatrixXi
faces(faces_total_size, 3);
2115 Eigen::MatrixXd mapped, tmp;
2116 int seg_index = 0, pts_index = 0, face_index = 0;
2117 for (
size_t i = 0; i < current_bases.size(); ++i)
2119 const auto &bs = current_bases[i];
2123 bs.eval_geom_mapping(sampler.simplex_points(), mapped);
2124 edges.block(seg_index, 0, sampler.simplex_edges().rows(), edges.cols()) = sampler.simplex_edges().array() + pts_index;
2125 seg_index += sampler.simplex_edges().rows();
2127 faces.block(face_index, 0, sampler.simplex_faces().rows(), 3) = sampler.simplex_faces().array() + pts_index;
2128 face_index += sampler.simplex_faces().rows();
2130 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2131 pts_index += mapped.rows();
2135 bs.eval_geom_mapping(sampler.cube_points(), mapped);
2136 edges.block(seg_index, 0, sampler.cube_edges().rows(), edges.cols()) = sampler.cube_edges().array() + pts_index;
2137 seg_index += sampler.cube_edges().rows();
2139 faces.block(face_index, 0, sampler.cube_faces().rows(), 3) = sampler.cube_faces().array() + pts_index;
2140 face_index += sampler.cube_faces().rows();
2142 points.block(pts_index, 0, mapped.rows(), points.cols()) = mapped;
2143 pts_index += mapped.rows();
2148 sampler.sample_polyhedron(state.
polys_3d.at(i).first, state.
polys_3d.at(i).second, vis_pts_poly, vis_faces_poly, vis_edges_poly);
2150 sampler.sample_polygon(state.
polys.at(i), vis_pts_poly, vis_faces_poly, vis_edges_poly);
2152 edges.block(seg_index, 0, vis_edges_poly.rows(), edges.cols()) = vis_edges_poly.array() + pts_index;
2153 seg_index += vis_edges_poly.rows();
2155 faces.block(face_index, 0, vis_faces_poly.rows(), 3) = vis_faces_poly.array() + pts_index;
2156 face_index += vis_faces_poly.rows();
2158 points.block(pts_index, 0, vis_pts_poly.rows(), points.cols()) = vis_pts_poly;
2159 pts_index += vis_pts_poly.rows();
2163 assert(pts_index == points.rows());
2164 assert(face_index ==
faces.rows());
2169 for (
long i = 0; i <
faces.rows(); ++i)
2171 const int v0 =
faces(i, 0);
2172 const int v1 =
faces(i, 1);
2173 const int v2 =
faces(i, 2);
2175 int tmpc =
faces(i, 2);
2182 Eigen::Matrix2d mmat;
2183 for (
long i = 0; i <
faces.rows(); ++i)
2185 const int v0 =
faces(i, 0);
2186 const int v1 =
faces(i, 1);
2187 const int v2 =
faces(i, 2);
2189 mmat.row(0) = points.row(v2) - points.row(v0);
2190 mmat.row(1) = points.row(v1) - points.row(v0);
2192 if (mmat.determinant() > 0)
2194 int tmpc =
faces(i, 2);
2201 Eigen::MatrixXd fun;
2205 pts_index, sol, fun,
true,
false);
2207 Eigen::MatrixXd exact_fun, err;
2211 problem.
exact(points, t, exact_fun);
2212 err = (fun - exact_fun).eval().rowwise().norm();
2215 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2217 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2219 tmpw = std::make_shared<paraviewo::VTUWriter>();
2220 paraviewo::ParaviewWriter &writer = *tmpw;
2224 writer.add_field(
"exact", exact_fun);
2225 writer.add_field(
"error", err);
2228 if (fun.cols() != 1)
2230 std::vector<assembler::Assembler::NamedMatrix> scalar_val;
2236 for (
const auto &v : scalar_val)
2237 writer.add_field(v.first, v.second);
2240 writer.add_field(
"solution", fun);
2242 writer.write_mesh(name, points, edges);
2246 const std::string &path,
2248 const Eigen::MatrixXd &sol,
2250 std::vector<SolutionFrame> &solution_frames)
const
2261 Eigen::MatrixXd fun(dirichlet_nodes_position.size(), actual_dim);
2262 Eigen::MatrixXd b_sidesets(dirichlet_nodes_position.size(), 1);
2263 b_sidesets.setZero();
2264 Eigen::MatrixXd points(dirichlet_nodes_position.size(), mesh.
dimension());
2265 std::vector<std::vector<int>> cells(dirichlet_nodes_position.size());
2267 for (
int i = 0; i < dirichlet_nodes_position.size(); ++i)
2269 const int n_id = dirichlet_nodes[i];
2273 b_sidesets(i) = s_id;
2276 for (
int j = 0; j < actual_dim; ++j)
2278 fun(i, j) = sol(n_id * actual_dim + j);
2281 points.row(i) = dirichlet_nodes_position[i];
2282 cells[i].push_back(i);
2285 std::shared_ptr<paraviewo::ParaviewWriter> tmpw;
2287 tmpw = std::make_shared<paraviewo::HDF5VTUWriter>();
2289 tmpw = std::make_shared<paraviewo::VTUWriter>();
2290 paraviewo::ParaviewWriter &writer = *tmpw;
2294 writer.add_field(
"sidesets", b_sidesets);
2296 writer.add_field(
"solution", fun);
2297 writer.write_mesh(path, points, cells,
false,
false);
2302 const std::string &name,
2303 const std::function<std::string(
int)> &vtu_names,
2304 int time_steps,
double t0,
double dt,
int skip_frame)
const
2306 paraviewo::PVDWriter::save_pvd(name, vtu_names, time_steps, t0, dt, skip_frame);
2322 const int nx = delta[0] / spacing + 1;
2323 const int ny = delta[1] / spacing + 1;
2324 const int nz = delta.cols() >= 3 ? (delta[2] / spacing + 1) : 1;
2325 const int n = nx * ny * nz;
2329 for (
int i = 0; i < nx; ++i)
2331 const double x = (delta[0] / (nx - 1)) * i + min[0];
2333 for (
int j = 0; j < ny; ++j)
2335 const double y = (delta[1] / (ny - 1)) * j + min[1];
2337 if (delta.cols() <= 2)
2343 for (
int k = 0; k < nz; ++k)
2345 const double z = (delta[2] / (nz - 1)) * k + min[2];
2354 std::vector<std::array<Eigen::Vector3d, 2>> boxes;
2360 const double eps = 1e-6;
2369 const Eigen::Vector3d min(
2374 const Eigen::Vector3d max(
2379 std::vector<unsigned int> candidates;
2381 bvh.intersect_box(min, max, candidates);
2383 for (
const auto cand : candidates)
2387 logger().warn(
"Element {} is not simplex, skipping", cand);
2391 Eigen::MatrixXd coords;
2394 for (
int d = 0; d < coords.size(); ++d)
2396 if (fabs(coords(d)) < 1e-8)
2398 else if (fabs(coords(d) - 1) < 1e-8)
2402 if (coords.array().minCoeff() >= 0 && coords.array().maxCoeff() <= 1)
2414 Eigen::MatrixXd samples_simplex, samples_cube, mapped, p0, p1, p;
2417 average_edge_length = 0;
2418 min_edge_length = std::numeric_limits<double>::max();
2420 if (!use_curved_mesh_size)
2424 min_edge_length = p.rowwise().norm().minCoeff();
2425 average_edge_length = p.rowwise().norm().mean();
2426 mesh_size = p.rowwise().norm().maxCoeff();
2428 logger().info(
"hmin: {}", min_edge_length);
2429 logger().info(
"hmax: {}", mesh_size);
2430 logger().info(
"havg: {}", average_edge_length);
2447 for (
size_t i = 0; i < bases_in.size(); ++i)
2456 bases_in[i].eval_geom_mapping(samples_simplex, mapped);
2461 bases_in[i].eval_geom_mapping(samples_cube, mapped);
2464 for (
int j = 0; j < n_edges; ++j)
2466 double current_edge = 0;
2467 for (
int k = 0; k < n_samples - 1; ++k)
2469 p0 = mapped.row(j * n_samples + k);
2470 p1 = mapped.row(j * n_samples + k + 1);
2473 current_edge += p.norm();
2476 mesh_size = std::max(current_edge, mesh_size);
2477 min_edge_length = std::min(current_edge, min_edge_length);
2478 average_edge_length += current_edge;
2483 average_edge_length /= n;
2485 logger().info(
"hmin: {}", min_edge_length);
2486 logger().info(
"hmax: {}", mesh_size);
2487 logger().info(
"havg: {}", average_edge_length);
2501 using namespace mesh;
2503 logger().info(
"Counting flipped elements...");
2507 for (
size_t i = 0; i < gbases.size(); ++i)
2513 if (!
vals.is_geom_mapping_positive(mesh.
is_volume(), gbases[i]))
2517 static const std::vector<std::string> element_type_names{{
2519 "RegularInteriorCube",
2520 "RegularBoundaryCube",
2521 "SimpleSingularInteriorCube",
2522 "MultiSingularInteriorCube",
2523 "SimpleSingularBoundaryCube",
2525 "MultiSingularBoundaryCube",
2531 log_and_throw_error(
"element {} is flipped, type {}", i, element_type_names[
static_cast<int>(els_tag[i])]);
2546 const std::vector<polyfem::basis::ElementBases> &bases,
2547 const std::vector<polyfem::basis::ElementBases> &gbases,
2551 const Eigen::MatrixXd &sol)
2555 logger().error(
"Build the bases first!");
2558 if (sol.size() <= 0)
2560 logger().error(
"Solve the problem first!");
2570 logger().info(
"Computing errors...");
2573 const int n_el = int(bases.size());
2575 Eigen::MatrixXd v_exact, v_approx;
2576 Eigen::MatrixXd v_exact_grad(0, 0), v_approx_grad;
2586 static const int p = 8;
2591 for (
int e = 0; e < n_el; ++e)
2601 v_approx.resize(
vals.val.rows(), actual_dim);
2604 v_approx_grad.resize(
vals.val.rows(), mesh.
dimension() * actual_dim);
2605 v_approx_grad.setZero();
2607 const int n_loc_bases = int(
vals.basis_values.size());
2609 for (
int i = 0; i < n_loc_bases; ++i)
2611 const auto &
val =
vals.basis_values[i];
2613 for (
size_t ii = 0; ii <
val.global.size(); ++ii)
2615 for (
int d = 0; d < actual_dim; ++d)
2617 v_approx.col(d) +=
val.global[ii].val * sol(
val.global[ii].index * actual_dim + d) *
val.val;
2618 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;
2623 const auto err = problem.
has_exact_sol() ? (v_exact - v_approx).eval().rowwise().norm().eval() : (v_approx).eval().rowwise().norm().eval();
2624 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();
2629 linf_err = std::max(linf_err, err.maxCoeff());
2630 grad_max_err = std::max(linf_err, err_grad.maxCoeff());
2672 l2_err += (err.array() * err.array() *
vals.det.array() *
vals.quadrature.weights.array()).sum();
2673 h1_err += (err_grad.array() * err_grad.array() *
vals.det.array() *
vals.quadrature.weights.array()).sum();
2674 lp_err += (err.array().pow(p) *
vals.det.array() *
vals.quadrature.weights.array()).sum();
2677 h1_semi_err = sqrt(fabs(h1_err));
2678 h1_err = sqrt(fabs(l2_err) + fabs(h1_err));
2679 l2_err = sqrt(fabs(l2_err));
2681 lp_err = pow(fabs(lp_err), 1. / p);
2686 const double computing_errors_time = timer.getElapsedTime();
2687 logger().info(
" took {}s", computing_errors_time);
2689 logger().info(
"-- L2 error: {}", l2_err);
2690 logger().info(
"-- Lp error: {}", lp_err);
2691 logger().info(
"-- H1 error: {}", h1_err);
2692 logger().info(
"-- H1 semi error: {}", h1_semi_err);
2695 logger().info(
"-- Linf error: {}", linf_err);
2696 logger().info(
"-- grad max error: {}", grad_max_err);
2711 regular_boundary_count = 0;
2712 simple_singular_count = 0;
2713 multi_singular_count = 0;
2715 non_regular_boundary_count = 0;
2716 non_regular_count = 0;
2717 undefined_count = 0;
2718 multi_singular_boundary_count = 0;
2722 for (
size_t i = 0; i < els_tag.size(); ++i)
2728 case ElementType::SIMPLEX:
2731 case ElementType::REGULAR_INTERIOR_CUBE:
2734 case ElementType::REGULAR_BOUNDARY_CUBE:
2735 regular_boundary_count++;
2737 case ElementType::SIMPLE_SINGULAR_INTERIOR_CUBE:
2738 simple_singular_count++;
2740 case ElementType::MULTI_SINGULAR_INTERIOR_CUBE:
2741 multi_singular_count++;
2743 case ElementType::SIMPLE_SINGULAR_BOUNDARY_CUBE:
2746 case ElementType::INTERFACE_CUBE:
2747 case ElementType::MULTI_SINGULAR_BOUNDARY_CUBE:
2748 multi_singular_boundary_count++;
2750 case ElementType::BOUNDARY_POLYTOPE:
2751 non_regular_boundary_count++;
2753 case ElementType::INTERIOR_POLYTOPE:
2754 non_regular_count++;
2756 case ElementType::UNDEFINED:
2760 throw std::runtime_error(
"Unknown element type");
2764 logger().info(
"simplex_count: \t{}", simplex_count);
2765 logger().info(
"regular_count: \t{}", regular_count);
2766 logger().info(
"regular_boundary_count: \t{}", regular_boundary_count);
2767 logger().info(
"simple_singular_count: \t{}", simple_singular_count);
2768 logger().info(
"multi_singular_count: \t{}", multi_singular_count);
2769 logger().info(
"boundary_count: \t{}", boundary_count);
2770 logger().info(
"multi_singular_boundary_count: \t{}", multi_singular_boundary_count);
2771 logger().info(
"non_regular_count: \t{}", non_regular_count);
2772 logger().info(
"non_regular_boundary_count: \t{}", non_regular_boundary_count);
2773 logger().info(
"undefined_count: \t{}", undefined_count);
2778 const nlohmann::json &args,
2779 const int n_bases,
const int n_pressure_bases,
2780 const Eigen::MatrixXd &sol,
2782 const Eigen::VectorXi &disc_orders,
2785 const std::string &formulation,
2786 const bool isoparametric,
2787 const int sol_at_node_id,
2793 j[
"geom_order"] = mesh.
orders().size() > 0 ? mesh.
orders().maxCoeff() : 1;
2794 j[
"geom_order_min"] = mesh.
orders().size() > 0 ? mesh.
orders().minCoeff() : 1;
2795 j[
"discr_order_min"] = disc_orders.minCoeff();
2796 j[
"discr_order_max"] = disc_orders.maxCoeff();
2797 j[
"iso_parametric"] = isoparametric;
2798 j[
"problem"] = problem.
name();
2799 j[
"mat_size"] = mat_size;
2800 j[
"num_bases"] = n_bases;
2801 j[
"num_pressure_bases"] = n_pressure_bases;
2802 j[
"num_non_zero"] = nn_zero;
2803 j[
"num_flipped"] = n_flipped;
2804 j[
"num_dofs"] = num_dofs;
2808 j[
"num_p1"] = (disc_orders.array() == 1).count();
2809 j[
"num_p2"] = (disc_orders.array() == 2).count();
2810 j[
"num_p3"] = (disc_orders.array() == 3).count();
2811 j[
"num_p4"] = (disc_orders.array() == 4).count();
2812 j[
"num_p5"] = (disc_orders.array() == 5).count();
2814 j[
"mesh_size"] = mesh_size;
2815 j[
"max_angle"] = max_angle;
2817 j[
"sigma_max"] = sigma_max;
2818 j[
"sigma_min"] = sigma_min;
2819 j[
"sigma_avg"] = sigma_avg;
2821 j[
"min_edge_length"] = min_edge_length;
2822 j[
"average_edge_length"] = average_edge_length;
2824 j[
"err_l2"] = l2_err;
2825 j[
"err_h1"] = h1_err;
2826 j[
"err_h1_semi"] = h1_semi_err;
2827 j[
"err_linf"] = linf_err;
2828 j[
"err_linf_grad"] = grad_max_err;
2829 j[
"err_lp"] = lp_err;
2831 j[
"spectrum"] = {spectrum(0), spectrum(1), spectrum(2), spectrum(3)};
2832 j[
"spectrum_condest"] = std::abs(spectrum(3)) / std::abs(spectrum(0));
2845 j[
"solver_info"] = solver_info;
2847 j[
"count_simplex"] = simplex_count;
2848 j[
"count_regular"] = regular_count;
2849 j[
"count_regular_boundary"] = regular_boundary_count;
2850 j[
"count_simple_singular"] = simple_singular_count;
2851 j[
"count_multi_singular"] = multi_singular_count;
2852 j[
"count_boundary"] = boundary_count;
2853 j[
"count_non_regular_boundary"] = non_regular_boundary_count;
2854 j[
"count_non_regular"] = non_regular_count;
2855 j[
"count_undefined"] = undefined_count;
2856 j[
"count_multi_singular_boundary"] = multi_singular_boundary_count;
2858 j[
"is_simplicial"] = mesh.
n_elements() == simplex_count;
2860 j[
"peak_memory"] =
getPeakRSS() / (1024 * 1024);
2864 std::vector<double> mmin(actual_dim);
2865 std::vector<double> mmax(actual_dim);
2867 for (
int d = 0; d < actual_dim; ++d)
2869 mmin[d] = std::numeric_limits<double>::max();
2870 mmax[d] = -std::numeric_limits<double>::max();
2873 for (
int i = 0; i < sol.size(); i += actual_dim)
2875 for (
int d = 0; d < actual_dim; ++d)
2877 mmin[d] = std::min(mmin[d], sol(i + d));
2878 mmax[d] = std::max(mmax[d], sol(i + d));
2882 std::vector<double> sol_at_node(actual_dim);
2884 if (sol_at_node_id >= 0)
2886 const int node_id = sol_at_node_id;
2888 for (
int d = 0; d < actual_dim; ++d)
2890 sol_at_node[d] = sol(node_id * actual_dim + d);
2894 j[
"sol_at_node"] = sol_at_node;
2895 j[
"sol_min"] = mmin;
2896 j[
"sol_max"] = mmax;
2898#if defined(POLYFEM_WITH_CPP_THREADS)
2900#elif defined(POLYFEM_WITH_TBB)
2903 j[
"num_threads"] = 1;
2906 j[
"formulation"] = formulation;
2912 : file(path), solve_data(solve_data)
2917 file << name <<
",";
2919 file <<
"total_energy" << std::endl;
2936 file << ((form && form->enabled()) ? form->value(sol) : 0) / s <<
",";
2943 : file(path), state(state), t0(t0), dt(dt)
2945 file <<
"step,time,forward,remeshing,global_relaxation,peak_mem,#V,#T" << std::endl;
2969 const double peak_mem =
getPeakRSS() / double(1 << 30);
2972 file << fmt::format(
2973 "{},{},{},{},{},{},{},{}\n",
2974 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 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_points(const std::string &path, const State &state, const Eigen::MatrixXd &sol, const ExportOptions &opts, std::vector< SolutionFrame > &solution_frames) const
saves the nodal values
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 save_wire(const std::string &name, const State &state, const Eigen::MatrixXd &sol, const double t, const ExportOptions &opts, std::vector< SolutionFrame > &solution_frames) const
saves the wireframe
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, std::vector< SolutionFrame > &solution_frames) const
saves the surface vtu file for for surface quantites, eg traction forces
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 build_grid(const polyfem::mesh::Mesh &mesh, const double spacing)
builds the grid to export the solution
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, std::vector< SolutionFrame > &solution_frames) const
saves the volume vtu file
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_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 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, std::vector< SolutionFrame > &solution_frames) const
exports everytihng, txt, vtu, etc
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, std::vector< SolutionFrame > &solution_frames) const
saves the vtu file for time t
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, std::vector< SolutionFrame > &solution_frames) const
saves the surface vtu file for for constact quantites, eg contact or friction forces
void init_sampler(const polyfem::mesh::Mesh &mesh, const double vismesh_rel_area)
unitalize the ref element sampler
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
bool is_linear() const
check if the mesh is linear
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::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
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)
std::string file_extension() const
return the extension of the output paraview files depending on use_hdf5
ExportOptions(const json &args, const bool is_mesh_linear, const bool is_problem_scalar, const bool solve_export_to_file)
initialize the flags based on the input args
bool discretization_order
bool solve_export_to_file