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