44#include <polysolve/linear/FEMSolver.hpp>
65 using namespace assembler;
68 using namespace utils;
73 void build_in_node_to_in_primitive(
const Mesh &mesh,
const MeshNodes &mesh_nodes,
74 Eigen::VectorXi &in_node_to_in_primitive,
75 Eigen::VectorXi &in_node_offset)
77 const int num_vertex_nodes = mesh_nodes.num_vertex_nodes();
78 const int num_edge_nodes = mesh_nodes.num_edge_nodes();
79 const int num_face_nodes = mesh_nodes.num_face_nodes();
80 const int num_cell_nodes = mesh_nodes.num_cell_nodes();
82 const int num_nodes = num_vertex_nodes + num_edge_nodes + num_face_nodes + num_cell_nodes;
84 const long n_vertices = num_vertex_nodes;
85 const int num_in_primitives = n_vertices + mesh.n_edges() + mesh.n_faces() + mesh.n_cells();
86 const int num_primitives = mesh.n_vertices() + mesh.n_edges() + mesh.n_faces() + mesh.n_cells();
88 in_node_to_in_primitive.resize(num_nodes);
89 in_node_offset.resize(num_nodes);
92 in_node_to_in_primitive.head(num_vertex_nodes).setLinSpaced(num_vertex_nodes, 0, num_vertex_nodes - 1);
93 in_node_offset.head(num_vertex_nodes).setZero();
95 int prim_offset = n_vertices;
96 int node_offset = num_vertex_nodes;
97 auto foo = [&](
const int num_prims,
const int num_prim_nodes) {
98 if (num_prims <= 0 || num_prim_nodes <= 0)
100 const Eigen::VectorXi range = Eigen::VectorXi::LinSpaced(num_prim_nodes, 0, num_prim_nodes - 1);
102 const int node_per_prim = num_prim_nodes / num_prims;
104 in_node_to_in_primitive.segment(node_offset, num_prim_nodes) =
105 range.array() / node_per_prim + prim_offset;
107 in_node_offset.segment(node_offset, num_prim_nodes) =
108 range.unaryExpr([&](
const int x) {
return x % node_per_prim; });
110 prim_offset += num_prims;
111 node_offset += num_prim_nodes;
114 foo(mesh.n_edges(), num_edge_nodes);
115 foo(mesh.n_faces(), num_face_nodes);
116 foo(mesh.n_cells(), num_cell_nodes);
119 bool build_in_primitive_to_primitive(
120 const Mesh &mesh,
const MeshNodes &mesh_nodes,
121 const Eigen::VectorXi &in_ordered_vertices,
122 const Eigen::MatrixXi &in_ordered_edges,
123 const Eigen::MatrixXi &in_ordered_faces,
124 Eigen::VectorXi &in_primitive_to_primitive)
127 const int num_vertex_nodes = mesh_nodes.num_vertex_nodes();
128 const int num_edge_nodes = mesh_nodes.num_edge_nodes();
129 const int num_face_nodes = mesh_nodes.num_face_nodes();
130 const int num_cell_nodes = mesh_nodes.num_cell_nodes();
131 const int num_nodes = num_vertex_nodes + num_edge_nodes + num_face_nodes + num_cell_nodes;
133 const long n_vertices = num_vertex_nodes;
134 const int num_in_primitives = n_vertices + mesh.n_edges() + mesh.n_faces() + mesh.n_cells();
135 const int num_primitives = mesh.n_vertices() + mesh.n_edges() + mesh.n_faces() + mesh.n_cells();
137 in_primitive_to_primitive.setLinSpaced(num_in_primitives, 0, num_in_primitives - 1);
145 if (in_ordered_vertices.rows() != n_vertices)
147 logger().warn(
"Node ordering disabled, in_ordered_vertices != n_vertices, {} != {}", in_ordered_vertices.rows(), n_vertices);
151 in_primitive_to_primitive.head(n_vertices) = in_ordered_vertices;
153 int in_offset = n_vertices;
154 int offset = mesh.n_vertices();
160 logger().trace(
"Building Mesh edges to IDs...");
162 const auto edges_to_ids = mesh.edges_to_ids();
163 if (in_ordered_edges.rows() != edges_to_ids.size())
165 logger().warn(
"Node ordering disabled, in_ordered_edges != edges_to_ids, {} != {}", in_ordered_edges.rows(), edges_to_ids.size());
169 logger().trace(
"Done (took {}s)", timer.getElapsedTime());
171 logger().trace(
"Building in-edge to edge mapping...");
173 for (
int in_ei = 0; in_ei < in_ordered_edges.rows(); in_ei++)
175 const std::pair<int, int> in_edge(
176 in_ordered_edges.row(in_ei).minCoeff(),
177 in_ordered_edges.row(in_ei).maxCoeff());
178 in_primitive_to_primitive[in_offset + in_ei] =
179 offset + edges_to_ids.at(in_edge);
182 logger().trace(
"Done (took {}s)", timer.getElapsedTime());
184 in_offset += mesh.n_edges();
185 offset += mesh.n_edges();
191 if (mesh.is_volume())
193 logger().trace(
"Building Mesh faces to IDs...");
195 const auto faces_to_ids = mesh.faces_to_ids();
196 if (in_ordered_faces.rows() != faces_to_ids.size())
198 logger().warn(
"Node ordering disabled, in_ordered_faces != faces_to_ids, {} != {}", in_ordered_faces.rows(), faces_to_ids.size());
202 logger().trace(
"Done (took {}s)", timer.getElapsedTime());
204 logger().trace(
"Building in-face to face mapping...");
206 for (
int in_fi = 0; in_fi < in_ordered_faces.rows(); in_fi++)
208 std::vector<int> in_face(in_ordered_faces.cols());
209 for (
int i = 0; i < in_face.size(); i++)
210 in_face[i] = in_ordered_faces(in_fi, i);
211 std::sort(in_face.begin(), in_face.end());
213 in_primitive_to_primitive[in_offset + in_fi] =
214 offset + faces_to_ids.at(in_face);
217 logger().trace(
"Done (took {}s)", timer.getElapsedTime());
219 in_offset += mesh.n_faces();
220 offset += mesh.n_faces();
230 indices.resize(
mesh->n_vertices());
237 std::vector<int> indices;
239 for (
int i = 0; i < p2n.size(); i++)
246 if (
args[
"space"][
"basis_type"] ==
"Spline")
248 logger().warn(
"Node ordering disabled, it dosent work for splines!");
254 logger().warn(
"Node ordering disabled, it works only for p < 4 and uniform order!");
258 if (!
mesh->is_conforming())
260 logger().warn(
"Node ordering disabled, not supported for non-conforming meshes!");
264 if (
mesh->has_poly())
266 logger().warn(
"Node ordering disabled, not supported for polygonal meshes!");
270 if (
mesh->in_ordered_vertices().size() <= 0 ||
mesh->in_ordered_edges().size() <= 0 || (
mesh->is_volume() &&
mesh->in_ordered_faces().size() <= 0))
272 logger().warn(
"Node ordering disabled, input vertices/edges/faces not computed!");
276 const int num_vertex_nodes =
mesh_nodes->num_vertex_nodes();
277 const int num_edge_nodes =
mesh_nodes->num_edge_nodes();
278 const int num_face_nodes =
mesh_nodes->num_face_nodes();
279 const int num_cell_nodes =
mesh_nodes->num_cell_nodes();
281 const int num_nodes = num_vertex_nodes + num_edge_nodes + num_face_nodes + num_cell_nodes;
283 const long n_vertices = num_vertex_nodes;
284 const int num_in_primitives = n_vertices +
mesh->n_edges() +
mesh->n_faces() +
mesh->n_cells();
285 const int num_primitives =
mesh->n_vertices() +
mesh->n_edges() +
mesh->n_faces() +
mesh->n_cells();
289 logger().trace(
"Building in-node to in-primitive mapping...");
291 Eigen::VectorXi in_node_to_in_primitive;
292 Eigen::VectorXi in_node_offset;
293 build_in_node_to_in_primitive(*
mesh, *
mesh_nodes, in_node_to_in_primitive, in_node_offset);
295 logger().trace(
"Done (took {}s)", timer.getElapsedTime());
297 logger().trace(
"Building in-primitive to primitive mapping...");
299 bool ok = build_in_primitive_to_primitive(
301 mesh->in_ordered_vertices(),
302 mesh->in_ordered_edges(),
303 mesh->in_ordered_faces(),
306 logger().trace(
"Done (took {}s)", timer.getElapsedTime());
315 const auto primitive_offset = [&](
int node) {
319 return mesh->n_vertices();
321 return mesh->n_vertices() +
mesh->n_edges();
323 return mesh->n_vertices() +
mesh->n_edges() +
mesh->n_faces();
324 throw std::runtime_error(
"Invalid node ID!");
327 logger().trace(
"Building primitive to node mapping...");
329 std::vector<std::vector<int>> primitive_to_nodes(num_primitives);
330 const std::vector<int> &grouped_nodes =
mesh_nodes->primitive_to_node();
332 for (
int i = 0; i < grouped_nodes.size(); i++)
334 int node = grouped_nodes[i];
335 assert(node < num_nodes);
338 int primitive =
mesh_nodes->node_to_primitive_gid().at(node) + primitive_offset(i);
339 assert(primitive < num_primitives);
340 primitive_to_nodes[primitive].push_back(node);
344 assert(node_count == num_nodes);
346 logger().trace(
"Done (took {}s)", timer.getElapsedTime());
348 logger().trace(
"Combining mappings...");
351 for (
int i = 0; i < num_nodes; i++)
354 const std::vector<int> &possible_nodes =
357 if (possible_nodes.size() == 1)
361 assert(possible_nodes.size() > 1);
367 assert(e_id < mesh->n_edges());
376 assert(
mesh->edge_vertex(e_id, 0) <
mesh->edge_vertex(e_id, 1));
377 int offset = (a - v0).squaredNorm() < (b - v0).squaredNorm()
379 : (possible_nodes.size() - in_node_offset[i] - 1);
384 logger().trace(
"Done (took {}s)", timer.getElapsedTime());
389 if (
args[
"materials"].is_null())
391 logger().error(
"specify some 'materials'");
392 assert(!
args[
"materials"].is_null());
393 throw "invalid input";
396 if (
args[
"materials"].is_array())
398 std::string current =
"";
399 for (
const auto &m :
args[
"materials"])
401 const std::string tmp = m[
"type"];
404 else if (current != tmp)
406 if (current ==
"LinearElasticity"
407 || current ==
"NeoHookean"
408 || current ==
"SaintVenant"
409 || current ==
"HookeLinearElasticity"
410 || current ==
"MooneyRivlin"
411 || current ==
"MooneyRivlin3Param"
412 || current ==
"UnconstrainedOgden"
413 || current ==
"IncompressibleOgden"
414 || current ==
"MultiModels")
416 if (tmp ==
"LinearElasticity"
417 || tmp ==
"NeoHookean"
418 || tmp ==
"SaintVenant"
419 || tmp ==
"HookeLinearElasticity"
420 || tmp ==
"MooneyRivlin"
421 || current ==
"MooneyRivlin3Param"
422 || tmp ==
"UnconstrainedOgden"
423 || tmp ==
"IncompressibleOgden")
424 current =
"MultiModels";
427 logger().error(
"Current material is {}, new material is {}, multimaterial supported only for LinearElasticity and NeoHookean", current, tmp);
428 throw "invalid input";
433 logger().error(
"Current material is {}, new material is {}, multimaterial supported only for LinearElasticity and NeoHookean", current, tmp);
434 throw "invalid input";
442 return args[
"materials"][
"type"];
449 logger().error(
"No pressure bases defined!");
454 Eigen::MatrixXd tmp = sol;
466 const std::vector<basis::ElementBases> &bases,
467 const std::vector<basis::ElementBases> &gbases,
468 Eigen::MatrixXd &basis_integrals)
472 logger().error(
"Works only on volumetric meshes!");
477 basis_integrals.resize(n_bases, 9);
478 basis_integrals.setZero();
479 Eigen::MatrixXd rhs(n_bases, 9);
483 for (
int e = 0; e < n_elements; ++e)
495 for (
int j = 0; j < n_local_bases; ++j)
512 for (
size_t ii = 0; ii < v.
global.size(); ++ii)
514 basis_integrals(v.
global[ii].index, 0) += integral_100 * v.
global[ii].val;
515 basis_integrals(v.
global[ii].index, 1) += integral_010 * v.
global[ii].val;
516 basis_integrals(v.
global[ii].index, 2) += integral_001 * v.
global[ii].val;
518 basis_integrals(v.
global[ii].index, 3) += integral_110 * v.
global[ii].val;
519 basis_integrals(v.
global[ii].index, 4) += integral_011 * v.
global[ii].val;
520 basis_integrals(v.
global[ii].index, 5) += integral_101 * v.
global[ii].val;
522 basis_integrals(v.
global[ii].index, 6) += integral_200 * v.
global[ii].val;
523 basis_integrals(v.
global[ii].index, 7) += integral_020 * v.
global[ii].val;
524 basis_integrals(v.
global[ii].index, 8) += integral_002 * v.
global[ii].val;
526 rhs(v.
global[ii].index, 6) += -2.0 * area * v.
global[ii].val;
527 rhs(v.
global[ii].index, 7) += -2.0 * area * v.
global[ii].val;
528 rhs(v.
global[ii].index, 8) += -2.0 * area * v.
global[ii].val;
533 basis_integrals -= rhs;
538 if (
mesh->has_poly())
541 if (
args[
"space"][
"basis_type"] ==
"Spline")
544 if (
mesh->is_rational())
547 if (
args[
"space"][
"use_p_ref"])
556 if (
mesh->orders().size() <= 0)
558 if (
args[
"space"][
"discr_order"] == 1)
561 return args[
"space"][
"advanced"][
"isoparametric"];
564 if (
mesh->orders().minCoeff() !=
mesh->orders().maxCoeff())
567 if (
args[
"space"][
"discr_order"] ==
mesh->orders().minCoeff())
574 return args[
"space"][
"advanced"][
"isoparametric"];
581 logger().error(
"Load the mesh first!");
585 mesh->prepare_mesh();
604 assert(
args[
"materials"].is_array());
606 std::vector<std::string> materials(
mesh->n_elements());
608 std::map<int, std::string> mats;
610 for (
const auto &m :
args[
"materials"])
611 mats[m[
"id"].get<int>()] = m[
"type"];
613 for (
int i = 0; i < materials.size(); ++i)
614 materials[i] = mats.at(
mesh->get_body_id(i));
616 mm->init_multimodels(materials);
628 logger().info(
"Building {} basis...", (
iso_parametric() ?
"isoparametric" :
"not isoparametric"));
629 const bool has_polys =
mesh->has_poly();
635 std::map<int, basis::InterfaceData> poly_edge_to_data_geom;
637 const auto &tmp_json =
args[
"space"][
"discr_order"];
638 if (tmp_json.is_number_integer())
642 else if (tmp_json.is_string())
644 const std::string discr_orders_path = tmp_json;
648 assert(tmp.cols() == 1);
651 else if (tmp_json.is_array())
653 const auto b_discr_orders = tmp_json;
655 std::map<int, int> b_orders;
656 for (
size_t i = 0; i < b_discr_orders.size(); ++i)
658 assert(b_discr_orders[i][
"id"].is_array() || b_discr_orders[i][
"id"].is_number_integer());
660 const int order = b_discr_orders[i][
"order"];
661 for (
const int id : json_as_array<int>(b_discr_orders[i][
"id"]))
663 b_orders[id] = order;
664 logger().trace(
"bid {}, discr {}",
id, order);
668 for (
int e = 0; e <
mesh->n_elements(); ++e)
670 const int bid =
mesh->get_body_id(e);
671 const auto order = b_orders.find(bid);
672 if (order == b_orders.end())
674 logger().debug(
"Missing discretization order for body {}; using 1", bid);
686 logger().error(
"space/discr_order must be either a number a path or an array");
687 throw "invalid json";
691 Eigen::MatrixXi geom_disc_orders;
694 if (
mesh->orders().size() <= 0)
697 geom_disc_orders.setConstant(1);
700 geom_disc_orders =
mesh->orders();
705 if (
args[
"space"][
"use_p_ref"])
709 args[
"space"][
"advanced"][
"B"],
710 args[
"space"][
"advanced"][
"h1_formula"],
711 args[
"space"][
"discr_order"],
712 args[
"space"][
"advanced"][
"discr_order_max"],
719 int quadrature_order =
args[
"space"][
"advanced"][
"quadrature_order"].get<
int>();
720 const int mass_quadrature_order =
args[
"space"][
"advanced"][
"mass_quadrature_order"].get<
int>();
726 logger().error(
"p refinement not supported in mixed formulation!");
733 const bool use_continuous_gbasis =
true;
735 if (
mesh->is_volume())
738 if (
args[
"space"][
"basis_type"] ==
"Spline")
755 n_geom_bases =
basis::LagrangeBasis3d::build_bases(tmp_mesh,
assembler->name(), quadrature_order, mass_quadrature_order, geom_disc_orders,
false, has_polys, !use_continuous_gbasis,
geom_bases_,
local_boundary, poly_edge_to_data_geom,
geom_mesh_nodes);
757 n_bases =
basis::LagrangeBasis3d::build_bases(tmp_mesh,
assembler->name(), quadrature_order, mass_quadrature_order,
disc_orders,
args[
"space"][
"basis_type"] ==
"Serendipity", has_polys,
false,
bases,
local_boundary,
poly_edge_to_data,
mesh_nodes);
763 n_pressure_bases =
basis::LagrangeBasis3d::build_bases(tmp_mesh,
assembler->name(), quadrature_order, mass_quadrature_order,
int(
args[
"space"][
"pressure_discr_order"]),
false, has_polys,
false,
pressure_bases,
local_boundary, poly_edge_to_data_geom,
pressure_mesh_nodes);
769 if (
args[
"space"][
"basis_type"] ==
"Spline")
787 n_geom_bases =
basis::LagrangeBasis2d::build_bases(tmp_mesh,
assembler->name(), quadrature_order, mass_quadrature_order, geom_disc_orders,
false, has_polys, !use_continuous_gbasis,
geom_bases_,
local_boundary, poly_edge_to_data_geom,
geom_mesh_nodes);
789 n_bases =
basis::LagrangeBasis2d::build_bases(tmp_mesh,
assembler->name(), quadrature_order, mass_quadrature_order,
disc_orders,
args[
"space"][
"basis_type"] ==
"Serendipity", has_polys,
false,
bases,
local_boundary,
poly_edge_to_data,
mesh_nodes);
795 n_pressure_bases =
basis::LagrangeBasis2d::build_bases(tmp_mesh,
assembler->name(), quadrature_order, mass_quadrature_order,
int(
args[
"space"][
"pressure_discr_order"]),
false, has_polys,
false,
pressure_bases,
local_boundary, poly_edge_to_data_geom,
pressure_mesh_nodes);
805 bases[i].compute_quadrature(b_quad);
821 std::map<std::array<int, 2>,
double> pairs;
822 for (
int e = 0; e < gbases.size(); e++)
824 const auto &gbs = gbases[e].bases;
825 const auto &bs =
bases[e].bases;
827 Eigen::MatrixXd local_pts;
828 const int order = bs.front().order();
829 if (
mesh->is_volume())
831 if (
mesh->is_simplex(e))
838 if (
mesh->is_simplex(e))
847 for (
int i = 0; i < bs.size(); i++)
849 for (
int j = 0; j < gbs.size(); j++)
851 if (std::abs(
vals.basis_values[j].val(i)) > 1e-7)
853 std::array<int, 2> index = {{gbs[j].global()[0].index, bs[i].global()[0].index}};
854 pairs.insert({index,
vals.basis_values[j].val(i)});
860 const int dim =
mesh->dimension();
861 std::vector<Eigen::Triplet<double>> coeffs;
863 for (
const auto &iter : pairs)
864 for (
int d = 0; d < dim; d++)
865 coeffs.emplace_back(iter.first[0] * dim + d, iter.first[1] * dim + d, iter.second);
874 const int dim =
mesh->dimension();
875 const int problem_dim =
problem->is_scalar() ? 1 : dim;
881 json directions =
args[
"boundary_conditions"][
"periodic_boundary"][
"correspondence"];
882 Eigen::MatrixXd tile_offset = Eigen::MatrixXd::Identity(dim, dim);
884 if (directions.size() > 0)
887 for (
int d = 0; d < dim; d++)
890 if (tmp.size() != dim)
892 tile_offset.col(d) = tmp;
901 if (
args[
"space"][
"advanced"][
"count_flipped_els"])
904 const int prev_bases =
n_bases;
909 logger().debug(
"Building node mapping...");
915 logger().debug(
"Done (took {}s)", timer2.getElapsedTime());
918 logger().info(
"Building collision mesh...");
942 for (
const auto &bs :
bases)
944 for (
const auto &b : bs.bases)
946 for (
const auto &lg : b.global())
948 if (lg.index == n_id)
972 for (
const auto &bs :
bases)
974 for (
const auto &b : bs.bases)
976 for (
const auto &lg : b.global())
978 if (lg.index == n_id)
1001 for (
int i = prev_bases; i <
n_bases; ++i)
1003 for (
int d = 0; d < problem_dim; ++d)
1015 std::vector<bool> isboundary(
n_bases,
false);
1018 const int e = lb.element_id();
1019 for (
int i = 0; i < lb.size(); ++i)
1021 const auto nodes =
bases[e].local_nodes_for_primitive(lb.global_primitive_id(i), *
mesh);
1024 isboundary[
bases[e].
bases[n].global()[0].index] =
true;
1033 const int actual_dim =
problem->is_scalar() ? 1 :
mesh->dimension();
1034 for (
int d = 0; d < actual_dim; d++)
1038 logger().info(
"Fix solution at node {} to remove singularity due to periodic BC", i);
1042 const int n_samples = 10;
1065 args[
"contact"][
"dhat"] = dhat;
1070 logger().info(
"dhat set to {}",
double(
args[
"contact"][
"dhat"]));
1091 if (
n_bases <=
args[
"solver"][
"advanced"][
"cache_size"])
1094 logger().info(
"Building cache...");
1100 logger().info(
" took {}s", timer.getElapsedTime());
1108 logger().warn(
"(Quasi-)Static problem without Dirichlet nodes, will fix solution at one node to find a unique solution!");
1118 logger().error(
"Load the mesh first!");
1132 logger().info(
"Computing polygonal basis...");
1143 if (
mesh->is_volume())
1145 if (
args[
"space"][
"poly_basis_type"] ==
"MeanValue" ||
args[
"space"][
"poly_basis_type"] ==
"Wachspress")
1146 logger().error(
"Barycentric bases not supported in 3D");
1150 args[
"space"][
"advanced"][
"n_harmonic_samples"],
1153 args[
"space"][
"advanced"][
"quadrature_order"],
1154 args[
"space"][
"advanced"][
"mass_quadrature_order"],
1155 args[
"space"][
"advanced"][
"integral_constraints"],
1163 if (
args[
"space"][
"poly_basis_type"] ==
"MeanValue")
1170 args[
"space"][
"advanced"][
"quadrature_order"],
1171 args[
"space"][
"advanced"][
"mass_quadrature_order"],
1174 else if (
args[
"space"][
"poly_basis_type"] ==
"Wachspress")
1181 args[
"space"][
"advanced"][
"quadrature_order"],
1182 args[
"space"][
"advanced"][
"mass_quadrature_order"],
1190 args[
"space"][
"advanced"][
"n_harmonic_samples"],
1193 args[
"space"][
"advanced"][
"quadrature_order"],
1194 args[
"space"][
"advanced"][
"mass_quadrature_order"],
1195 args[
"space"][
"advanced"][
"integral_constraints"],
1205 if (
mesh->is_volume())
1207 if (
args[
"space"][
"poly_basis_type"] ==
"MeanValue" ||
args[
"space"][
"poly_basis_type"] ==
"Wachspress")
1209 logger().error(
"Barycentric bases not supported in 3D");
1210 throw "not implemented";
1217 args[
"space"][
"advanced"][
"n_harmonic_samples"],
1220 args[
"space"][
"advanced"][
"quadrature_order"],
1221 args[
"space"][
"advanced"][
"mass_quadrature_order"],
1222 args[
"space"][
"advanced"][
"integral_constraints"],
1231 if (
args[
"space"][
"poly_basis_type"] ==
"MeanValue")
1237 n_bases,
args[
"space"][
"advanced"][
"quadrature_order"],
1238 args[
"space"][
"advanced"][
"mass_quadrature_order"],
1241 else if (
args[
"space"][
"poly_basis_type"] ==
"Wachspress")
1247 n_bases,
args[
"space"][
"advanced"][
"quadrature_order"],
1248 args[
"space"][
"advanced"][
"mass_quadrature_order"],
1256 args[
"space"][
"advanced"][
"n_harmonic_samples"],
1259 args[
"space"][
"advanced"][
"quadrature_order"],
1260 args[
"space"][
"advanced"][
"mass_quadrature_order"],
1261 args[
"space"][
"advanced"][
"integral_constraints"],
1279 assert(!
mesh->is_volume());
1280 const int dim =
mesh->dimension();
1281 const int n_tiles = 2;
1283 if (
mesh->dimension() != 2)
1287 for (
const auto &bs :
bases)
1288 for (
const auto &b : bs.bases)
1289 for (
const auto &g : b.global())
1290 V.row(g.index) = g.node;
1293 for (
int i = 0; i < E.size(); i++)
1296 Eigen::MatrixXd bbox(
V.cols(), 2);
1297 bbox.col(0) =
V.colwise().minCoeff();
1298 bbox.col(1) =
V.colwise().maxCoeff();
1302 std::vector<int> ind;
1303 for (
int i = 0; i < E.rows(); i++)
1308 E = E(ind, Eigen::all).eval();
1311 Eigen::MatrixXd Vtmp, Vnew;
1312 Eigen::MatrixXi Etmp, Enew;
1313 Vtmp.setZero(
V.rows() * n_tiles * n_tiles,
V.cols());
1314 Etmp.setZero(E.rows() * n_tiles * n_tiles, E.cols());
1316 Eigen::MatrixXd tile_offset =
periodic_bc->get_affine_matrix();
1318 for (
int i = 0, idx = 0; i < n_tiles; i++)
1320 for (
int j = 0; j < n_tiles; j++)
1322 Eigen::Vector2d block_id;
1325 Vtmp.middleRows(idx *
V.rows(),
V.rows()) =
V;
1328 for (
int vid = 0; vid <
V.rows(); vid++)
1329 Vtmp.block(idx *
V.rows() + vid, 0, 1, 2) += (tile_offset * block_id).transpose();
1331 Etmp.middleRows(idx * E.rows(), E.rows()) = E.array() + idx *
V.rows();
1337 Eigen::VectorXi indices;
1339 std::vector<int> tmp;
1340 for (
int i = 0; i <
V.rows(); i++)
1346 indices.resize(tmp.size() * n_tiles * n_tiles);
1347 for (
int i = 0; i < n_tiles * n_tiles; i++)
1349 indices.segment(i * tmp.size(), tmp.size()) = Eigen::Map<Eigen::VectorXi, Eigen::Unaligned>(tmp.data(), tmp.size());
1350 indices.segment(i * tmp.size(), tmp.size()).array() += i *
V.rows();
1354 Eigen::VectorXi potentially_duplicate_mask(Vtmp.rows());
1355 potentially_duplicate_mask.setZero();
1356 potentially_duplicate_mask(indices).array() = 1;
1357 Eigen::MatrixXd candidates = Vtmp(indices, Eigen::all);
1359 Eigen::VectorXi SVI;
1360 std::vector<int> SVJ;
1361 SVI.setConstant(Vtmp.rows(), -1);
1363 const double eps = (bbox.col(1) - bbox.col(0)).maxCoeff() *
args[
"boundary_conditions"][
"periodic_boundary"][
"tolerance"].get<
double>();
1364 for (
int i = 0; i < Vtmp.rows(); i++)
1370 if (potentially_duplicate_mask(i))
1372 Eigen::VectorXd diffs = (candidates.rowwise() - Vtmp.row(i)).rowwise().norm();
1373 for (
int j = 0; j < diffs.size(); j++)
1375 SVI[indices[j]] =
id;
1380 Vnew = Vtmp(SVJ, Eigen::all);
1382 Enew.resizeLike(Etmp);
1383 for (
int d = 0; d < Etmp.cols(); d++)
1384 Enew.col(d) = SVI(Etmp.col(d));
1386 std::vector<bool> is_on_surface = ipc::CollisionMesh::construct_is_on_surface(Vnew.rows(), Enew);
1388 Eigen::MatrixXi boundary_triangles;
1389 Eigen::SparseMatrix<double> displacement_map;
1399 for (
int i = 0; i <
V.rows(); i++)
1400 for (
int j = 0; j < n_tiles * n_tiles; j++)
1418 const std::vector<basis::ElementBases> &bases,
1419 const std::vector<basis::ElementBases> &geom_bases,
1420 const std::vector<mesh::LocalBoundary> &total_local_boundary,
1423 const std::function<std::string(
const std::string &)> &resolve_input_path,
1424 const Eigen::VectorXi &in_node_to_node,
1425 ipc::CollisionMesh &collision_mesh)
1427 Eigen::MatrixXd collision_vertices;
1428 Eigen::VectorXi collision_codim_vids;
1429 Eigen::MatrixXi collision_edges, collision_triangles;
1430 std::vector<Eigen::Triplet<double>> displacement_map_entries;
1432 if (
args.contains(
"/contact/collision_mesh"_json_pointer)
1433 &&
args.at(
"/contact/collision_mesh/enabled"_json_pointer).get<
bool>())
1435 const json collision_mesh_args =
args.at(
"/contact/collision_mesh"_json_pointer);
1436 if (collision_mesh_args.contains(
"linear_map"))
1438 assert(displacement_map_entries.empty());
1439 assert(collision_mesh_args.contains(
"mesh"));
1440 const std::string
root_path = utils::json_value<std::string>(
args,
"root_path",
"");
1446 in_node_to_node, transformation, collision_vertices, collision_codim_vids,
1447 collision_edges, collision_triangles, displacement_map_entries);
1451 assert(collision_mesh_args.contains(
"max_edge_length"));
1453 "Building collision proxy with max edge length={} ...",
1454 collision_mesh_args[
"max_edge_length"].get<double>());
1459 collision_mesh_args[
"max_edge_length"], collision_vertices,
1460 collision_triangles, displacement_map_entries,
1461 collision_mesh_args[
"tessellation_type"]);
1462 if (collision_triangles.size())
1463 igl::edges(collision_triangles, collision_edges);
1465 logger().debug(fmt::format(
1466 std::locale(
"en_US.UTF-8"),
1467 "Done (took {:g}s, {:L} vertices, {:L} triangles)",
1468 timer.getElapsedTime(),
1469 collision_vertices.rows(), collision_triangles.rows()));
1476 collision_vertices, collision_edges, collision_triangles, displacement_map_entries);
1481 const int num_fe_collision_vertices = collision_vertices.rows();
1482 assert(collision_edges.size() == 0 || collision_edges.maxCoeff() < num_fe_collision_vertices);
1483 assert(collision_triangles.size() == 0 || collision_triangles.maxCoeff() < num_fe_collision_vertices);
1493 if (!displacement_map_entries.empty())
1495 displacement_map_entries.reserve(displacement_map_entries.size() +
obstacle.
n_vertices());
1498 displacement_map_entries.emplace_back(num_fe_collision_vertices + i, num_fe_nodes + i, 1.0);
1503 std::vector<bool> is_on_surface = ipc::CollisionMesh::construct_is_on_surface(
1504 collision_vertices.rows(), collision_edges);
1505 for (
const int vid : collision_codim_vids)
1507 is_on_surface[vid] =
true;
1510 Eigen::SparseMatrix<double> displacement_map;
1511 if (!displacement_map_entries.empty())
1513 displacement_map.resize(collision_vertices.rows(),
n_bases);
1514 displacement_map.setFromTriplets(displacement_map_entries.begin(), displacement_map_entries.end());
1518 is_on_surface, collision_vertices, collision_edges, collision_triangles,
1523 return collision_mesh.to_full_vertex_id(vi) < num_fe_collision_vertices
1524 ||
collision_mesh.to_full_vertex_id(vj) < num_fe_collision_vertices;
1534 logger().error(
"Load the mesh first!");
1539 logger().error(
"Build the bases first!");
1542 if (
assembler->name() ==
"OperatorSplitting")
1549 if (!
problem->is_time_dependent())
1560 logger().info(
"Assembling mass mat...");
1567 std::vector<Eigen::Triplet<double>> mass_blocks;
1568 mass_blocks.reserve(velocity_mass.nonZeros());
1570 for (
int k = 0; k < velocity_mass.outerSize(); ++k)
1572 for (StiffnessMatrix::InnerIterator it(velocity_mass, k); it; ++it)
1574 mass_blocks.emplace_back(it.row(), it.col(), it.value());
1579 mass.setFromTriplets(mass_blocks.begin(), mass_blocks.end());
1580 mass.makeCompressed();
1587 assert(
mass.size() > 0);
1590 for (
int k = 0; k <
mass.outerSize(); ++k)
1593 for (StiffnessMatrix::InnerIterator it(
mass, k); it; ++it)
1595 assert(it.col() == k);
1603 if (
args[
"solver"][
"advanced"][
"lump_mass_matrix"])
1620 const std::vector<basis::ElementBases> &bases_,
1623 json rhs_solver_params =
args[
"solver"][
"linear"];
1624 if (!rhs_solver_params.contains(
"Pardiso"))
1625 rhs_solver_params[
"Pardiso"] = {};
1626 rhs_solver_params[
"Pardiso"][
"mtype"] = -2;
1628 const int size =
problem->is_scalar() ? 1 :
mesh->dimension();
1630 return std::make_shared<RhsAssembler>(
1635 args[
"space"][
"advanced"][
"bc_method"],
1641 const std::vector<basis::ElementBases> &bases_)
const
1643 const int size =
problem->is_scalar() ? 1 :
mesh->dimension();
1645 return std::make_shared<PressureAssembler>(
1658 logger().error(
"Load the mesh first!");
1663 logger().error(
"Build the bases first!");
1670 p_params[
"formulation"] =
assembler->name();
1674 mesh->bounding_box(min, max);
1675 delta = (max - min) / 2. + min;
1676 if (
mesh->is_volume())
1677 p_params[
"bbox_center"] = {delta(0), delta(1), delta(2)};
1679 p_params[
"bbox_center"] = {delta(0), delta(1)};
1681 problem->set_parameters(p_params);
1686 logger().info(
"Assigning rhs...");
1695 const int prev_size =
rhs.size();
1697 rhs.conservativeResize(prev_size + n_larger,
rhs.cols());
1698 if (
assembler->name() ==
"OperatorSplitting")
1706 rhs.block(prev_size, 0, n_larger,
rhs.cols()).setZero();
1717 rhs.block(prev_size, 0, n_larger,
rhs.cols()) = tmp;
1730 logger().error(
"Load the mesh first!");
1735 logger().error(
"Build the bases first!");
1755 if (
problem->is_time_dependent())
1757 const double t0 =
args[
"time"][
"t0"];
1758 const int time_steps =
args[
"time"][
"time_steps"];
1759 const double dt =
args[
"time"][
"dt"];
1762 if (
args[
"output"][
"advanced"][
"save_time_sequence"])
1764 logger().info(
"Time sequence of simulation will be written to: \"{}\"",
1768 if (
assembler->name() ==
"NavierStokes")
1770 else if (
assembler->name() ==
"OperatorSplitting")
1777 throw std::runtime_error(
"Nonlinear scalar problems are not supported yet!");
1783 if (
assembler->name() ==
"NavierStokes")
1795 throw std::runtime_error(
"Nonlinear scalar problems are not supported yet!");
1804 if (!state_path.empty())
ElementAssemblyValues vals
DECLARE_DIFFSCALAR_BASE()
io::OutRuntimeData timings
runtime statistics
assembler::MacroStrainValue macro_strain_constraint
std::shared_ptr< utils::PeriodicBoundary > periodic_bc
periodic BC and periodic mesh utils
int n_bases
number of bases
void cache_transient_adjoint_quantities(const int current_step, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &disp_grad)
int n_pressure_bases
number of pressure bases
std::string formulation() const
return the formulation (checks if the problem is scalar or not and deals with multiphysics)
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
int n_boundary_samples() const
quadrature used for projecting boundary conditions
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
double starting_max_edge_length
std::vector< mesh::LocalBoundary > local_pressure_boundary
mapping from elements to nodes for pressure boundary conditions
void sol_to_pressure(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
splits the solution in solution and pressure for mixed problems
std::string root_path() const
Get the root path for the state (e.g., args["root_path"] or ".")
void assemble_mass_mat()
assemble mass, step 4 of solve build mass matrix based on defined basis modifies mass (and maybe more...
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
void build_polygonal_basis()
builds bases for polygons, called inside build_basis
mesh::Obstacle obstacle
Obstacles used in collisions.
std::shared_ptr< assembler::Assembler > assembler
assemblers
ipc::CollisionMesh periodic_collision_mesh
IPC collision mesh under periodic BC.
void build_periodic_collision_mesh()
io::OutGeometryData out_geom
visualization stuff
bool use_avg_pressure
use average pressure for stokes problem to fix the additional dofs, true by default if false,...
std::string resolve_output_path(const std::string &path) const
Resolve output path relative to output_dir if the path is not absolute.
std::vector< int > neumann_nodes
per node neumann
ipc::CollisionMesh collision_mesh
IPC collision mesh.
std::vector< basis::ElementBases > geom_bases_
Geometric mapping bases, if the elements are isoparametric, this list is empty.
solver::CacheLevel optimization_enabled
Eigen::VectorXi in_primitive_to_primitive
maps in vertices/edges/faces/cells to polyfem vertices/edges/faces/cells
std::map< int, basis::InterfaceData > poly_edge_to_data
nodes on the boundary of polygonal elements, used for harmonic bases
std::vector< int > pressure_boundary_nodes
list of neumann boundary nodes
std::vector< basis::ElementBases > pressure_bases
FE pressure bases for mixed elements, the size is #elements.
bool is_homogenization() const
void solve_transient_tensor_nonlinear(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol)
solves transient tensor nonlinear problem
double starting_min_edge_length
std::shared_ptr< assembler::Mass > mass_matrix_assembler
void build_node_mapping()
build the mapping from input nodes to polyfem nodes
std::vector< int > primitive_to_node() const
bool has_dhat
stores if input json contains dhat
std::shared_ptr< polyfem::mesh::MeshNodes > pressure_mesh_nodes
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
std::shared_ptr< polyfem::mesh::MeshNodes > mesh_nodes
Mapping from input nodes to FE nodes.
StiffnessMatrix mass
Mass matrix, it is computed only for time dependent problems.
void solve_homogenization(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol)
std::vector< int > dirichlet_nodes
per node dirichlet
StiffnessMatrix gbasis_nodes_to_basis_nodes
bool has_periodic_bc() const
void solve_problem(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves the problems
void build_basis()
builds the bases step 2 of solve modifies bases, pressure_bases, geom_bases_, boundary_nodes,...
void solve_navier_stokes(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves a navier stokes
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
std::vector< int > node_to_primitive() const
json args
main input arguments containing all defaults
io::OutStatsData stats
Other statistics.
std::vector< RowVectorNd > neumann_nodes_position
void assemble_rhs()
compute rhs, step 3 of solve build rhs vector based on defined basis and given rhs of the problem mod...
double min_boundary_edge_length
Eigen::VectorXi periodic_collision_mesh_to_basis
index mapping from periodic 2x2 collision mesh to FE periodic mesh
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
void solve_transient_navier_stokes_split(const int time_steps, const double dt, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves transient navier stokes with operator splitting
int n_geom_bases
number of geometric bases
assembler::AssemblyValsCache pressure_ass_vals_cache
used to store assembly values for pressure for small problems
std::vector< mesh::LocalBoundary > local_boundary
mapping from elements to nodes for dirichlet boundary conditions
double avg_mass
average system mass, used for contact with IPC
bool iso_parametric() const
check if using iso parametric bases
std::shared_ptr< assembler::RhsAssembler > build_rhs_assembler() const
build a RhsAssembler for the problem
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
void solve_transient_navier_stokes(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves transient navier stokes with FEM
void init_solve(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
initialize solver
void init_linear_solve(Eigen::MatrixXd &sol, const double t=1.0)
initialize the linear solve
assembler::AssemblyValsCache mass_ass_vals_cache
void solve_tensor_nonlinear(Eigen::MatrixXd &sol, const int t=0, const bool init_lagging=true)
solves nonlinear problems
std::string resolve_input_path(const std::string &path, const bool only_if_exists=false) const
Resolve input path relative to root_path() if the path is not absolute.
std::vector< int > boundary_nodes
list of boundary nodes
solver::SolveData solve_data
timedependent stuff cached
std::shared_ptr< polyfem::mesh::MeshNodes > geom_mesh_nodes
void solve_linear(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves a linear problem
void init_nonlinear_tensor_solve(Eigen::MatrixXd &sol, const double t=1.0, const bool init_time_integrator=true)
initialize the nonlinear solver
bool is_contact_enabled() const
does the simulation has contact
void build_collision_mesh()
extracts the boundary mesh for collision, called in build_basis
Eigen::VectorXi disc_orders
vector of discretization orders, used when not all elements have the same degree, one per element
std::vector< mesh::LocalBoundary > local_neumann_boundary
mapping from elements to nodes for neumann boundary conditions
std::shared_ptr< assembler::PressureAssembler > build_pressure_assembler() const
void solve_transient_linear(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves transient linear problem
bool is_problem_linear() const
Returns whether the system is linear. Collisions and pressure add nonlinearity to the problem.
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.
std::unordered_map< int, std::vector< mesh::LocalBoundary > > local_pressure_cavity
mapping from elements to nodes for pressure boundary conditions
std::string velocity() const
static double convert(const json &val, const std::string &unit_type)
const std::string & length() const
Caches basis evaluation and geometric mapping at every element.
void init(const bool is_volume, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const bool is_mass=false)
computes the basis evaluation and geometric mapping for each of the given ElementBases in bases initi...
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,...
quadrature::Quadrature quadrature
assemble matrix based on the local assembler local assembler is eg Laplace, LinearElasticity etc
void init(const int dim, const json ¶m)
static int build_bases(const mesh::Mesh2D &mesh, const std::string &assembler, const int quadrature_order, const int mass_quadrature_order, const int discr_order, const bool serendipity, const bool has_polys, const bool is_geom_bases, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, InterfaceData > &poly_edge_to_data, std::shared_ptr< mesh::MeshNodes > &mesh_nodes)
Builds FE basis functions over the entire mesh (P1, P2 over triangles, Q1, Q2 over quads).
static int build_bases(const mesh::Mesh3D &mesh, const std::string &assembler, const int quadrature_order, const int mass_quadrature_order, const int discr_order, const bool serendipity, const bool has_polys, const bool is_geom_bases, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, InterfaceData > &poly_face_to_data, std::shared_ptr< mesh::MeshNodes > &mesh_nodes)
Builds FE basis functions over the entire mesh (P1, P2 over tets, Q1, Q2 over hes).
static int build_bases(const std::string &assembler_name, const int dim, const mesh::Mesh2D &mesh, const int n_bases, const int quadrature_order, const int mass_quadrature_order, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, Eigen::MatrixXd > &mapped_boundary)
static int build_bases(const assembler::LinearAssembler &assembler, const int n_samples_per_edge, const mesh::Mesh2D &mesh, const int n_bases, const int quadrature_order, const int mass_quadrature_order, const int integral_constraints, std::vector< ElementBases > &bases, const std::vector< ElementBases > &gbases, const std::map< int, InterfaceData > &poly_edge_to_data, std::map< int, Eigen::MatrixXd > &mapped_boundary)
Build bases over the remaining polygons of a mesh.
static int build_bases(const assembler::LinearAssembler &assembler, const int n_samples_per_edge, const mesh::Mesh3D &mesh, const int n_bases, const int quadrature_order, const int mass_quadrature_order, const int integral_constraints, std::vector< ElementBases > &bases, const std::vector< ElementBases > &gbases, const std::map< int, InterfaceData > &poly_face_to_data, std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &mapped_boundary)
Build bases over the remaining polygons of a mesh.
static int build_bases(const mesh::Mesh2D &mesh, const std::string &assembler, const int quadrature_order, const int mass_quadrature_order, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, InterfaceData > &poly_edge_to_data)
static int build_bases(const mesh::Mesh3D &mesh, const std::string &assembler, const int quadrature_order, const int mass_quadrature_order, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, InterfaceData > &poly_face_to_data)
static int build_bases(const std::string &assembler_name, const int dim, const mesh::Mesh2D &mesh, const int n_bases, const int quadrature_order, const int mass_quadrature_order, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, Eigen::MatrixXd > &mapped_boundary)
void build_grid(const polyfem::mesh::Mesh &mesh, const double spacing)
builds the grid to export the solution
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
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
int n_flipped
number of flipped elements, compute only when using count_flipped_els (false by default)
Eigen::Vector4d spectrum
spectrum of the stiffness matrix, enable only if POLYSOLVE_WITH_SPECTRA is ON (off by default)
void count_flipped_elements(const polyfem::mesh::Mesh &mesh, const std::vector< polyfem::basis::ElementBases > &gbases)
counts the number of flipped elements
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...
long long nn_zero
non zeros and sytem matrix size num dof is the total dof in the system
double mesh_size
max edge lenght
void reset()
clears all stats
double min_edge_length
min edge lenght
bool is_volume() const override
checks if mesh is volume
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
const Eigen::MatrixXi & e() const
const Eigen::MatrixXi & f() const
const Eigen::MatrixXd & v() const
const Eigen::VectorXi & codim_v() const
static void p_refine(const mesh::Mesh &mesh, const double B, const bool h1_formula, const int base_p, const int discr_order_max, io::OutStatsData &stats, Eigen::VectorXi &disc_orders)
compute a priori prefinement
std::shared_ptr< assembler::RhsAssembler > rhs_assembler
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)
bool read_matrix(const std::string &path, Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &mat)
Reads a matrix from a file. Determines the file format based on the path's extension.
bool write_matrix(const std::string &path, const Mat &mat)
Writes a matrix to a file. Determines the file format based on the path's extension.
void load_collision_proxy(const std::string &mesh_filename, const std::string &weights_filename, const Eigen::VectorXi &in_node_to_node, const json &transformation, Eigen::MatrixXd &vertices, Eigen::VectorXi &codim_vertices, Eigen::MatrixXi &edges, Eigen::MatrixXi &faces, std::vector< Eigen::Triplet< double > > &displacement_map_entries)
Load a collision proxy mesh and displacement map from files.
void build_collision_proxy(const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &geom_bases, const std::vector< LocalBoundary > &total_local_boundary, const int n_bases, const int dim, const double max_edge_length, Eigen::MatrixXd &proxy_vertices, Eigen::MatrixXi &proxy_faces, std::vector< Eigen::Triplet< double > > &displacement_map_entries, const CollisionProxyTessellation tessellation)
Eigen::SparseMatrix< double > lump_matrix(const Eigen::SparseMatrix< double > &M)
Lump each row of a matrix into the diagonal.
std::string resolve_path(const std::string &path, const std::string &input_file_path, const bool only_if_exists=false)
std::vector< T > json_as_array(const json &j)
Return the value of a json object as an array.
void append_rows(DstMat &dst, const SrcMat &src)
spdlog::logger & logger()
Retrieves the current logger.
void compute_integral_constraints(const Mesh3D &mesh, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, Eigen::MatrixXd &basis_integrals)
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 > VectorNd
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
void log_and_throw_error(const std::string &msg)
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix