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