PolyFEM
Loading...
Searching...
No Matches
GeometryReader.cpp
Go to the documentation of this file.
1#include "GeometryReader.hpp"
2
3#include <polyfem/mesh/Mesh.hpp>
4#include <polyfem/mesh/MeshUtils.hpp>
5#include <polyfem/io/MshReader.hpp>
6#include <polyfem/utils/StringUtils.hpp>
7
8#include <polyfem/utils/JSONUtils.hpp>
9#include <polyfem/utils/Selection.hpp>
10#include <polyfem/utils/Logger.hpp>
11
12#include <Eigen/Core>
13
14#include <igl/edges.h>
15#include <igl/boundary_facets.h>
16
17#include <strnatcmp.h>
18#include <glob/glob.h>
19#include <filesystem>
20
21namespace polyfem::mesh
22{
23 using namespace polyfem::utils;
24
25 void apply_geometry_selection(
26 Mesh &mesh,
27 const json &geometry_selection,
28 const std::string &root_path)
29 {
30 if (geometry_selection.is_object() && geometry_selection.contains("same_as_volume"))
31 {
32 if (!geometry_selection["same_as_volume"].get<bool>())
33 return;
34
35 std::vector<int> geometry_ids(mesh.n_elements());
36 for (int e = 0; e < mesh.n_elements(); ++e)
37 geometry_ids[e] = mesh.get_body_id(e);
38 mesh.set_geometry_ids(geometry_ids);
39 return;
40 }
41
42 Selection::BBox bbox;
43 mesh.bounding_box(bbox[0], bbox[1]);
44 const auto selections = Selection::build_selections(geometry_selection, bbox, root_path);
45
46 Eigen::MatrixXd barycenters;
47 mesh.compute_element_barycenters(barycenters);
48 std::vector<int> geometry_ids(mesh.n_elements(), 0);
49 for (int e = 0; e < mesh.n_elements(); ++e)
50 {
51 const auto vertices = mesh.element_vertices(e);
52 for (const auto &selection : selections)
53 {
54 if (selection->inside(e, vertices, barycenters.row(e)))
55 {
56 geometry_ids[e] = selection->id(e, vertices, barycenters.row(e));
57 break;
58 }
59 }
60 }
61 mesh.set_geometry_ids(geometry_ids);
62 }
63
64 std::unique_ptr<Mesh> read_fem_mesh(
65 const Units &units,
66 const json &j_mesh,
67 const std::string &root_path,
68 const bool non_conforming)
69 {
70 if (!is_param_valid(j_mesh, "mesh"))
71 log_and_throw_error("Mesh {} is mising a \"mesh\" field!", j_mesh);
72
73 if (j_mesh["extract"].get<std::string>() != "volume")
74 log_and_throw_error("Only volumetric elements are implemented for FEM meshes!");
75
76 std::unique_ptr<Mesh> mesh = Mesh::create(resolve_path(j_mesh["mesh"], root_path), non_conforming);
77
78 // --------------------------------------------------------------------
79
80 // NOTE: Normaliziation is done before transformations are applied and/or any selection operators
81 if (j_mesh["advanced"]["normalize_mesh"])
82 mesh->normalize();
83
84 // --------------------------------------------------------------------
85
86 Selection::BBox bbox;
87 mesh->bounding_box(bbox[0], bbox[1]);
88
89 const std::string unit = j_mesh["unit"];
90 double unit_scale = 1;
91 if (!unit.empty())
92 unit_scale = Units::convert(1, unit, units.length());
93
94 {
95 MatrixNd A;
96 VectorNd b;
97 construct_affine_transformation(
98 unit_scale,
99 j_mesh["transformation"],
100 (bbox[1] - bbox[0]).cwiseAbs().transpose(),
101 A, b);
102 mesh->apply_affine_transformation(A, b);
103 }
104
105 mesh->bounding_box(bbox[0], bbox[1]);
106
107 // --------------------------------------------------------------------
108 std::vector<std::shared_ptr<Selection>> surface_selections =
109 is_param_valid(j_mesh, "surface_selection") ? Selection::build_selections(j_mesh["surface_selection"], bbox, root_path) : std::vector<std::shared_ptr<Selection>>();
110
111 // --------------------------------------------------------------------
112
113 const int n_refs = j_mesh["n_refs"];
114 const double refinement_location = j_mesh["advanced"]["refinement_location"];
115 // TODO: renable this
116 // if (n_refs <= 0 && args["poly_bases"] == "MFSHarmonic" && mesh->has_poly())
117 // {
118 // if (args["force_no_ref_for_harmonic"])
119 // logger().warn("Using harmonic bases without refinement");
120 // else
121 // n_refs = 1;
122 // }
123 if (n_refs > 0)
124 {
125 if (mesh->has_boundary_ids() && surface_selections.empty())
126 log_and_throw_error("Unable to refine a mesh with stored surface selections; provide an explicit surface_selection to recompute them after refinement.");
127
128 // Check if the stored volume selection is uniform.
129 assert(mesh->n_elements() > 0);
130 const int uniform_value = mesh->get_body_id(0);
131 for (int i = 1; i < mesh->n_elements(); ++i)
132 if (mesh->get_body_id(i) != uniform_value)
133 log_and_throw_error("Unable to apply stored nonuniform volume_selection because n_refs={} > 0!", n_refs);
134
135 logger().info("Performing global h-refinement with {} refinements", n_refs);
136 mesh->refine(n_refs, refinement_location);
137 mesh->set_body_ids(std::vector<int>(mesh->n_elements(), uniform_value));
138 }
139
140 // --------------------------------------------------------------------
141
142 if (j_mesh["advanced"]["min_component"].get<int>() != -1)
143 log_and_throw_error("Option \"min_component\" in geometry not implement yet!");
144 // TODO:
145 // if (args["min_component"] > 0) {
146 // Eigen::SparseMatrix<int> adj;
147 // igl::facet_adjacency_matrix(boundary_triangles, adj);
148 // Eigen::MatrixXi C, counts;
149 // igl::connected_components(adj, C, counts);
150 // std::vector<int> valid;
151 // const int min_count = args["min_component"];
152 // for (int i = 0; i < counts.size(); ++i) {
153 // if (counts(i) >= min_count) {
154 // valid.push_back(i);
155 // }
156 // }
157 // tris.clear();
158 // for (int i = 0; i < C.size(); ++i) {
159 // for (int v : valid) {
160 // if (v == C(i)) {
161 // tris.emplace_back(boundary_triangles(i, 0), boundary_triangles(i, 1), boundary_triangles(i, 2));
162 // break;
163 // }
164 // }
165 // }
166 // boundary_triangles.resize(tris.size(), 3);
167 // for (int i = 0; i < tris.size(); ++i) {
168 // boundary_triangles.row(i) << std::get<0>(tris[i]), std::get<1>(tris[i]), std::get<2>(tris[i]);
169 // }
170 // }
171
172 // --------------------------------------------------------------------
173
174 if (j_mesh["advanced"]["force_linear_geometry"].get<bool>())
175 log_and_throw_error("Option \"force_linear_geometry\" in geometry not implement yet!");
176 // TODO:
177 // if (!iso_parametric()) {
178 // if (args["force_linear_geometry"] || mesh->orders().size() <= 0) {
179 // geom_disc_orders.resizeLike(disc_orders);
180 // geom_disc_orders.setConstant(1);
181 // } else {
182 // geom_disc_orders = mesh->orders();
183 // }
184 // }
185
186 // --------------------------------------------------------------------
187
188 const std::vector<std::shared_ptr<Selection>> node_selections =
189 is_param_valid(j_mesh, "point_selection") ? Selection::build_selections(j_mesh["point_selection"], bbox, root_path) : std::vector<std::shared_ptr<Selection>>();
190
191 if (!node_selections.empty())
192 {
193 bool boundary_only = true;
194 for (const auto &selection : node_selections)
195 {
196 if (!selection->boundary_only())
197 {
198 boundary_only = false;
199 break;
200 }
201 }
202
203 mesh->compute_node_ids([&](const size_t n_id, const RowVectorNd &p, bool is_boundary) {
204 if (boundary_only && !is_boundary)
205 return -1;
206
207 const std::vector<int> tmp = {int(n_id)};
208 for (const auto &selection : node_selections)
209 {
210 if (selection->boundary_only() && !is_boundary)
211 continue;
212
213 if (selection->inside(n_id, tmp, p))
214 return selection->id(n_id, tmp, p);
215 }
216 return std::numeric_limits<int>::max(); // default for no selected boundary
217 });
218 }
219
220 if (!j_mesh["curve_selection"].is_null())
221 log_and_throw_error("Geometry curve selections are not implemented!");
222
223 // --------------------------------------------------------------------
224
225 if (!surface_selections.empty())
226 {
227 bool boundary_only = true;
228 for (const auto &selection : surface_selections)
229 {
230 if (!selection->boundary_only())
231 {
232 boundary_only = false;
233 break;
234 }
235 }
236
237 mesh->compute_boundary_ids([&](const size_t p_id, const std::vector<int> &vs, const RowVectorNd &p, bool is_boundary) {
238 if (boundary_only && !is_boundary)
239 return -1;
240
241 for (const auto &selection : surface_selections)
242 {
243 if (selection->boundary_only() && !is_boundary)
244 continue;
245
246 if (selection->inside(p_id, vs, p))
247 return selection->id(p_id, vs, p);
248 }
249 return std::numeric_limits<int>::max(); // default for no selected boundary
250 });
251 }
252
253 // --------------------------------------------------------------------
254
255 // If the selection is of the form {"id_offset": ...}
256 const json volume_selection = j_mesh["volume_selection"];
257 if (volume_selection.is_object()
258 && volume_selection.size() == 1
259 && volume_selection.contains("id_offset"))
260 {
261 const int id_offset = volume_selection["id_offset"].get<int>();
262 if (id_offset != 0)
263 {
264 const int n_body_ids = mesh->n_elements();
265 std::vector<int> body_ids(n_body_ids);
266 for (int i = 0; i < n_body_ids; ++i)
267 body_ids[i] = mesh->get_body_id(i) + id_offset;
268 mesh->set_body_ids(body_ids);
269 }
270 }
271 else
272 {
273 // Specified volume selection has priority over mesh's stored ids
274 std::vector<std::shared_ptr<Selection>> volume_selections =
275 Selection::build_selections(volume_selection, bbox, root_path);
276
277 // Append the mesh's stored ids to the volume selection as a lowest priority selection
278 if (mesh->has_body_ids())
279 volume_selections.push_back(std::make_shared<SpecifiedSelection>(mesh->get_body_ids()));
280
281 mesh->compute_body_ids([&](const size_t cell_id, const std::vector<int> &vs, const RowVectorNd &p) -> int {
282 for (const auto &selection : volume_selections)
283 {
284 // TODO: add vs to compute_body_ids
285 if (selection->inside(cell_id, vs, p))
286 return selection->id(cell_id, vs, p);
287 }
288 return 0;
289 });
290 }
291
292 // --------------------------------------------------------------------
293
294 if (is_param_valid(j_mesh, "geometry_selection"))
295 apply_geometry_selection(*mesh, j_mesh["geometry_selection"], root_path);
296
297 // --------------------------------------------------------------------
298
299 return mesh;
300 }
301
302 // ========================================================================
303
304 std::unique_ptr<Mesh> read_fem_geometry(
305 const Units &units,
306 const json &geometry,
307 const std::string &root_path,
308 const std::vector<std::string> &_names,
309 const std::vector<Eigen::MatrixXd> &_vertices,
310 const std::vector<Eigen::MatrixXi> &_cells,
311 const bool non_conforming)
312 {
313 // TODO: fix me for hdf5
314 // {
315 // int index = -1;
316 // for (int i = 0; i < names.size(); ++i)
317 // {
318 // if (names[i] == args["meshes"])
319 // {
320 // index = i;
321 // break;
322 // }
323 // }
324 // assert(index >= 0);
325 // if (vertices[index].cols() == 2)
326 // mesh = std::make_unique<polyfem::CMesh2D>();
327 // else
328 // mesh = std::make_unique<polyfem::Mesh3D>();
329 // mesh->build_from_matrices(vertices[index], cells[index]);
330 // }
331 assert(_names.empty());
332 assert(_vertices.empty());
333 assert(_cells.empty());
334
335 // --------------------------------------------------------------------
336
337 if (geometry.empty())
338 log_and_throw_error("Provided geometry is empty!");
339
340 std::vector<json> geometries = utils::json_as_array(geometry);
341
342 // --------------------------------------------------------------------
343
344 std::unique_ptr<Mesh> mesh = nullptr;
345
346 for (const json &geometry : geometries)
347 {
348 if (!geometry["enabled"].get<bool>() || geometry["is_obstacle"].get<bool>())
349 continue;
350
351 if (geometry["type"] != "mesh" && geometry["type"] != "mesh_array")
352 log_and_throw_error("Invalid geometry type \"{}\" for FEM mesh!", geometry["type"]);
353
354 const std::unique_ptr<Mesh> tmp_mesh = read_fem_mesh(units, geometry, root_path, non_conforming);
355
356 if (mesh == nullptr)
357 mesh = tmp_mesh->copy();
358 else
359 mesh->append(tmp_mesh);
360
361 if (geometry["type"] == "mesh_array")
362 {
363 Selection::BBox bbox;
364 tmp_mesh->bounding_box(bbox[0], bbox[1]);
365
366 const long dim = tmp_mesh->dimension();
367 const bool is_offset_relative = geometry["array"]["relative"];
368 const double offset = geometry["array"]["offset"];
369 const VectorNd dimensions = (bbox[1] - bbox[0]);
370 const VectorNi size = geometry["array"]["size"];
371
372 for (int i = 0; i < size[0]; ++i)
373 {
374 for (int j = 0; j < size[1]; ++j)
375 {
376 for (int k = 0; k < (size.size() > 2 ? size[2] : 1); ++k)
377 {
378 if (i == 0 && j == 0 && k == 0)
379 continue;
380
381 RowVectorNd translation = offset * Eigen::RowVector3d(i, j, k).head(dim);
382 if (is_offset_relative)
383 translation.array() *= dimensions.array();
384
385 const std::unique_ptr<Mesh> copy_mesh = tmp_mesh->copy();
386 copy_mesh->apply_affine_transformation(MatrixNd::Identity(dim, dim), translation);
387 mesh->append(copy_mesh);
388 }
389 }
390 }
391 }
392 }
393
394 // --------------------------------------------------------------------
395
396 return mesh;
397 }
398
399 // ========================================================================
400
401 void read_obstacle_mesh(
402 const Units &units,
403 const json &j_mesh,
404 const std::string &root_path,
405 const int dim,
406 Eigen::MatrixXd &vertices,
407 Eigen::VectorXi &codim_vertices,
408 Eigen::MatrixXi &codim_edges,
409 Eigen::MatrixXi &faces)
410 {
411 if (!is_param_valid(j_mesh, "mesh"))
412 log_and_throw_error("Mesh obstacle {} is mising a \"mesh\" field!", j_mesh);
413
414 const std::string mesh_path = resolve_path(j_mesh["mesh"], root_path);
415
416 bool read_success = read_surface_mesh(
417 mesh_path, vertices, codim_vertices, codim_edges, faces);
418
419 if (!read_success)
420 // error already logged in read_surface_mesh()
421 throw std::runtime_error(fmt::format("Unable to read mesh: {}", mesh_path));
422
423 const int prev_dim = vertices.cols();
424 vertices.conservativeResize(vertices.rows(), dim);
425 if (prev_dim < dim)
426 vertices.rightCols(dim - prev_dim).setZero();
427
428 // --------------------------------------------------------------------
429
430 {
431 const std::string unit = j_mesh["unit"];
432 double unit_scale = 1;
433 if (!unit.empty())
434 unit_scale = Units::convert(1, unit, units.length());
435
436 const VectorNd mesh_dimensions = (vertices.colwise().maxCoeff() - vertices.colwise().minCoeff()).cwiseAbs();
437 MatrixNd A;
438 VectorNd b;
439 construct_affine_transformation(unit_scale, j_mesh["transformation"], mesh_dimensions, A, b);
440 vertices = vertices * A.transpose();
441 vertices.rowwise() += b.transpose();
442 }
443
444 std::string extract = j_mesh["extract"];
445 // Default: "volume" clashes with defaults for non obstacle, here assume volume is suface
446 if (extract == "volume")
447 extract = "surface";
448
449 if (extract == "points")
450 {
451 // points -> vertices (drop edges and faces)
452 codim_edges.resize(0, 0);
453 faces.resize(0, 0);
454 codim_vertices.resize(vertices.rows());
455 for (int i = 0; i < codim_vertices.size(); ++i)
456 codim_vertices[i] = i;
457 }
458 else if (extract == "edges" && faces.size() != 0)
459 {
460 // edges -> edges (drop faces)
461 Eigen::MatrixXi edges;
462 igl::edges(faces, edges);
463 faces.resize(0, 0);
464 codim_edges.conservativeResize(codim_edges.rows() + edges.rows(), 2);
465 codim_edges.bottomRows(edges.rows()) = edges;
466 }
467 else if (extract == "surface" && dim == 2 && faces.size() != 0)
468 {
469 // surface (2D) -> boundary edges (drop faces and interior edges)
470 Eigen::MatrixXi boundary_edges;
471 igl::boundary_facets(faces, boundary_edges);
472 codim_edges.conservativeResize(codim_edges.rows() + boundary_edges.rows(), 2);
473 codim_edges.bottomRows(boundary_edges.rows()) = boundary_edges;
474 faces.resize(0, 0); // Clear faces
475 }
476 // surface (3D) -> boundary faces
477 // No need to do anything for (extract == "surface" && dim == 3) since we used read_surface_mesh
478 else if (extract == "volume")
479 {
480 // volume -> undefined
481 log_and_throw_error("Volumetric elements not supported for collision obstacles!");
482 }
483
484 if (j_mesh["n_refs"].get<int>() != 0)
485 {
486 if (faces.size() != 0)
487 log_and_throw_error("Option \"n_refs\" for triangle obstacles not implement yet!");
488
489 const int n_refs = j_mesh["n_refs"];
490 const double refinement_location = j_mesh["advanced"]["refinement_location"];
491 for (int i = 0; i < n_refs; i++)
492 {
493 const size_t n_vertices = vertices.rows();
494 const size_t n_edges = codim_edges.rows();
495 vertices.conservativeResize(n_vertices + n_edges, vertices.cols());
496 codim_edges.conservativeResize(2 * n_edges, codim_edges.cols());
497 for (size_t ei = 0; ei < n_edges; ei++)
498 {
499 const int v0i = codim_edges(ei, 0);
500 const int v1i = codim_edges(ei, 1);
501 const int v2i = n_vertices + ei;
502 vertices.row(v2i) = (vertices.row(v1i) - vertices.row(v0i)) * refinement_location + vertices.row(v0i);
503 codim_edges.row(ei) << v0i, v2i;
504 codim_edges.row(n_edges + ei) << v2i, v1i;
505 }
506 }
507 }
508 }
509
510 // ========================================================================
511
512 Obstacle read_obstacle_geometry(
513 const Units &units,
514 const json &geometry,
515 const std::vector<json> &displacements,
516 const std::vector<json> &dirichlets,
517 const std::string &root_path,
518 const int dim,
519 const std::vector<std::string> &_names,
520 const std::vector<Eigen::MatrixXd> &_vertices,
521 const std::vector<Eigen::MatrixXi> &_cells,
522 const bool non_conforming)
523 {
524 // TODO: fix me for hdf5
525 // {
526 // int index = -1;
527 // for (int i = 0; i < names.size(); ++i)
528 // {
529 // if (names[i] == args["meshes"])
530 // {
531 // index = i;
532 // break;
533 // }
534 // }
535 // assert(index >= 0);
536 // if (vertices[index].cols() == 2)
537 // mesh = std::make_unique<polyfem::CMesh2D>();
538 // else
539 // mesh = std::make_unique<polyfem::Mesh3D>();
540 // mesh->build_from_matrices(vertices[index], cells[index]);
541 // }
542 assert(_names.empty());
543 assert(_vertices.empty());
544 assert(_cells.empty());
545
546 Obstacle obstacle;
547
548 if (geometry.empty())
549 return obstacle;
550
551 std::vector<json> geometries = utils::json_as_array(geometry);
552
553 for (const json &geometry : geometries)
554 {
555
556 if (!geometry["is_obstacle"].get<bool>())
557 continue;
558
559 if (!geometry["enabled"].get<bool>())
560 continue;
561
562 if (geometry["type"] == "mesh" || geometry["type"] == "mesh_array")
563 {
564 Eigen::MatrixXd vertices;
565 Eigen::VectorXi codim_vertices;
566 Eigen::MatrixXi codim_edges;
567 Eigen::MatrixXi faces;
568 read_obstacle_mesh(units,
569 geometry, root_path, dim, vertices, codim_vertices,
570 codim_edges, faces);
571
572 if (geometry["type"] == "mesh_array")
573 {
574 const Selection::BBox bbox{{vertices.colwise().minCoeff(), vertices.colwise().maxCoeff()}};
575
576 const bool is_offset_relative = geometry["array"]["relative"];
577 const double offset = geometry["array"]["offset"];
578 const VectorNd dimensions = (bbox[1] - bbox[0]);
579 const VectorNi size = geometry["array"]["size"];
580
581 const int N = size.head(dim).prod();
582 const int nV = vertices.rows(), nCV = codim_vertices.rows(), nCE = codim_edges.rows(), nF = faces.rows();
583
584 vertices.conservativeResize(N * nV, Eigen::NoChange);
585 codim_vertices.conservativeResize(N * nCV, Eigen::NoChange);
586 codim_edges.conservativeResize(N * nCE, Eigen::NoChange);
587 faces.conservativeResize(N * nF, Eigen::NoChange);
588
589 for (int i = 0; i < size[0]; ++i)
590 {
591 for (int j = 0; j < size[1]; ++j)
592 {
593 for (int k = 0; k < (size.size() > 2 ? size[2] : 1); ++k)
594 {
595 RowVectorNd translation = offset * Eigen::RowVector3d(i, j, k).head(vertices.cols());
596 if (is_offset_relative)
597 translation.array() *= dimensions.array();
598
599 int n = i * size[1] + j;
600 if (size.size() > 2)
601 n = n * size[2] + k;
602 if (n == 0)
603 continue;
604
605 vertices.middleRows(n * nV, nV) = vertices.topRows(nV).rowwise() + translation;
606 if (nCV)
607 codim_vertices.segment(n * nV, nV) = codim_vertices.head(nV).array() + n * nV;
608 if (nCE)
609 codim_edges.middleRows(n * nCE, nCE) = codim_edges.topRows(nCE).array() + n * nV;
610 if (nF)
611 faces.middleRows(n * nF, nF) = faces.topRows(nF).array() + n * nV;
612 }
613 }
614 }
615 }
616
617 json displacement = "{\"value\":[0, 0, 0]}"_json;
618 if (is_param_valid(geometry, "surface_selection"))
619 {
620 if (!geometry["surface_selection"].is_number())
621 log_and_throw_error("Invalid surface_selection for obstacle, needs to be an integer!");
622
623 const int id = geometry["surface_selection"];
624 for (const json &disp : dirichlets)
625 {
626 if ((disp["id"].is_string() && disp["id"].get<std::string>() == "all")
627 || (disp["id"].is_number_integer() && disp["id"].get<int>() == id))
628 {
629 displacement = disp;
630 break;
631 }
632 else if (disp["id"].is_array())
633 {
634 for (const json &disp_id : disp["id"])
635 {
636 assert(disp_id.is_number_integer());
637 if (disp_id.get<int>() == id)
638 {
639 displacement = disp;
640 break;
641 }
642 }
643 }
644 }
645 for (const json &disp : displacements)
646 {
647 if ((disp["id"].is_string() && disp["id"].get<std::string>() == "all")
648 || (disp["id"].is_number_integer() && disp["id"].get<int>() == id))
649 {
650 displacement = disp;
651 break;
652 }
653 else if (disp["id"].is_array())
654 {
655 for (const json &disp_id : disp["id"])
656 {
657 assert(disp_id.is_number_integer());
658 if (disp_id.get<int>() == id)
659 {
660 displacement = disp;
661 break;
662 }
663 }
664 }
665 }
666 }
667
668 obstacle.append_mesh(
669 vertices, codim_vertices, codim_edges, faces, displacement, root_path);
670 }
671 else if (geometry["type"] == "plane")
672 {
673 obstacle.append_plane(geometry["point"], geometry["normal"]);
674 }
675 else if (geometry["type"] == "ground")
676 {
677 VectorNd gravity = VectorNd::Zero(dim); // TODO: Expose as parameter
678 gravity[1] = -9.81;
679 const double height = geometry["height"];
680 assert(gravity.norm() != 0);
681 const VectorNd normal = -gravity.normalized();
682 const VectorNd point = height * normal; // origin + height * normal
683 obstacle.append_plane(point, normal);
684 }
685 else if (geometry["type"] == "mesh_sequence")
686 {
687 namespace fs = std::filesystem;
688 std::vector<fs::path> mesh_files;
689 if (geometry["mesh_sequence"].is_array())
690 {
691 mesh_files = geometry["mesh_sequence"].get<std::vector<fs::path>>();
692 }
693 else
694 {
695 assert(geometry["mesh_sequence"].is_string());
696 const fs::path meshes(resolve_path(geometry["mesh_sequence"], root_path));
697
698 if (fs::is_directory(meshes))
699 {
700 for (const auto &entry : std::filesystem::directory_iterator(meshes))
701 {
702 if (entry.is_regular_file())
703 mesh_files.push_back(entry.path());
704 }
705 }
706 else
707 {
708 mesh_files = glob::rglob(meshes.string());
709 }
710 // Sort the file names naturally
711 std::sort(mesh_files.begin(), mesh_files.end(), [](const fs::path &p1, const fs::path &p2) {
712 return strnatcmp(p1.string().c_str(), p2.string().c_str()) < 0;
713 });
714 }
715
716 std::vector<Eigen::MatrixXd> vertices(mesh_files.size());
717 Eigen::VectorXi codim_vertices;
718 Eigen::MatrixXi codim_edges;
719 Eigen::MatrixXi faces;
720
721 for (int i = 0; i < mesh_files.size(); ++i)
722 {
723 json jmesh = geometry;
724 jmesh["mesh"] = mesh_files[i];
725 jmesh["n_refs"] = 0;
726
727 Eigen::VectorXi tmp_codim_vertices;
728 Eigen::MatrixXi tmp_codim_edges;
729 Eigen::MatrixXi tmp_faces;
730 read_obstacle_mesh(units,
731 jmesh, root_path, dim, vertices[i],
732 tmp_codim_vertices, tmp_codim_edges, tmp_faces);
733 if (i == 0)
734 {
735 codim_vertices = tmp_codim_vertices;
736 codim_edges = tmp_codim_edges;
737 faces = tmp_faces;
738 }
739 else
740 {
741 assert((codim_vertices.array() == tmp_codim_vertices.array()).all());
742 assert((codim_edges.array() == tmp_codim_edges.array()).all());
743 assert((faces.array() == tmp_faces.array()).all());
744 }
745 }
746
747 obstacle.append_mesh_sequence(
748 vertices, codim_vertices, codim_edges, faces, geometry["fps"]);
749 }
750 else
751 {
752 log_and_throw_error("Invalid geometry type \"{}\" for obstacle!", geometry["type"]);
753 }
754 }
755
756 obstacle.set_units(units);
757 return obstacle;
758 }
759
760 // ========================================================================
761
762 void construct_affine_transformation(
763 const double unit_scale,
764 const json &transform,
765 const VectorNd &mesh_dimensions,
766 MatrixNd &A,
767 VectorNd &b)
768 {
769 const int dim = mesh_dimensions.size();
770
771 // -----
772 // Scale
773 // -----
774
775 RowVectorNd scale;
776 if (transform["dimensions"].is_array()) // default is nullptr
777 {
778 VectorNd modified_dimensions =
779 (mesh_dimensions.array() == 0).select(1, mesh_dimensions);
780
781 scale = transform["dimensions"];
782 const int scale_size = scale.size();
783 scale.conservativeResize(dim);
784 if (scale_size < dim)
785 scale.tail(dim - scale_size).setZero();
786
787 scale.array() /= modified_dimensions.array();
788 }
789 else if (transform["scale"].is_number())
790 {
791 scale.setConstant(dim, transform["scale"].get<double>());
792 }
793 else
794 {
795 assert(transform["scale"].is_array());
796 scale = transform["scale"];
797 const int scale_size = scale.size();
798 scale.conservativeResize(dim);
799 if (scale_size < dim)
800 scale.tail(dim - scale_size).setZero();
801
802 if (scale_size == 0)
803 scale.setOnes();
804 }
805
806 A = (unit_scale * scale).asDiagonal();
807
808 // ------
809 // Rotate
810 // ------
811
812 // Rotate around the models origin NOT the bodies center of mass.
813 // We could expose this choice as a "rotate_around" field.
814 MatrixNd R = MatrixNd::Identity(dim, dim);
815 if (!transform["rotation"].is_null())
816 {
817 if (dim == 2)
818 {
819 if (transform["rotation"].is_number())
820 R = Eigen::Rotation2Dd(deg2rad(transform["rotation"].get<double>()))
821 .toRotationMatrix();
822 else if (!transform["rotation"].is_array() || !transform["rotation"].empty())
823 log_and_throw_error("Invalid 2D rotation; 2D rotations can only be a angle in degrees.");
824 }
825 else if (dim == 3)
826 {
827 R = to_rotation_matrix(transform["rotation"], transform["rotation_mode"]);
828 }
829 }
830
831 A = R * A; // Scale first, then rotate
832
833 // ---------
834 // Translate
835 // ---------
836
837 b = transform["translation"];
838 const int translation_size = b.size();
839 b.conservativeResize(dim);
840 if (translation_size < dim)
841 b.tail(dim - translation_size).setZero();
842 }
843
844} // namespace polyfem::mesh
point
Eigen::RowVectorXd point
Definition NavierStokesFSI.cpp:61
faces
std::vector< Eigen::VectorXi > faces
Definition PressureAssembler.cpp:52
polyfem::Units::convert
static double convert(const json &val, const std::string &unit_type)
Definition Units.cpp:35
polyfem::Units::length
const std::string & length() const
Definition Units.hpp:19
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:49
polyfem::mesh::Mesh::n_elements
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
Definition Mesh.hpp:174
polyfem::mesh::Mesh::get_body_id
virtual int get_body_id(const int primitive) const
Get the volume selection of an element (cell in 3d, face in 2d)
Definition Mesh.hpp:525
polyfem::mesh::Mesh::compute_element_barycenters
virtual void compute_element_barycenters(Eigen::MatrixXd &barycenters) const =0
utility for 2d/3d.
polyfem::mesh::Mesh::bounding_box
virtual void bounding_box(RowVectorNd &min, RowVectorNd &max) const =0
computes the bbox of the mesh
polyfem::mesh::Mesh::set_geometry_ids
void set_geometry_ids(const std::vector< int > &geometry_ids)
Set the geometry selection, one ID per element.
Definition Mesh.hpp:540
polyfem::mesh::Mesh::element_vertices
std::vector< int > element_vertices(const int el_id) const
list of vids of an element
Definition Mesh.hpp:242
polyfem::mesh::Mesh::create
static std::unique_ptr< Mesh > create(const std::string &path, const bool non_conforming=false)
factory to build the proper mesh
Definition Mesh.cpp:228
polyfem::mesh::Obstacle::append_mesh
void append_mesh(const Eigen::MatrixXd &vertices, const Eigen::VectorXi &codim_vertices, const Eigen::MatrixXi &codim_edges, const Eigen::MatrixXi &faces, const json &displacement, const std::string &root_path)
Definition Obstacle.cpp:96
polyfem::mesh::Obstacle::append_plane
void append_plane(const VectorNd &point, const VectorNd &normal)
Definition Obstacle.cpp:164
polyfem::mesh::Obstacle::append_mesh_sequence
void append_mesh_sequence(const std::vector< Eigen::MatrixXd > &vertices, const Eigen::VectorXi &codim_vertices, const Eigen::MatrixXi &codim_edges, const Eigen::MatrixXi &faces, const int fps)
Definition Obstacle.cpp:125
polyfem::mesh::Obstacle::set_units
void set_units(const Units &units)
Definition Obstacle.cpp:254
polyfem::utils::Selection::BBox
std::array< RowVectorNd, 2 > BBox
Definition Selection.hpp:13
polyfem::utils::Selection::build_selections
static std::vector< std::shared_ptr< utils::Selection > > build_selections(const json &j_selections, const BBox &mesh_bbox, const std::string &root_path)
Build a vector of selection objects from a JSON selection(s).
Definition Selection.cpp:49
polyfem::mesh::read_obstacle_mesh
void read_obstacle_mesh(const Units &units, const json &j_mesh, const std::string &root_path, const int dim, Eigen::MatrixXd &vertices, Eigen::VectorXi &codim_vertices, Eigen::MatrixXi &codim_edges, Eigen::MatrixXi &faces)
read a obstacle mesh from a geometry JSON
Definition GeometryReader.cpp:401
polyfem::mesh::construct_affine_transformation
void construct_affine_transformation(const double unit_scale, const json &transform, const VectorNd &mesh_dimensions, MatrixNd &A, VectorNd &b)
Construct an affine transformation .
Definition GeometryReader.cpp:762
polyfem::mesh::read_obstacle_geometry
Obstacle read_obstacle_geometry(const Units &units, const json &geometry, const std::vector< json > &displacements, const std::vector< json > &dirichlets, const std::string &root_path, const int dim, const std::vector< std::string > &_names, const std::vector< Eigen::MatrixXd > &_vertices, const std::vector< Eigen::MatrixXi > &_cells, const bool non_conforming)
read a FEM mesh from a geometry JSON
Definition GeometryReader.cpp:512
polyfem::mesh::read_fem_geometry
std::unique_ptr< Mesh > read_fem_geometry(const Units &units, const json &geometry, const std::string &root_path, const std::vector< std::string > &_names, const std::vector< Eigen::MatrixXd > &_vertices, const std::vector< Eigen::MatrixXi > &_cells, const bool non_conforming)
read FEM meshes from a geometry JSON array (or single)
Definition GeometryReader.cpp:304
polyfem::mesh::read_fem_mesh
std::unique_ptr< Mesh > read_fem_mesh(const Units &units, const json &j_mesh, const std::string &root_path, const bool non_conforming)
read a FEM mesh from a geometry JSON
Definition GeometryReader.cpp:64
polyfem::mesh::read_surface_mesh
bool read_surface_mesh(const std::string &mesh_path, Eigen::MatrixXd &vertices, Eigen::VectorXi &codim_vertices, Eigen::MatrixXi &codim_edges, Eigen::MatrixXi &faces)
read a surface mesh
Definition MeshUtils.cpp:1292
polyfem::mesh::apply_geometry_selection
void apply_geometry_selection(Mesh &mesh, const json &geometry_selection, const std::string &root_path)
Apply a geometry selection to a FEM mesh.
Definition GeometryReader.cpp:25
polyfem::utils::resolve_path
std::string resolve_path(const std::string &path, const std::string &input_file_path, const bool only_if_exists=false)
Definition StringUtils.cpp:103
polyfem::utils::to_rotation_matrix
Eigen::Matrix3d to_rotation_matrix(const json &jr, std::string mode)
Definition JSONUtils.cpp:61
polyfem::utils::json_as_array
std::vector< T > json_as_array(const json &j)
Return the value of a json object as an array.
Definition JSONUtils.hpp:38
polyfem::utils::deg2rad
T deg2rad(T deg)
Definition JSONUtils.hpp:18
polyfem::utils::is_param_valid
bool is_param_valid(const json &params, const std::string &key)
Determine if a key exists and is non-null in a json object.
Definition JSONUtils.cpp:131
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
polyfem::VectorNd
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 > VectorNd
Definition Types.hpp:11
polyfem::json
nlohmann::json json
Definition Common.hpp:9
polyfem::MatrixNd
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor, 3, 3 > MatrixNd
Definition Types.hpp:14
polyfem::VectorNi
Eigen::Matrix< int, Eigen::Dynamic, 1, 0, 3, 1 > VectorNi
Definition Types.hpp:12
polyfem::RowVectorNd
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13
polyfem::log_and_throw_error
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73