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