Loading [MathJax]/jax/input/TeX/config.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
State.cpp
Go to the documentation of this file.
1#include <polyfem/State.hpp>
2#include <polyfem/Common.hpp>
3
7
10
19
22
24
27
30
33
36
41
44
45#include <polysolve/linear/FEMSolver.hpp>
46
48
49#include <igl/edges.h>
50#include <igl/Timer.h>
51
52#include <iostream>
53#include <algorithm>
54#include <memory>
55#include <filesystem>
56
58
61
62using namespace Eigen;
63
64namespace polyfem
65{
66 using namespace assembler;
67 using namespace mesh;
68 using namespace io;
69 using namespace utils;
70
71 namespace
72 {
74 void build_in_node_to_in_primitive(const Mesh &mesh, const MeshNodes &mesh_nodes,
75 Eigen::VectorXi &in_node_to_in_primitive,
76 Eigen::VectorXi &in_node_offset)
77 {
78 const int num_vertex_nodes = mesh_nodes.num_vertex_nodes();
79 const int num_edge_nodes = mesh_nodes.num_edge_nodes();
80 const int num_face_nodes = mesh_nodes.num_face_nodes();
81 const int num_cell_nodes = mesh_nodes.num_cell_nodes();
82
83 const int num_nodes = num_vertex_nodes + num_edge_nodes + num_face_nodes + num_cell_nodes;
84
85 const long n_vertices = num_vertex_nodes;
86 const int num_in_primitives = n_vertices + mesh.n_edges() + mesh.n_faces() + mesh.n_cells();
87 const int num_primitives = mesh.n_vertices() + mesh.n_edges() + mesh.n_faces() + mesh.n_cells();
88
89 in_node_to_in_primitive.resize(num_nodes);
90 in_node_offset.resize(num_nodes);
91
92 // Only one node per vertex, so this is an identity map.
93 in_node_to_in_primitive.head(num_vertex_nodes).setLinSpaced(num_vertex_nodes, 0, num_vertex_nodes - 1); // vertex nodes
94 in_node_offset.head(num_vertex_nodes).setZero();
95
96 int prim_offset = n_vertices;
97 int node_offset = num_vertex_nodes;
98 auto foo = [&](const int num_prims, const int num_prim_nodes) {
99 if (num_prims <= 0 || num_prim_nodes <= 0)
100 return;
101 const Eigen::VectorXi range = Eigen::VectorXi::LinSpaced(num_prim_nodes, 0, num_prim_nodes - 1);
102 // TODO: This assumes isotropic degree of element.
103 const int node_per_prim = num_prim_nodes / num_prims;
104
105 in_node_to_in_primitive.segment(node_offset, num_prim_nodes) =
106 range.array() / node_per_prim + prim_offset;
107
108 in_node_offset.segment(node_offset, num_prim_nodes) =
109 range.unaryExpr([&](const int x) { return x % node_per_prim; });
110
111 prim_offset += num_prims;
112 node_offset += num_prim_nodes;
113 };
114
115 foo(mesh.n_edges(), num_edge_nodes);
116 foo(mesh.n_faces(), num_face_nodes);
117 foo(mesh.n_cells(), num_cell_nodes);
118 }
119
120 bool build_in_primitive_to_primitive(
121 const Mesh &mesh, const MeshNodes &mesh_nodes,
122 const Eigen::VectorXi &in_ordered_vertices,
123 const Eigen::MatrixXi &in_ordered_edges,
124 const Eigen::MatrixXi &in_ordered_faces,
125 Eigen::VectorXi &in_primitive_to_primitive)
126 {
127 // NOTE: Assume in_cells_to_cells is identity
128 const int num_vertex_nodes = mesh_nodes.num_vertex_nodes();
129 const int num_edge_nodes = mesh_nodes.num_edge_nodes();
130 const int num_face_nodes = mesh_nodes.num_face_nodes();
131 const int num_cell_nodes = mesh_nodes.num_cell_nodes();
132 const int num_nodes = num_vertex_nodes + num_edge_nodes + num_face_nodes + num_cell_nodes;
133
134 const long n_vertices = num_vertex_nodes;
135 const int num_in_primitives = n_vertices + mesh.n_edges() + mesh.n_faces() + mesh.n_cells();
136 const int num_primitives = mesh.n_vertices() + mesh.n_edges() + mesh.n_faces() + mesh.n_cells();
137
138 in_primitive_to_primitive.setLinSpaced(num_in_primitives, 0, num_in_primitives - 1);
139
140 igl::Timer timer;
141
142 // ------------
143 // Map vertices
144 // ------------
145
146 if (in_ordered_vertices.rows() != n_vertices)
147 {
148 logger().warn("Node ordering disabled, in_ordered_vertices != n_vertices, {} != {}", in_ordered_vertices.rows(), n_vertices);
149 return false;
150 }
151
152 in_primitive_to_primitive.head(n_vertices) = in_ordered_vertices;
153
154 int in_offset = n_vertices;
155 int offset = mesh.n_vertices();
156
157 // ---------
158 // Map edges
159 // ---------
160
161 logger().trace("Building Mesh edges to IDs...");
162 timer.start();
163 const auto edges_to_ids = mesh.edges_to_ids();
164 if (in_ordered_edges.rows() != edges_to_ids.size())
165 {
166 logger().warn("Node ordering disabled, in_ordered_edges != edges_to_ids, {} != {}", in_ordered_edges.rows(), edges_to_ids.size());
167 return false;
168 }
169 timer.stop();
170 logger().trace("Done (took {}s)", timer.getElapsedTime());
171
172 logger().trace("Building in-edge to edge mapping...");
173 timer.start();
174 for (int in_ei = 0; in_ei < in_ordered_edges.rows(); in_ei++)
175 {
176 const std::pair<int, int> in_edge(
177 in_ordered_edges.row(in_ei).minCoeff(),
178 in_ordered_edges.row(in_ei).maxCoeff());
179 in_primitive_to_primitive[in_offset + in_ei] =
180 offset + edges_to_ids.at(in_edge); // offset edge ids
181 }
182 timer.stop();
183 logger().trace("Done (took {}s)", timer.getElapsedTime());
184
185 in_offset += mesh.n_edges();
186 offset += mesh.n_edges();
187
188 // ---------
189 // Map faces
190 // ---------
191
192 if (mesh.is_volume())
193 {
194 logger().trace("Building Mesh faces to IDs...");
195 timer.start();
196 const auto faces_to_ids = mesh.faces_to_ids();
197 if (in_ordered_faces.rows() != faces_to_ids.size())
198 {
199 logger().warn("Node ordering disabled, in_ordered_faces != faces_to_ids, {} != {}", in_ordered_faces.rows(), faces_to_ids.size());
200 return false;
201 }
202 timer.stop();
203 logger().trace("Done (took {}s)", timer.getElapsedTime());
204
205 logger().trace("Building in-face to face mapping...");
206 timer.start();
207 for (int in_fi = 0; in_fi < in_ordered_faces.rows(); in_fi++)
208 {
209 std::vector<int> in_face(in_ordered_faces.cols());
210 for (int i = 0; i < in_face.size(); i++)
211 in_face[i] = in_ordered_faces(in_fi, i);
212 std::sort(in_face.begin(), in_face.end());
213
214 in_primitive_to_primitive[in_offset + in_fi] =
215 offset + faces_to_ids.at(in_face); // offset face ids
216 }
217 timer.stop();
218 logger().trace("Done (took {}s)", timer.getElapsedTime());
219
220 in_offset += mesh.n_faces();
221 offset += mesh.n_faces();
222 }
223
224 return true;
225 }
226 } // namespace
227
228 std::vector<int> State::primitive_to_node() const
229 {
230 auto indices = iso_parametric() ? mesh_nodes->primitive_to_node() : geom_mesh_nodes->primitive_to_node();
231 indices.resize(mesh->n_vertices());
232 return indices;
233 }
234
235 std::vector<int> State::node_to_primitive() const
236 {
237 auto p2n = primitive_to_node();
238 std::vector<int> indices;
239 indices.resize(n_geom_bases);
240 for (int i = 0; i < p2n.size(); i++)
241 indices[p2n[i]] = i;
242 return indices;
243 }
244
246 {
247 if (args["space"]["basis_type"] == "Spline")
248 {
249 logger().warn("Node ordering disabled, it dosent work for splines!");
250 return;
251 }
252
253 if (disc_orders.maxCoeff() >= 4 || disc_orders.maxCoeff() != disc_orders.minCoeff())
254 {
255 logger().warn("Node ordering disabled, it works only for p < 4 and uniform order!");
256 return;
257 }
258
259 if (!mesh->is_conforming())
260 {
261 logger().warn("Node ordering disabled, not supported for non-conforming meshes!");
262 return;
263 }
264
265 if (mesh->has_poly())
266 {
267 logger().warn("Node ordering disabled, not supported for polygonal meshes!");
268 return;
269 }
270
271 if (mesh->in_ordered_vertices().size() <= 0 || mesh->in_ordered_edges().size() <= 0 || (mesh->is_volume() && mesh->in_ordered_faces().size() <= 0))
272 {
273 logger().warn("Node ordering disabled, input vertices/edges/faces not computed!");
274 return;
275 }
276
277 const int num_vertex_nodes = mesh_nodes->num_vertex_nodes();
278 const int num_edge_nodes = mesh_nodes->num_edge_nodes();
279 const int num_face_nodes = mesh_nodes->num_face_nodes();
280 const int num_cell_nodes = mesh_nodes->num_cell_nodes();
281
282 const int num_nodes = num_vertex_nodes + num_edge_nodes + num_face_nodes + num_cell_nodes;
283
284 const long n_vertices = num_vertex_nodes;
285 const int num_in_primitives = n_vertices + mesh->n_edges() + mesh->n_faces() + mesh->n_cells();
286 const int num_primitives = mesh->n_vertices() + mesh->n_edges() + mesh->n_faces() + mesh->n_cells();
287
288 igl::Timer timer;
289
290 logger().trace("Building in-node to in-primitive mapping...");
291 timer.start();
292 Eigen::VectorXi in_node_to_in_primitive;
293 Eigen::VectorXi in_node_offset;
294 build_in_node_to_in_primitive(*mesh, *mesh_nodes, in_node_to_in_primitive, in_node_offset);
295 timer.stop();
296 logger().trace("Done (took {}s)", timer.getElapsedTime());
297
298 logger().trace("Building in-primitive to primitive mapping...");
299 timer.start();
300 bool ok = build_in_primitive_to_primitive(
301 *mesh, *mesh_nodes,
302 mesh->in_ordered_vertices(),
303 mesh->in_ordered_edges(),
304 mesh->in_ordered_faces(),
306 timer.stop();
307 logger().trace("Done (took {}s)", timer.getElapsedTime());
308
309 if (!ok)
310 {
311 in_node_to_node.resize(0);
313 return;
314 }
315 const auto &tmp = mesh_nodes->in_ordered_vertices();
316 int max_tmp = -1;
317 for (auto v : tmp)
318 max_tmp = std::max(max_tmp, v);
319
320 in_node_to_node.resize(max_tmp + 1);
321 for (int i = 0; i < tmp.size(); ++i)
322 {
323 if (tmp[i] >= 0)
324 in_node_to_node[tmp[i]] = i;
325 }
326 }
327
328 std::string State::formulation() const
329 {
330 if (args["materials"].is_null())
331 {
332 logger().error("specify some 'materials'");
333 assert(!args["materials"].is_null());
334 throw std::runtime_error("invalid input");
335 }
336
337 if (args["materials"].is_array())
338 {
339 std::string current = "";
340 for (const auto &m : args["materials"])
341 {
342 const std::string tmp = m["type"];
343 if (current.empty())
344 current = tmp;
345 else if (current != tmp)
346 {
348 {
350 current = "MultiModels";
351 else
352 {
353 logger().error("Current material is {}, new material is {}, multimaterial supported only for LinearElasticity and NeoHookean", current, tmp);
354 throw std::runtime_error("invalid input");
355 }
356 }
357 else
358 {
359 logger().error("Current material is {}, new material is {}, multimaterial supported only for LinearElasticity and NeoHookean", current, tmp);
360 throw std::runtime_error("invalid input");
361 }
362 }
363 }
364
365 return current;
366 }
367 else
368 return args["materials"]["type"];
369 }
370
371 void State::sol_to_pressure(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
372 {
373 if (n_pressure_bases <= 0)
374 {
375 logger().error("No pressure bases defined!");
376 return;
377 }
378
379 assert(mixed_assembler != nullptr);
380 Eigen::MatrixXd tmp = sol;
381
382 int fluid_offset = use_avg_pressure ? (assembler->is_fluid() ? 1 : 0) : 0;
383 sol = tmp.topRows(tmp.rows() - n_pressure_bases - fluid_offset);
384 assert(sol.size() == n_bases * (problem->is_scalar() ? 1 : mesh->dimension()));
385 pressure = tmp.middleRows(tmp.rows() - n_pressure_bases - fluid_offset, n_pressure_bases);
386 assert(pressure.size() == n_pressure_bases);
387 }
388
390 const Mesh3D &mesh,
391 const int n_bases,
392 const std::vector<basis::ElementBases> &bases,
393 const std::vector<basis::ElementBases> &gbases,
394 Eigen::MatrixXd &basis_integrals)
395 {
396 if (!mesh.is_volume())
397 {
398 logger().error("Works only on volumetric meshes!");
399 return;
400 }
401 assert(mesh.is_volume());
402
403 basis_integrals.resize(n_bases, 9);
404 basis_integrals.setZero();
405 Eigen::MatrixXd rhs(n_bases, 9);
406 rhs.setZero();
407
408 const int n_elements = mesh.n_elements();
409 for (int e = 0; e < n_elements; ++e)
410 {
411 // if (mesh.is_polytope(e)) {
412 // continue;
413 // }
414 // ElementAssemblyValues vals = values[e];
415 // const ElementAssemblyValues &gvals = gvalues[e];
417 vals.compute(e, mesh.is_volume(), bases[e], gbases[e]);
418
419 // Computes the discretized integral of the PDE over the element
420 const int n_local_bases = int(vals.basis_values.size());
421 for (int j = 0; j < n_local_bases; ++j)
422 {
423 const AssemblyValues &v = vals.basis_values[j];
424 const double integral_100 = (v.grad_t_m.col(0).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
425 const double integral_010 = (v.grad_t_m.col(1).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
426 const double integral_001 = (v.grad_t_m.col(2).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
427
428 const double integral_110 = ((vals.val.col(1).array() * v.grad_t_m.col(0).array() + vals.val.col(0).array() * v.grad_t_m.col(1).array()) * vals.det.array() * vals.quadrature.weights.array()).sum();
429 const double integral_011 = ((vals.val.col(2).array() * v.grad_t_m.col(1).array() + vals.val.col(1).array() * v.grad_t_m.col(2).array()) * vals.det.array() * vals.quadrature.weights.array()).sum();
430 const double integral_101 = ((vals.val.col(0).array() * v.grad_t_m.col(2).array() + vals.val.col(2).array() * v.grad_t_m.col(0).array()) * vals.det.array() * vals.quadrature.weights.array()).sum();
431
432 const double integral_200 = 2 * (vals.val.col(0).array() * v.grad_t_m.col(0).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
433 const double integral_020 = 2 * (vals.val.col(1).array() * v.grad_t_m.col(1).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
434 const double integral_002 = 2 * (vals.val.col(2).array() * v.grad_t_m.col(2).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
435
436 const double area = (v.val.array() * vals.det.array() * vals.quadrature.weights.array()).sum();
437
438 for (size_t ii = 0; ii < v.global.size(); ++ii)
439 {
440 basis_integrals(v.global[ii].index, 0) += integral_100 * v.global[ii].val;
441 basis_integrals(v.global[ii].index, 1) += integral_010 * v.global[ii].val;
442 basis_integrals(v.global[ii].index, 2) += integral_001 * v.global[ii].val;
443
444 basis_integrals(v.global[ii].index, 3) += integral_110 * v.global[ii].val;
445 basis_integrals(v.global[ii].index, 4) += integral_011 * v.global[ii].val;
446 basis_integrals(v.global[ii].index, 5) += integral_101 * v.global[ii].val;
447
448 basis_integrals(v.global[ii].index, 6) += integral_200 * v.global[ii].val;
449 basis_integrals(v.global[ii].index, 7) += integral_020 * v.global[ii].val;
450 basis_integrals(v.global[ii].index, 8) += integral_002 * v.global[ii].val;
451
452 rhs(v.global[ii].index, 6) += -2.0 * area * v.global[ii].val;
453 rhs(v.global[ii].index, 7) += -2.0 * area * v.global[ii].val;
454 rhs(v.global[ii].index, 8) += -2.0 * area * v.global[ii].val;
455 }
456 }
457 }
458
459 basis_integrals -= rhs;
460 }
461
463 {
464 if (mesh->has_poly())
465 return true;
466
467 if (args["space"]["basis_type"] == "Bernstein")
468 return false;
469
470 if (args["space"]["basis_type"] == "Spline")
471 return true;
472
473 if (mesh->is_rational())
474 return false;
475
476 if (args["space"]["use_p_ref"])
477 return false;
478
479 if (has_periodic_bc())
480 return false;
481
483 return false;
484
485 if (mesh->orders().size() <= 0)
486 {
487 if (args["space"]["discr_order"] == 1)
488 return true;
489 else
490 return args["space"]["advanced"]["isoparametric"];
491 }
492
493 if (mesh->orders().minCoeff() != mesh->orders().maxCoeff())
494 return false;
495
496 if (args["space"]["discr_order"] == mesh->orders().minCoeff())
497 return true;
498
499 // TODO:
500 // if (args["space"]["discr_order"] == 1 && args["force_linear_geometry"])
501 // return true;
502
503 return args["space"]["advanced"]["isoparametric"];
504 }
505
507 {
508 if (!mesh)
509 {
510 logger().error("Load the mesh first!");
511 return;
512 }
513
514 mesh->prepare_mesh();
515
516 bases.clear();
517 pressure_bases.clear();
518 geom_bases_.clear();
519 boundary_nodes.clear();
520 dirichlet_nodes.clear();
521 neumann_nodes.clear();
522 local_boundary.clear();
523 total_local_boundary.clear();
526 local_pressure_cavity.clear();
527 polys.clear();
528 poly_edge_to_data.clear();
529 rhs.resize(0, 0);
530 basis_nodes_to_gbasis_nodes.resize(0, 0);
531
532 if (assembler::MultiModel *mm = dynamic_cast<assembler::MultiModel *>(assembler.get()))
533 {
534 assert(args["materials"].is_array());
535
536 std::vector<std::string> materials(mesh->n_elements());
537
538 std::map<int, std::string> mats;
539
540 for (const auto &m : args["materials"])
541 mats[m["id"].get<int>()] = m["type"];
542
543 for (int i = 0; i < materials.size(); ++i)
544 materials[i] = mats.at(mesh->get_body_id(i));
545
546 mm->init_multimodels(materials);
547 }
548
549 n_bases = 0;
550 n_geom_bases = 0;
552
553 stats.reset();
554
555 disc_orders.resize(mesh->n_elements());
556
557 problem->init(*mesh);
558 logger().info("Building {} basis...", (iso_parametric() ? "isoparametric" : "not isoparametric"));
559 const bool has_polys = mesh->has_poly();
560
561 local_boundary.clear();
564 local_pressure_cavity.clear();
565 std::map<int, basis::InterfaceData> poly_edge_to_data_geom; // temp dummy variable
566
567 const auto &tmp_json = args["space"]["discr_order"];
568 if (tmp_json.is_number_integer())
569 {
570 disc_orders.setConstant(tmp_json);
571 }
572 else if (tmp_json.is_string())
573 {
574 const std::string discr_orders_path = utils::resolve_path(tmp_json, root_path());
575 Eigen::MatrixXi tmp;
576 read_matrix(discr_orders_path, tmp);
577 assert(tmp.size() == disc_orders.size());
578 assert(tmp.cols() == 1);
579 disc_orders = tmp;
580 }
581 else if (tmp_json.is_array())
582 {
583 const auto b_discr_orders = tmp_json;
584
585 std::map<int, int> b_orders;
586 for (size_t i = 0; i < b_discr_orders.size(); ++i)
587 {
588 assert(b_discr_orders[i]["id"].is_array() || b_discr_orders[i]["id"].is_number_integer());
589
590 const int order = b_discr_orders[i]["order"];
591 for (const int id : json_as_array<int>(b_discr_orders[i]["id"]))
592 {
593 b_orders[id] = order;
594 logger().trace("bid {}, discr {}", id, order);
595 }
596 }
597
598 for (int e = 0; e < mesh->n_elements(); ++e)
599 {
600 const int bid = mesh->get_body_id(e);
601 const auto order = b_orders.find(bid);
602 if (order == b_orders.end())
603 {
604 logger().debug("Missing discretization order for body {}; using 1", bid);
605 b_orders[bid] = 1;
606 disc_orders[e] = 1;
607 }
608 else
609 {
610 disc_orders[e] = order->second;
611 }
612 }
613 }
614 else
615 {
616 logger().error("space/discr_order must be either a number a path or an array");
617 throw std::runtime_error("invalid json");
618 }
619 // TODO: same for pressure!
620
621#ifdef POLYFEM_WITH_BEZIER
622 if (!mesh->is_simplicial())
623#else
624 if constexpr (true)
625#endif
626 {
627 args["space"]["advanced"]["count_flipped_els_continuous"] = false;
628 args["output"]["paraview"]["options"]["jacobian_validity"] = false;
629 args["solver"]["advanced"]["check_inversion"] = "Discrete";
630 }
631 else if (args["solver"]["advanced"]["check_inversion"] != "Discrete")
632 {
633 args["space"]["advanced"]["use_corner_quadrature"] = true;
634 }
635
636 Eigen::MatrixXi geom_disc_orders;
637 if (!iso_parametric())
638 {
639 if (mesh->orders().size() <= 0)
640 {
641 geom_disc_orders.resizeLike(disc_orders);
642 geom_disc_orders.setConstant(1);
643 }
644 else
645 geom_disc_orders = mesh->orders();
646 }
647
648 igl::Timer timer;
649 timer.start();
650 if (args["space"]["use_p_ref"])
651 {
653 *mesh,
654 args["space"]["advanced"]["B"],
655 args["space"]["advanced"]["h1_formula"],
656 args["space"]["discr_order"],
657 args["space"]["advanced"]["discr_order_max"],
658 stats,
660
661 logger().info("min p: {} max p: {}", disc_orders.minCoeff(), disc_orders.maxCoeff());
662 }
663
664 int quadrature_order = args["space"]["advanced"]["quadrature_order"].get<int>();
665 const int mass_quadrature_order = args["space"]["advanced"]["mass_quadrature_order"].get<int>();
666 if (mixed_assembler != nullptr)
667 {
668 const int disc_order = disc_orders.maxCoeff();
669 if (disc_order - disc_orders.minCoeff() != 0)
670 {
671 logger().error("p refinement not supported in mixed formulation!");
672 return;
673 }
674 }
675
676 // shape optimization needs continuous geometric basis
677 // const bool use_continuous_gbasis = optimization_enabled == solver::CacheLevel::Derivatives;
678 const bool use_continuous_gbasis = true;
679 const bool use_corner_quadrature = args["space"]["advanced"]["use_corner_quadrature"];
680
681 if (mesh->is_volume())
682 {
683 const Mesh3D &tmp_mesh = *dynamic_cast<Mesh3D *>(mesh.get());
684 if (args["space"]["basis_type"] == "Spline")
685 {
686 // if (!iso_parametric())
687 // {
688 // logger().error("Splines must be isoparametric, ignoring...");
689 // // LagrangeBasis3d::build_bases(tmp_mesh, quadrature_order, geom_disc_orders, has_polys, geom_bases_, local_boundary, poly_edge_to_data_geom, mesh_nodes);
690 // SplineBasis3d::build_bases(tmp_mesh, quadrature_order, geom_bases_, local_boundary, poly_edge_to_data);
691 // }
692
693 n_bases = basis::SplineBasis3d::build_bases(tmp_mesh, assembler->name(), quadrature_order, mass_quadrature_order, bases, local_boundary, poly_edge_to_data);
694
695 // if (iso_parametric() && args["fit_nodes"])
696 // SplineBasis3d::fit_nodes(tmp_mesh, n_bases, bases);
697 }
698 else
699 {
700 if (!iso_parametric())
701 n_geom_bases = basis::LagrangeBasis3d::build_bases(tmp_mesh, assembler->name(), quadrature_order, mass_quadrature_order, geom_disc_orders, false, false, has_polys, !use_continuous_gbasis, use_corner_quadrature, geom_bases_, local_boundary, poly_edge_to_data_geom, geom_mesh_nodes);
702
703 n_bases = basis::LagrangeBasis3d::build_bases(tmp_mesh, assembler->name(), quadrature_order, mass_quadrature_order, disc_orders, args["space"]["basis_type"] == "Bernstein", args["space"]["basis_type"] == "Serendipity", has_polys, false, use_corner_quadrature, bases, local_boundary, poly_edge_to_data, mesh_nodes);
704 }
705
706 // if(problem->is_mixed())
707 if (mixed_assembler != nullptr)
708 {
709 n_pressure_bases = basis::LagrangeBasis3d::build_bases(tmp_mesh, assembler->name(), quadrature_order, mass_quadrature_order, int(args["space"]["pressure_discr_order"]), args["space"]["basis_type"] == "Bernstein", false, has_polys, false, use_corner_quadrature, pressure_bases, local_boundary, poly_edge_to_data_geom, pressure_mesh_nodes);
710 }
711 }
712 else
713 {
714 const Mesh2D &tmp_mesh = *dynamic_cast<Mesh2D *>(mesh.get());
715 if (args["space"]["basis_type"] == "Spline")
716 {
717 // TODO:
718 // if (!iso_parametric())
719 // {
720 // logger().error("Splines must be isoparametric, ignoring...");
721 // // LagrangeBasis2d::build_bases(tmp_mesh, quadrature_order, disc_orders, has_polys, geom_bases_, local_boundary, poly_edge_to_data_geom, mesh_nodes);
722 // n_bases = SplineBasis2d::build_bases(tmp_mesh, quadrature_order, geom_bases_, local_boundary, poly_edge_to_data);
723 // }
724
725 n_bases = basis::SplineBasis2d::build_bases(tmp_mesh, assembler->name(), quadrature_order, mass_quadrature_order, bases, local_boundary, poly_edge_to_data);
726
727 // if (iso_parametric() && args["fit_nodes"])
728 // SplineBasis2d::fit_nodes(tmp_mesh, n_bases, bases);
729 }
730 else
731 {
732 if (!iso_parametric())
733 n_geom_bases = basis::LagrangeBasis2d::build_bases(tmp_mesh, assembler->name(), quadrature_order, mass_quadrature_order, geom_disc_orders, false, false, has_polys, !use_continuous_gbasis, use_corner_quadrature, geom_bases_, local_boundary, poly_edge_to_data_geom, geom_mesh_nodes);
734
735 n_bases = basis::LagrangeBasis2d::build_bases(tmp_mesh, assembler->name(), quadrature_order, mass_quadrature_order, disc_orders, args["space"]["basis_type"] == "Bernstein", args["space"]["basis_type"] == "Serendipity", has_polys, false, use_corner_quadrature, bases, local_boundary, poly_edge_to_data, mesh_nodes);
736 }
737
738 // if(problem->is_mixed())
739 if (mixed_assembler != nullptr)
740 {
741 n_pressure_bases = basis::LagrangeBasis2d::build_bases(tmp_mesh, assembler->name(), quadrature_order, mass_quadrature_order, int(args["space"]["pressure_discr_order"]), args["space"]["basis_type"] == "Bernstein", false, has_polys, false, use_corner_quadrature, pressure_bases, local_boundary, poly_edge_to_data_geom, pressure_mesh_nodes);
742 }
743 }
744
745 if (mixed_assembler != nullptr)
746 {
747 assert(bases.size() == pressure_bases.size());
748 for (int i = 0; i < pressure_bases.size(); ++i)
749 {
751 bases[i].compute_quadrature(b_quad);
752 pressure_bases[i].set_quadrature([b_quad](quadrature::Quadrature &quad) { quad = b_quad; });
753 }
754 }
755
756 timer.stop();
757
759
760 if (n_geom_bases == 0)
762
763 auto &gbases = geom_bases();
764
766 {
767 std::map<std::array<int, 2>, double> pairs;
768 for (int e = 0; e < gbases.size(); e++)
769 {
770 const auto &gbs = gbases[e].bases;
771 const auto &bs = bases[e].bases;
772
773 Eigen::MatrixXd local_pts;
774 const int order = bs.front().order();
775 if (mesh->is_volume())
776 {
777 if (mesh->is_simplex(e))
778 autogen::p_nodes_3d(order, local_pts);
779 else
780 autogen::q_nodes_3d(order, local_pts);
781 }
782 else
783 {
784 if (mesh->is_simplex(e))
785 autogen::p_nodes_2d(order, local_pts);
786 else
787 autogen::q_nodes_2d(order, local_pts);
788 }
789
791 vals.compute(e, mesh->is_volume(), local_pts, gbases[e], gbases[e]);
792
793 for (int i = 0; i < bs.size(); i++)
794 {
795 for (int j = 0; j < gbs.size(); j++)
796 {
797 if (std::abs(vals.basis_values[j].val(i)) > 1e-7)
798 {
799 std::array<int, 2> index = {{gbs[j].global()[0].index, bs[i].global()[0].index}};
800 pairs.insert({index, vals.basis_values[j].val(i)});
801 }
802 }
803 }
804 }
805
806 const int dim = mesh->dimension();
807 std::vector<Eigen::Triplet<double>> coeffs;
808 coeffs.clear();
809 for (const auto &iter : pairs)
810 for (int d = 0; d < dim; d++)
811 coeffs.emplace_back(iter.first[0] * dim + d, iter.first[1] * dim + d, iter.second);
812
813 basis_nodes_to_gbasis_nodes.resize(n_geom_bases * mesh->dimension(), n_bases * mesh->dimension());
814 basis_nodes_to_gbasis_nodes.setFromTriplets(coeffs.begin(), coeffs.end());
815 }
816
817 for (const auto &lb : local_boundary)
818 total_local_boundary.emplace_back(lb);
819
820 const int dim = mesh->dimension();
821 const int problem_dim = problem->is_scalar() ? 1 : dim;
822
823 // handle periodic bc
824 if (has_periodic_bc())
825 {
826 // collect periodic directions
827 json directions = args["boundary_conditions"]["periodic_boundary"]["correspondence"];
828 Eigen::MatrixXd tile_offset = Eigen::MatrixXd::Identity(dim, dim);
829
830 if (directions.size() > 0)
831 {
832 Eigen::VectorXd tmp;
833 for (int d = 0; d < dim; d++)
834 {
835 tmp = directions[d];
836 if (tmp.size() != dim)
837 log_and_throw_error("Invalid size of periodic directions!");
838 tile_offset.col(d) = tmp;
839 }
840 }
841
842 periodic_bc = std::make_shared<PeriodicBoundary>(problem->is_scalar(), n_bases, bases, mesh_nodes, tile_offset, args["boundary_conditions"]["periodic_boundary"]["tolerance"].get<double>());
843
844 macro_strain_constraint.init(dim, args["boundary_conditions"]["periodic_boundary"]);
845 }
846
847 if (args["space"]["advanced"]["count_flipped_els"])
849
850 const int prev_bases = n_bases;
852
853 {
854 igl::Timer timer2;
855 logger().debug("Building node mapping...");
856 timer2.start();
858 problem->update_nodes(in_node_to_node);
859 mesh->update_nodes(in_node_to_node);
860 timer2.stop();
861 logger().debug("Done (took {}s)", timer2.getElapsedTime());
862 }
863
864 logger().info("Building collision mesh...");
866 if (periodic_bc && args["contact"]["periodic"])
868 logger().info("Done!");
869
870 const int prev_b_size = local_boundary.size();
871 problem->setup_bc(*mesh, n_bases - obstacle.n_vertices(),
880
881 // setp nodal values
882 {
884 for (int n = 0; n < dirichlet_nodes.size(); ++n)
885 {
886 const int n_id = dirichlet_nodes[n];
887 bool found = false;
888 for (const auto &bs : bases)
889 {
890 for (const auto &b : bs.bases)
891 {
892 for (const auto &lg : b.global())
893 {
894 if (lg.index == n_id)
895 {
896 dirichlet_nodes_position[n] = lg.node;
897 found = true;
898 break;
899 }
900 }
901
902 if (found)
903 break;
904 }
905
906 if (found)
907 break;
908 }
909
910 assert(found);
911 }
912
914 for (int n = 0; n < neumann_nodes.size(); ++n)
915 {
916 const int n_id = neumann_nodes[n];
917 bool found = false;
918 for (const auto &bs : bases)
919 {
920 for (const auto &b : bs.bases)
921 {
922 for (const auto &lg : b.global())
923 {
924 if (lg.index == n_id)
925 {
926 neumann_nodes_position[n] = lg.node;
927 found = true;
928 break;
929 }
930 }
931
932 if (found)
933 break;
934 }
935
936 if (found)
937 break;
938 }
939
940 assert(found);
941 }
942 }
943
944 const bool has_neumann = local_neumann_boundary.size() > 0 || local_boundary.size() < prev_b_size;
945 use_avg_pressure = !has_neumann;
946
947 for (int i = prev_bases; i < n_bases; ++i)
948 {
949 for (int d = 0; d < problem_dim; ++d)
950 boundary_nodes.push_back(i * problem_dim + d);
951 }
952
953 std::sort(boundary_nodes.begin(), boundary_nodes.end());
954 auto it = std::unique(boundary_nodes.begin(), boundary_nodes.end());
955 boundary_nodes.resize(std::distance(boundary_nodes.begin(), it));
956
957 // for elastic pure periodic problem, find an internal node and force zero dirichlet
958 if ((!problem->is_time_dependent() || args["time"]["quasistatic"]) && boundary_nodes.size() == 0 && has_periodic_bc())
959 {
960 // find an internal node to force zero dirichlet
961 std::vector<bool> isboundary(n_bases, false);
962 for (const auto &lb : total_local_boundary)
963 {
964 const int e = lb.element_id();
965 for (int i = 0; i < lb.size(); ++i)
966 {
967 const auto nodes = bases[e].local_nodes_for_primitive(lb.global_primitive_id(i), *mesh);
968
969 for (int n : nodes)
970 isboundary[bases[e].bases[n].global()[0].index] = true;
971 }
972 }
973 int i = 0;
974 for (; i < n_bases; i++)
975 if (!isboundary[i]) // (!periodic_bc->is_periodic_dof(i))
976 break;
977 if (i >= n_bases)
978 log_and_throw_error("Failed to find a non-periodic node!");
979 const int actual_dim = problem->is_scalar() ? 1 : mesh->dimension();
980 for (int d = 0; d < actual_dim; d++)
981 {
982 boundary_nodes.push_back(i * actual_dim + d);
983 }
984 logger().info("Fix solution at node {} to remove singularity due to periodic BC", i);
985 }
986
987 const auto &curret_bases = geom_bases();
988 const int n_samples = 10;
989 stats.compute_mesh_size(*mesh, curret_bases, n_samples, args["output"]["advanced"]["curved_mesh_size"]);
991 {
993 }
995 {
997 }
998
999 if (is_contact_enabled())
1000 {
1001 min_boundary_edge_length = std::numeric_limits<double>::max();
1002 for (const auto &edge : collision_mesh.edges().rowwise())
1003 {
1004 const VectorNd v0 = collision_mesh.rest_positions().row(edge(0));
1005 const VectorNd v1 = collision_mesh.rest_positions().row(edge(1));
1006 min_boundary_edge_length = std::min(min_boundary_edge_length, (v1 - v0).norm());
1007 }
1008
1009 double dhat = Units::convert(args["contact"]["dhat"], units.length());
1010 args["contact"]["epsv"] = Units::convert(args["contact"]["epsv"], units.velocity());
1011 args["contact"]["dhat"] = dhat;
1012
1013 if (!has_dhat && dhat > min_boundary_edge_length)
1014 {
1015 args["contact"]["dhat"] = double(args["contact"]["dhat_percentage"]) * min_boundary_edge_length;
1016 logger().info("dhat set to {}", double(args["contact"]["dhat"]));
1017 }
1018 else
1019 {
1020 if (dhat > min_boundary_edge_length)
1021 logger().warn("dhat larger than min boundary edge, {} > {}", dhat, min_boundary_edge_length);
1022 }
1023 }
1024
1025 logger().info("n_bases {}", n_bases);
1026
1027 timings.building_basis_time = timer.getElapsedTime();
1028 logger().info(" took {}s", timings.building_basis_time);
1029
1030 logger().info("flipped elements {}", stats.n_flipped);
1031 logger().info("h: {}", stats.mesh_size);
1032 logger().info("n bases: {}", n_bases);
1033 logger().info("n pressure bases: {}", n_pressure_bases);
1034
1037 if (n_bases <= args["solver"]["advanced"]["cache_size"])
1038 {
1039 timer.start();
1040 logger().info("Building cache...");
1041 ass_vals_cache.init(mesh->is_volume(), bases, curret_bases);
1042 mass_ass_vals_cache.init(mesh->is_volume(), bases, curret_bases, true);
1043 if (mixed_assembler != nullptr)
1044 pressure_ass_vals_cache.init(mesh->is_volume(), pressure_bases, curret_bases);
1045
1046 logger().info(" took {}s", timer.getElapsedTime());
1047 }
1048
1049 out_geom.build_grid(*mesh, args["output"]["advanced"]["sol_on_grid"]);
1050
1051 if ((!problem->is_time_dependent() || args["time"]["quasistatic"]) && boundary_nodes.empty())
1052 {
1053 if (has_periodic_bc())
1054 logger().warn("(Quasi-)Static problem without Dirichlet nodes, will fix solution at one node to find a unique solution!");
1055 else
1056 {
1057 if (args["constraints"]["hard"].empty())
1058 log_and_throw_error("Static problem need to have some Dirichlet nodes!");
1059 else
1060 logger().warn("Relying on hard constraints to avoid infinite solutions");
1061 }
1062 }
1063 }
1064
1066 {
1067 if (!mesh)
1068 {
1069 logger().error("Load the mesh first!");
1070 return;
1071 }
1072
1073 rhs.resize(0, 0);
1074
1075 if (poly_edge_to_data.empty() && polys.empty())
1076 {
1078 return;
1079 }
1080
1081 igl::Timer timer;
1082 timer.start();
1083 logger().info("Computing polygonal basis...");
1084
1085 // std::sort(boundary_nodes.begin(), boundary_nodes.end());
1086
1087 // mixed not supports polygonal bases
1088 assert(n_pressure_bases == 0 || poly_edge_to_data.size() == 0);
1089
1090 int new_bases = 0;
1091
1092 if (iso_parametric())
1093 {
1094 if (mesh->is_volume())
1095 {
1096 if (args["space"]["poly_basis_type"] == "MeanValue" || args["space"]["poly_basis_type"] == "Wachspress")
1097 logger().error("Barycentric bases not supported in 3D");
1098 assert(assembler->is_linear());
1100 *dynamic_cast<LinearAssembler *>(assembler.get()),
1101 args["space"]["advanced"]["n_harmonic_samples"],
1102 *dynamic_cast<Mesh3D *>(mesh.get()),
1103 n_bases,
1104 args["space"]["advanced"]["quadrature_order"],
1105 args["space"]["advanced"]["mass_quadrature_order"],
1106 args["space"]["advanced"]["integral_constraints"],
1107 bases,
1108 bases,
1110 polys_3d);
1111 }
1112 else
1113 {
1114 if (args["space"]["poly_basis_type"] == "MeanValue")
1115 {
1117 assembler->name(),
1118 assembler->is_tensor() ? 2 : 1,
1119 *dynamic_cast<Mesh2D *>(mesh.get()),
1120 n_bases,
1121 args["space"]["advanced"]["quadrature_order"],
1122 args["space"]["advanced"]["mass_quadrature_order"],
1124 }
1125 else if (args["space"]["poly_basis_type"] == "Wachspress")
1126 {
1128 assembler->name(),
1129 assembler->is_tensor() ? 2 : 1,
1130 *dynamic_cast<Mesh2D *>(mesh.get()),
1131 n_bases,
1132 args["space"]["advanced"]["quadrature_order"],
1133 args["space"]["advanced"]["mass_quadrature_order"],
1135 }
1136 else
1137 {
1138 assert(assembler->is_linear());
1140 *dynamic_cast<LinearAssembler *>(assembler.get()),
1141 args["space"]["advanced"]["n_harmonic_samples"],
1142 *dynamic_cast<Mesh2D *>(mesh.get()),
1143 n_bases,
1144 args["space"]["advanced"]["quadrature_order"],
1145 args["space"]["advanced"]["mass_quadrature_order"],
1146 args["space"]["advanced"]["integral_constraints"],
1147 bases,
1148 bases,
1150 polys);
1151 }
1152 }
1153 }
1154 else
1155 {
1156 if (mesh->is_volume())
1157 {
1158 if (args["space"]["poly_basis_type"] == "MeanValue" || args["space"]["poly_basis_type"] == "Wachspress")
1159 {
1160 logger().error("Barycentric bases not supported in 3D");
1161 throw std::runtime_error("not implemented");
1162 }
1163 else
1164 {
1165 assert(assembler->is_linear());
1167 *dynamic_cast<LinearAssembler *>(assembler.get()),
1168 args["space"]["advanced"]["n_harmonic_samples"],
1169 *dynamic_cast<Mesh3D *>(mesh.get()),
1170 n_bases,
1171 args["space"]["advanced"]["quadrature_order"],
1172 args["space"]["advanced"]["mass_quadrature_order"],
1173 args["space"]["advanced"]["integral_constraints"],
1174 bases,
1177 polys_3d);
1178 }
1179 }
1180 else
1181 {
1182 if (args["space"]["poly_basis_type"] == "MeanValue")
1183 {
1185 assembler->name(),
1186 assembler->is_tensor() ? 2 : 1,
1187 *dynamic_cast<Mesh2D *>(mesh.get()),
1188 n_bases, args["space"]["advanced"]["quadrature_order"],
1189 args["space"]["advanced"]["mass_quadrature_order"],
1191 }
1192 else if (args["space"]["poly_basis_type"] == "Wachspress")
1193 {
1195 assembler->name(),
1196 assembler->is_tensor() ? 2 : 1,
1197 *dynamic_cast<Mesh2D *>(mesh.get()),
1198 n_bases, args["space"]["advanced"]["quadrature_order"],
1199 args["space"]["advanced"]["mass_quadrature_order"],
1201 }
1202 else
1203 {
1204 assert(assembler->is_linear());
1206 *dynamic_cast<LinearAssembler *>(assembler.get()),
1207 args["space"]["advanced"]["n_harmonic_samples"],
1208 *dynamic_cast<Mesh2D *>(mesh.get()),
1209 n_bases,
1210 args["space"]["advanced"]["quadrature_order"],
1211 args["space"]["advanced"]["mass_quadrature_order"],
1212 args["space"]["advanced"]["integral_constraints"],
1213 bases,
1216 polys);
1217 }
1218 }
1219 }
1220
1221 timer.stop();
1222 timings.computing_poly_basis_time = timer.getElapsedTime();
1223 logger().info(" took {}s", timings.computing_poly_basis_time);
1224
1225 n_bases += new_bases;
1226 }
1227
1229 {
1230 assert(!mesh->is_volume());
1231 const int dim = mesh->dimension();
1232 const int n_tiles = 2;
1233
1234 if (mesh->dimension() != 2)
1235 log_and_throw_error("Periodic collision mesh is only implemented in 2D!");
1236
1237 Eigen::MatrixXd V(n_bases, dim);
1238 for (const auto &bs : bases)
1239 for (const auto &b : bs.bases)
1240 for (const auto &g : b.global())
1241 V.row(g.index) = g.node;
1242
1243 Eigen::MatrixXi E = collision_mesh.edges();
1244 for (int i = 0; i < E.size(); i++)
1245 E(i) = collision_mesh.to_full_vertex_id(E(i));
1246
1247 Eigen::MatrixXd bbox(V.cols(), 2);
1248 bbox.col(0) = V.colwise().minCoeff();
1249 bbox.col(1) = V.colwise().maxCoeff();
1250
1251 // remove boundary edges on periodic BC, buggy
1252 {
1253 std::vector<int> ind;
1254 for (int i = 0; i < E.rows(); i++)
1255 {
1256 if (!periodic_bc->is_periodic_dof(E(i, 0)) || !periodic_bc->is_periodic_dof(E(i, 1)))
1257 ind.push_back(i);
1258 }
1259
1260 E = E(ind, Eigen::all).eval();
1261 }
1262
1263 Eigen::MatrixXd Vtmp, Vnew;
1264 Eigen::MatrixXi Etmp, Enew;
1265 Vtmp.setZero(V.rows() * n_tiles * n_tiles, V.cols());
1266 Etmp.setZero(E.rows() * n_tiles * n_tiles, E.cols());
1267
1268 Eigen::MatrixXd tile_offset = periodic_bc->get_affine_matrix();
1269
1270 for (int i = 0, idx = 0; i < n_tiles; i++)
1271 {
1272 for (int j = 0; j < n_tiles; j++)
1273 {
1274 Eigen::Vector2d block_id;
1275 block_id << i, j;
1276
1277 Vtmp.middleRows(idx * V.rows(), V.rows()) = V;
1278 // Vtmp.block(idx * V.rows(), 0, V.rows(), 1).array() += tile_offset(0) * i;
1279 // Vtmp.block(idx * V.rows(), 1, V.rows(), 1).array() += tile_offset(1) * j;
1280 for (int vid = 0; vid < V.rows(); vid++)
1281 Vtmp.block(idx * V.rows() + vid, 0, 1, 2) += (tile_offset * block_id).transpose();
1282
1283 Etmp.middleRows(idx * E.rows(), E.rows()) = E.array() + idx * V.rows();
1284 idx += 1;
1285 }
1286 }
1287
1288 // clean duplicated vertices
1289 Eigen::VectorXi indices;
1290 {
1291 std::vector<int> tmp;
1292 for (int i = 0; i < V.rows(); i++)
1293 {
1294 if (periodic_bc->is_periodic_dof(i))
1295 tmp.push_back(i);
1296 }
1297
1298 indices.resize(tmp.size() * n_tiles * n_tiles);
1299 for (int i = 0; i < n_tiles * n_tiles; i++)
1300 {
1301 indices.segment(i * tmp.size(), tmp.size()) = Eigen::Map<Eigen::VectorXi, Eigen::Unaligned>(tmp.data(), tmp.size());
1302 indices.segment(i * tmp.size(), tmp.size()).array() += i * V.rows();
1303 }
1304 }
1305
1306 Eigen::VectorXi potentially_duplicate_mask(Vtmp.rows());
1307 potentially_duplicate_mask.setZero();
1308 potentially_duplicate_mask(indices).array() = 1;
1309
1310 Eigen::MatrixXd candidates = Vtmp(indices, Eigen::all);
1311
1312 Eigen::VectorXi SVI;
1313 std::vector<int> SVJ;
1314 SVI.setConstant(Vtmp.rows(), -1);
1315 int id = 0;
1316 const double eps = (bbox.col(1) - bbox.col(0)).maxCoeff() * args["boundary_conditions"]["periodic_boundary"]["tolerance"].get<double>();
1317 for (int i = 0; i < Vtmp.rows(); i++)
1318 {
1319 if (SVI[i] < 0)
1320 {
1321 SVI[i] = id;
1322 SVJ.push_back(i);
1323 if (potentially_duplicate_mask(i))
1324 {
1325 Eigen::VectorXd diffs = (candidates.rowwise() - Vtmp.row(i)).rowwise().norm();
1326 for (int j = 0; j < diffs.size(); j++)
1327 if (diffs(j) < eps)
1328 SVI[indices[j]] = id;
1329 }
1330 id++;
1331 }
1332 }
1333
1334 Vnew = Vtmp(SVJ, Eigen::all);
1335
1336 Enew.resizeLike(Etmp);
1337 for (int d = 0; d < Etmp.cols(); d++)
1338 Enew.col(d) = SVI(Etmp.col(d));
1339
1340 std::vector<bool> is_on_surface = ipc::CollisionMesh::construct_is_on_surface(Vnew.rows(), Enew);
1341
1342 Eigen::MatrixXi boundary_triangles;
1343 Eigen::SparseMatrix<double> displacement_map;
1344 periodic_collision_mesh = ipc::CollisionMesh(is_on_surface,
1345 Vnew,
1346 Enew,
1347 boundary_triangles,
1348 displacement_map);
1349
1350 periodic_collision_mesh.init_area_jacobians();
1351
1352 periodic_collision_mesh_to_basis.setConstant(Vnew.rows(), -1);
1353 for (int i = 0; i < V.rows(); i++)
1354 for (int j = 0; j < n_tiles * n_tiles; j++)
1355 periodic_collision_mesh_to_basis(SVI[j * V.rows() + i]) = i;
1356
1357 if (periodic_collision_mesh_to_basis.maxCoeff() + 1 != V.rows())
1358 log_and_throw_error("Failed to tile mesh!");
1359 }
1360
1362 {
1365 args, [this](const std::string &p) { return resolve_input_path(p); },
1367 }
1368
1370 const mesh::Mesh &mesh,
1371 const int n_bases,
1372 const std::vector<basis::ElementBases> &bases,
1373 const std::vector<basis::ElementBases> &geom_bases,
1374 const std::vector<mesh::LocalBoundary> &total_local_boundary,
1375 const mesh::Obstacle &obstacle,
1376 const json &args,
1377 const std::function<std::string(const std::string &)> &resolve_input_path,
1378 const Eigen::VectorXi &in_node_to_node,
1379 ipc::CollisionMesh &collision_mesh)
1380 {
1381 Eigen::MatrixXd collision_vertices;
1382 Eigen::VectorXi collision_codim_vids;
1383 Eigen::MatrixXi collision_edges, collision_triangles;
1384 std::vector<Eigen::Triplet<double>> displacement_map_entries;
1385
1386 if (args.contains("/contact/collision_mesh"_json_pointer)
1387 && args.at("/contact/collision_mesh/enabled"_json_pointer).get<bool>())
1388 {
1389 const json collision_mesh_args = args.at("/contact/collision_mesh"_json_pointer);
1390 if (collision_mesh_args.contains("linear_map"))
1391 {
1392 assert(displacement_map_entries.empty());
1393 assert(collision_mesh_args.contains("mesh"));
1394 const std::string root_path = utils::json_value<std::string>(args, "root_path", "");
1395 // TODO: handle transformation per geometry
1396 const json transformation = json_as_array(args["geometry"])[0]["transformation"];
1398 utils::resolve_path(collision_mesh_args["mesh"], root_path),
1399 utils::resolve_path(collision_mesh_args["linear_map"], root_path),
1400 in_node_to_node, transformation, collision_vertices, collision_codim_vids,
1401 collision_edges, collision_triangles, displacement_map_entries);
1402 }
1403 else
1404 {
1405 assert(collision_mesh_args.contains("max_edge_length"));
1406 logger().debug(
1407 "Building collision proxy with max edge length={} ...",
1408 collision_mesh_args["max_edge_length"].get<double>());
1409 igl::Timer timer;
1410 timer.start();
1413 collision_mesh_args["max_edge_length"], collision_vertices,
1414 collision_triangles, displacement_map_entries,
1415 collision_mesh_args["tessellation_type"]);
1416 if (collision_triangles.size())
1417 igl::edges(collision_triangles, collision_edges);
1418 timer.stop();
1419 logger().debug(fmt::format(
1420 std::locale("en_US.UTF-8"),
1421 "Done (took {:g}s, {:L} vertices, {:L} triangles)",
1422 timer.getElapsedTime(),
1423 collision_vertices.rows(), collision_triangles.rows()));
1424 }
1425 }
1426 else
1427 {
1430 collision_vertices, collision_edges, collision_triangles, displacement_map_entries);
1431 }
1432
1433 // n_bases already contains the obstacle vertices
1434 const int num_fe_nodes = n_bases - obstacle.n_vertices();
1435 const int num_fe_collision_vertices = collision_vertices.rows();
1436 assert(collision_edges.size() == 0 || collision_edges.maxCoeff() < num_fe_collision_vertices);
1437 assert(collision_triangles.size() == 0 || collision_triangles.maxCoeff() < num_fe_collision_vertices);
1438
1439 // Append the obstacles to the collision mesh
1440 if (obstacle.n_vertices() > 0)
1441 {
1442 append_rows(collision_vertices, obstacle.v());
1443 append_rows(collision_codim_vids, obstacle.codim_v().array() + num_fe_collision_vertices);
1444 append_rows(collision_edges, obstacle.e().array() + num_fe_collision_vertices);
1445 append_rows(collision_triangles, obstacle.f().array() + num_fe_collision_vertices);
1446
1447 if (!displacement_map_entries.empty())
1448 {
1449 displacement_map_entries.reserve(displacement_map_entries.size() + obstacle.n_vertices());
1450 for (int i = 0; i < obstacle.n_vertices(); i++)
1451 {
1452 displacement_map_entries.emplace_back(num_fe_collision_vertices + i, num_fe_nodes + i, 1.0);
1453 }
1454 }
1455 }
1456
1457 std::vector<bool> is_on_surface = ipc::CollisionMesh::construct_is_on_surface(
1458 collision_vertices.rows(), collision_edges);
1459 for (const int vid : collision_codim_vids)
1460 {
1461 is_on_surface[vid] = true;
1462 }
1463
1464 Eigen::SparseMatrix<double> displacement_map;
1465 if (!displacement_map_entries.empty())
1466 {
1467 displacement_map.resize(collision_vertices.rows(), n_bases);
1468 displacement_map.setFromTriplets(displacement_map_entries.begin(), displacement_map_entries.end());
1469 }
1470
1471 collision_mesh = ipc::CollisionMesh(
1472 is_on_surface, collision_vertices, collision_edges, collision_triangles,
1473 displacement_map);
1474
1475 collision_mesh.can_collide = [&collision_mesh, num_fe_collision_vertices](size_t vi, size_t vj) {
1476 // obstacles do not collide with other obstacles
1477 return collision_mesh.to_full_vertex_id(vi) < num_fe_collision_vertices
1478 || collision_mesh.to_full_vertex_id(vj) < num_fe_collision_vertices;
1479 };
1480
1481 collision_mesh.init_area_jacobians();
1482 }
1483
1485 {
1486 if (!mesh)
1487 {
1488 logger().error("Load the mesh first!");
1489 return;
1490 }
1491 if (n_bases <= 0)
1492 {
1493 logger().error("Build the bases first!");
1494 return;
1495 }
1496 if (assembler->name() == "OperatorSplitting")
1497 {
1499 avg_mass = 1;
1500 return;
1501 }
1502
1503 if (!problem->is_time_dependent())
1504 {
1505 avg_mass = 1;
1507 return;
1508 }
1509
1510 mass.resize(0, 0);
1511
1512 igl::Timer timer;
1513 timer.start();
1514 logger().info("Assembling mass mat...");
1515
1516 if (mixed_assembler != nullptr)
1517 {
1518 StiffnessMatrix velocity_mass;
1519 mass_matrix_assembler->assemble(mesh->is_volume(), n_bases, bases, geom_bases(), mass_ass_vals_cache, 0, velocity_mass, true);
1520
1521 std::vector<Eigen::Triplet<double>> mass_blocks;
1522 mass_blocks.reserve(velocity_mass.nonZeros());
1523
1524 for (int k = 0; k < velocity_mass.outerSize(); ++k)
1525 {
1526 for (StiffnessMatrix::InnerIterator it(velocity_mass, k); it; ++it)
1527 {
1528 mass_blocks.emplace_back(it.row(), it.col(), it.value());
1529 }
1530 }
1531
1532 mass.resize(n_bases * assembler->size(), n_bases * assembler->size());
1533 mass.setFromTriplets(mass_blocks.begin(), mass_blocks.end());
1534 mass.makeCompressed();
1535 }
1536 else
1537 {
1538 mass_matrix_assembler->assemble(mesh->is_volume(), n_bases, bases, geom_bases(), mass_ass_vals_cache, 0, mass, true);
1539 }
1540
1541 assert(mass.size() > 0);
1542
1543 avg_mass = 0;
1544 for (int k = 0; k < mass.outerSize(); ++k)
1545 {
1546
1547 for (StiffnessMatrix::InnerIterator it(mass, k); it; ++it)
1548 {
1549 assert(it.col() == k);
1550 avg_mass += it.value();
1551 }
1552 }
1553
1554 avg_mass /= mass.rows();
1555 logger().info("average mass {}", avg_mass);
1556
1557 if (args["solver"]["advanced"]["lump_mass_matrix"])
1558 {
1560 }
1561
1562 timer.stop();
1563 timings.assembling_mass_mat_time = timer.getElapsedTime();
1564 logger().info(" took {}s", timings.assembling_mass_mat_time);
1565
1566 stats.nn_zero = mass.nonZeros();
1567 stats.num_dofs = mass.rows();
1568 stats.mat_size = (long long)mass.rows() * (long long)mass.cols();
1569 logger().info("sparsity: {}/{}", stats.nn_zero, stats.mat_size);
1570 }
1571
1572 std::shared_ptr<RhsAssembler> State::build_rhs_assembler(
1573 const int n_bases_,
1574 const std::vector<basis::ElementBases> &bases_,
1575 const assembler::AssemblyValsCache &ass_vals_cache_) const
1576 {
1577 json rhs_solver_params = args["solver"]["linear"];
1578 if (!rhs_solver_params.contains("Pardiso"))
1579 rhs_solver_params["Pardiso"] = {};
1580 rhs_solver_params["Pardiso"]["mtype"] = -2; // matrix type for Pardiso (2 = SPD)
1581
1582 const int size = problem->is_scalar() ? 1 : mesh->dimension();
1583
1584 return std::make_shared<RhsAssembler>(
1588 n_bases_, size, bases_, geom_bases(), ass_vals_cache_, *problem,
1589 args["space"]["advanced"]["bc_method"],
1590 rhs_solver_params);
1591 }
1592
1593 std::shared_ptr<PressureAssembler> State::build_pressure_assembler(
1594 const int n_bases_,
1595 const std::vector<basis::ElementBases> &bases_) const
1596 {
1597 const int size = problem->is_scalar() ? 1 : mesh->dimension();
1598
1599 return std::make_shared<PressureAssembler>(
1605 n_bases_, size, bases_, geom_bases(), *problem);
1606 }
1607
1609 {
1610 if (!mesh)
1611 {
1612 logger().error("Load the mesh first!");
1613 return;
1614 }
1615 if (n_bases <= 0)
1616 {
1617 logger().error("Build the bases first!");
1618 return;
1619 }
1620
1621 igl::Timer timer;
1622
1623 json p_params = {};
1624 p_params["formulation"] = assembler->name();
1625 p_params["root_path"] = root_path();
1626 {
1627 RowVectorNd min, max, delta;
1628 mesh->bounding_box(min, max);
1629 delta = (max - min) / 2. + min;
1630 if (mesh->is_volume())
1631 p_params["bbox_center"] = {delta(0), delta(1), delta(2)};
1632 else
1633 p_params["bbox_center"] = {delta(0), delta(1)};
1634 }
1635 problem->set_parameters(p_params);
1636
1637 rhs.resize(0, 0);
1638
1639 timer.start();
1640 logger().info("Assigning rhs...");
1641
1643 solve_data.rhs_assembler->assemble(mass_matrix_assembler->density(), rhs);
1644 rhs *= -1;
1645
1646 // if(problem->is_mixed())
1647 if (mixed_assembler != nullptr)
1648 {
1649 const int prev_size = rhs.size();
1650 const int n_larger = n_pressure_bases + (use_avg_pressure ? (assembler->is_fluid() ? 1 : 0) : 0);
1651 rhs.conservativeResize(prev_size + n_larger, rhs.cols());
1652 if (assembler->name() == "OperatorSplitting")
1653 {
1655 return;
1656 }
1657 // Divergence free rhs
1658 if (assembler->name() != "Bilaplacian" || local_neumann_boundary.empty())
1659 {
1660 rhs.block(prev_size, 0, n_larger, rhs.cols()).setZero();
1661 }
1662 else
1663 {
1664 Eigen::MatrixXd tmp(n_pressure_bases, 1);
1665 tmp.setZero();
1666
1667 std::shared_ptr<RhsAssembler> tmp_rhs_assembler = build_rhs_assembler(
1669
1670 tmp_rhs_assembler->set_bc(std::vector<LocalBoundary>(), std::vector<int>(), n_boundary_samples(), local_neumann_boundary, tmp);
1671 rhs.block(prev_size, 0, n_larger, rhs.cols()) = tmp;
1672 }
1673 }
1674
1675 timer.stop();
1676 timings.assigning_rhs_time = timer.getElapsedTime();
1677 logger().info(" took {}s", timings.assigning_rhs_time);
1678 }
1679
1680 void State::solve_problem(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
1681 {
1682 if (!mesh)
1683 {
1684 logger().error("Load the mesh first!");
1685 return;
1686 }
1687 if (n_bases <= 0)
1688 {
1689 logger().error("Build the bases first!");
1690 return;
1691 }
1692
1693 // if (rhs.size() <= 0)
1694 // {
1695 // logger().error("Assemble the rhs first!");
1696 // return;
1697 // }
1698
1699 // sol.resize(0, 0);
1700 // pressure.resize(0, 0);
1701 stats.spectrum.setZero();
1702
1703 igl::Timer timer;
1704 timer.start();
1705 logger().info("Solving {}", assembler->name());
1706
1707 init_solve(sol, pressure);
1708
1709 if (problem->is_time_dependent())
1710 {
1711 const double t0 = args["time"]["t0"];
1712 const int time_steps = args["time"]["time_steps"];
1713 const double dt = args["time"]["dt"];
1714
1715 // Pre log the output path for easier watching
1716 if (args["output"]["advanced"]["save_time_sequence"])
1717 {
1718 logger().info("Time sequence of simulation will be written to: \"{}\"",
1719 resolve_output_path(args["output"]["paraview"]["file_name"]));
1720 }
1721
1722 if (assembler->name() == "NavierStokes")
1723 solve_transient_navier_stokes(time_steps, t0, dt, sol, pressure);
1724 else if (assembler->name() == "OperatorSplitting")
1725 solve_transient_navier_stokes_split(time_steps, dt, sol, pressure);
1726 else if (is_homogenization())
1727 solve_homogenization(time_steps, t0, dt, sol);
1728 else if (is_problem_linear())
1729 solve_transient_linear(time_steps, t0, dt, sol, pressure);
1730 else if (!assembler->is_linear() && problem->is_scalar())
1731 throw std::runtime_error("Nonlinear scalar problems are not supported yet!");
1732 else
1733 solve_transient_tensor_nonlinear(time_steps, t0, dt, sol);
1734 }
1735 else
1736 {
1737 if (assembler->name() == "NavierStokes")
1738 solve_navier_stokes(sol, pressure);
1739 else if (is_homogenization())
1740 solve_homogenization(/* time steps */ 0, /* t0 */ 0, /* dt */ 0, sol);
1741 else if (is_problem_linear())
1742 {
1743 init_linear_solve(sol);
1744 solve_linear(sol, pressure);
1746 cache_transient_adjoint_quantities(0, sol, Eigen::MatrixXd::Zero(mesh->dimension(), mesh->dimension()));
1747 }
1748 else if (!assembler->is_linear() && problem->is_scalar())
1749 throw std::runtime_error("Nonlinear scalar problems are not supported yet!");
1750 else
1751 {
1755 cache_transient_adjoint_quantities(0, sol, Eigen::MatrixXd::Zero(mesh->dimension(), mesh->dimension()));
1756
1757 const std::string state_path = resolve_output_path(args["output"]["data"]["state"]);
1758 if (!state_path.empty())
1759 write_matrix(state_path, "u", sol);
1760 }
1761 }
1762
1763 timer.stop();
1764 timings.solving_time = timer.getElapsedTime();
1765 logger().info(" took {}s", timings.solving_time);
1766 }
1767
1768} // namespace polyfem
int V
ElementAssemblyValues vals
Definition Assembler.cpp:22
int x
DECLARE_DIFFSCALAR_BASE()
io::OutRuntimeData timings
runtime statistics
Definition State.hpp:590
assembler::MacroStrainValue macro_strain_constraint
Definition State.hpp:712
std::shared_ptr< utils::PeriodicBoundary > periodic_bc
periodic BC and periodic mesh utils
Definition State.hpp:386
int n_bases
number of bases
Definition State.hpp:178
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
Definition State.hpp:180
std::string formulation() const
return the formulation (checks if the problem is scalar or not and deals with multiphysics)
Definition State.cpp:328
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
Definition State.hpp:196
int n_boundary_samples() const
quadrature used for projecting boundary conditions
Definition State.hpp:263
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
double starting_max_edge_length
Definition State.hpp:594
std::vector< mesh::LocalBoundary > local_pressure_boundary
mapping from elements to nodes for pressure boundary conditions
Definition State.hpp:442
void sol_to_pressure(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
splits the solution in solution and pressure for mixed problems
Definition State.cpp:371
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...
Definition State.cpp:1484
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
Definition State.hpp:455
void build_polygonal_basis()
builds bases for polygons, called inside build_basis
Definition State.cpp:1065
mesh::Obstacle obstacle
Obstacles used in collisions.
Definition State.hpp:473
std::shared_ptr< assembler::Assembler > assembler
assemblers
Definition State.hpp:155
ipc::CollisionMesh periodic_collision_mesh
IPC collision mesh under periodic BC.
Definition State.hpp:523
void build_periodic_collision_mesh()
Definition State.cpp:1228
io::OutGeometryData out_geom
visualization stuff
Definition State.hpp:588
bool use_avg_pressure
use average pressure for stokes problem to fix the additional dofs, true by default if false,...
Definition State.hpp:211
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
Definition State.hpp:451
ipc::CollisionMesh collision_mesh
IPC collision mesh.
Definition State.hpp:520
std::vector< basis::ElementBases > geom_bases_
Geometric mapping bases, if the elements are isoparametric, this list is empty.
Definition State.hpp:175
solver::CacheLevel optimization_enabled
Definition State.hpp:654
Eigen::VectorXi in_primitive_to_primitive
maps in vertices/edges/faces/cells to polyfem vertices/edges/faces/cells
Definition State.hpp:457
std::map< int, basis::InterfaceData > poly_edge_to_data
nodes on the boundary of polygonal elements, used for harmonic bases
Definition State.hpp:446
std::vector< int > pressure_boundary_nodes
list of neumann boundary nodes
Definition State.hpp:434
std::vector< basis::ElementBases > pressure_bases
FE pressure bases for mixed elements, the size is #elements.
Definition State.hpp:173
bool is_homogenization() const
Definition State.hpp:718
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
Definition State.hpp:593
std::shared_ptr< assembler::Mass > mass_matrix_assembler
Definition State.hpp:157
void build_node_mapping()
build the mapping from input nodes to polyfem nodes
Definition State.cpp:245
std::vector< int > primitive_to_node() const
Definition State.cpp:228
bool has_dhat
stores if input json contains dhat
Definition State.hpp:579
std::shared_ptr< polyfem::mesh::MeshNodes > pressure_mesh_nodes
Definition State.hpp:193
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
std::shared_ptr< polyfem::mesh::MeshNodes > mesh_nodes
Mapping from input nodes to FE nodes.
Definition State.hpp:193
StiffnessMatrix mass
Mass matrix, it is computed only for time dependent problems.
Definition State.hpp:202
void solve_homogenization(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol)
std::vector< int > dirichlet_nodes
per node dirichlet
Definition State.hpp:448
bool has_periodic_bc() const
Definition State.hpp:387
void solve_problem(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves the problems
Definition State.cpp:1680
void build_basis()
builds the bases step 2 of solve modifies bases, pressure_bases, geom_bases_, boundary_nodes,...
Definition State.cpp:506
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
Definition State.hpp:168
std::vector< RowVectorNd > dirichlet_nodes_position
Definition State.hpp:449
std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > polys_3d
polyhedra, used since poly have no geom mapping
Definition State.hpp:187
std::vector< int > node_to_primitive() const
Definition State.cpp:235
json args
main input arguments containing all defaults
Definition State.hpp:101
io::OutStatsData stats
Other statistics.
Definition State.hpp:592
std::vector< RowVectorNd > neumann_nodes_position
Definition State.hpp:452
void assemble_rhs()
compute rhs, step 3 of solve build rhs vector based on defined basis and given rhs of the problem mod...
Definition State.cpp:1608
double min_boundary_edge_length
Definition State.hpp:595
Eigen::VectorXi periodic_collision_mesh_to_basis
index mapping from periodic 2x2 collision mesh to FE periodic mesh
Definition State.hpp:525
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:171
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
Definition State.hpp:182
assembler::AssemblyValsCache pressure_ass_vals_cache
used to store assembly values for pressure for small problems
Definition State.hpp:199
std::vector< mesh::LocalBoundary > local_boundary
mapping from elements to nodes for dirichlet boundary conditions
Definition State.hpp:438
double avg_mass
average system mass, used for contact with IPC
Definition State.hpp:204
bool iso_parametric() const
check if using iso parametric bases
Definition State.cpp:462
std::shared_ptr< assembler::RhsAssembler > build_rhs_assembler() const
build a RhsAssembler for the problem
Definition State.hpp:248
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
Definition State.hpp:436
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
Definition State.hpp:197
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
Definition State.hpp:432
solver::SolveData solve_data
timedependent stuff cached
Definition State.hpp:321
std::shared_ptr< polyfem::mesh::MeshNodes > geom_mesh_nodes
Definition State.hpp:193
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 have contact
Definition State.hpp:556
void build_collision_mesh()
extracts the boundary mesh for collision, called in build_basis
Definition State.cpp:1361
Eigen::VectorXi disc_orders
vector of discretization orders, used when not all elements have the same degree, one per element
Definition State.hpp:190
std::vector< mesh::LocalBoundary > local_neumann_boundary
mapping from elements to nodes for neumann boundary conditions
Definition State.hpp:440
std::shared_ptr< assembler::PressureAssembler > build_pressure_assembler() const
Definition State.hpp:256
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.
Definition State.hpp:407
std::shared_ptr< assembler::MixedAssembler > mixed_assembler
Definition State.hpp:159
std::map< int, Eigen::MatrixXd > polys
polygons, used since poly have no geom mapping
Definition State.hpp:185
Eigen::MatrixXd rhs
System right-hand side.
Definition State.hpp:207
std::unordered_map< int, std::vector< mesh::LocalBoundary > > local_pressure_cavity
mapping from elements to nodes for pressure boundary conditions
Definition State.hpp:444
StiffnessMatrix basis_nodes_to_gbasis_nodes
Definition State.hpp:706
std::string velocity() const
Definition Units.hpp:29
static double convert(const json &val, const std::string &unit_type)
Definition Units.cpp:31
const std::string & length() const
Definition Units.hpp:19
static bool is_elastic_material(const std::string &material)
utility to check if material is one of the elastic materials
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
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,...
assemble matrix based on the local assembler local assembler is eg Laplace, LinearElasticity etc
void init(const int dim, const json &param)
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 bernstein, const bool serendipity, const bool has_polys, const bool is_geom_bases, const bool use_corner_quadrature, 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 bernstein, const bool serendipity, const bool has_polys, const bool is_geom_bases, const bool use_corner_quadrature, 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
Definition OutData.cpp:2352
static void extract_boundary_mesh(const mesh::Mesh &mesh, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< mesh::LocalBoundary > &total_local_boundary, Eigen::MatrixXd &node_positions, Eigen::MatrixXi &boundary_edges, Eigen::MatrixXi &boundary_triangles, std::vector< Eigen::Triplet< double > > &displacement_map_entries)
extracts the boundary mesh
Definition OutData.cpp:145
double assembling_stiffness_mat_time
time to assembly
Definition OutData.hpp:350
double assigning_rhs_time
time to computing the rhs
Definition OutData.hpp:354
double assembling_mass_mat_time
time to assembly mass
Definition OutData.hpp:352
double building_basis_time
time to construct the basis
Definition OutData.hpp:344
double solving_time
time to solve
Definition OutData.hpp:356
double computing_poly_basis_time
time to build the polygonal/polyhedral bases
Definition OutData.hpp:348
int n_flipped
number of flipped elements, compute only when using count_flipped_els (false by default)
Definition OutData.hpp:397
Eigen::Vector4d spectrum
spectrum of the stiffness matrix, enable only if POLYSOLVE_WITH_SPECTRA is ON (off by default)
Definition OutData.hpp:371
void count_flipped_elements(const polyfem::mesh::Mesh &mesh, const std::vector< polyfem::basis::ElementBases > &gbases)
counts the number of flipped elements
Definition OutData.cpp:2537
void compute_mesh_size(const polyfem::mesh::Mesh &mesh_in, const std::vector< polyfem::basis::ElementBases > &bases_in, const int n_samples, const bool use_curved_mesh_size)
computes the mesh size, it samples every edges n_samples times uses curved_mesh_size (false by defaul...
Definition OutData.cpp:2450
long long nn_zero
non zeros and sytem matrix size num dof is the total dof in the system
Definition OutData.hpp:389
double mesh_size
max edge lenght
Definition OutData.hpp:378
void reset()
clears all stats
Definition OutData.cpp:2528
double min_edge_length
min edge lenght
Definition OutData.hpp:380
bool is_volume() const override
checks if mesh is volume
Definition Mesh3D.hpp:28
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:39
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
Definition Mesh.hpp:161
const Eigen::MatrixXi & e() const
Definition Obstacle.hpp:44
const Eigen::MatrixXi & f() const
Definition Obstacle.hpp:43
const Eigen::MatrixXd & v() const
Definition Obstacle.hpp:41
const Eigen::VectorXi & codim_v() const
Definition Obstacle.hpp:42
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
Definition APriori.cpp:242
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.
Definition MatrixIO.cpp:18
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.
Definition MatrixIO.cpp:42
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.
Definition JSONUtils.hpp:38
void append_rows(DstMat &dst, const SrcMat &src)
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
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)
Definition State.cpp:389
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 > VectorNd
Definition Types.hpp:11
nlohmann::json json
Definition Common.hpp:9
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22