GeometryLoader.cpp

Go to the documentation of this file.

#include "GeometryLoader.hpp"
 
#include <polyfem/mesh/Mesh.hpp>
#include <polyfem/mesh/MeshLoader.hpp>
#include <polyfem/mesh/MeshUtils.hpp>
#include <polyfem/utils/JSONUtils.hpp>
#include <polyfem/utils/Selection.hpp>
#include <polyfem/utils/Logger.hpp>
 
#include <Eigen/Core>
 
#include <igl/edges.h>
#include <igl/boundary_facets.h>
 
#include <strnatcmp.h>
 
namespace polyfem::mesh
{
    using namespace polyfem::utils;
 
    LoadedGeometry GeometryLoader::load(
        const json &geometry,
        const std::vector<json> &obstacle_displacements,
        const std::vector<json> &dirichlet_conditions,
        const bool non_conforming) const
    {
        LoadedGeometry result;
        result.fem = load_fem(geometry, non_conforming);
        if (!result.fem)
            log_and_throw_error("Configured geometry contains no enabled FEM mesh.");
        result.obstacle = load_obstacles(
            geometry, obstacle_displacements, dirichlet_conditions,
            result.fem->dimension(), non_conforming);
        return result;
    }
 
    void GeometryLoader::apply_geometry_selection(
        Mesh &mesh,
        const json &geometry_selection) const
    {
        if (geometry_selection.is_object() && geometry_selection.contains("same_as_volume"))
        {
            if (!geometry_selection["same_as_volume"].get<bool>())
                return;
 
            std::vector<int> geometry_ids(mesh.n_elements());
            for (int e = 0; e < mesh.n_elements(); ++e)
                geometry_ids[e] = mesh.get_body_id(e);
            mesh.set_geometry_ids(geometry_ids);
            return;
        }
 
        Selection::BBox bbox;
        mesh.bounding_box(bbox[0], bbox[1]);
        const auto selections = Selection::build_selections(geometry_selection, bbox, resources_);
 
        Eigen::MatrixXd barycenters;
        mesh.compute_element_barycenters(barycenters);
        std::vector<int> geometry_ids(mesh.n_elements(), 0);
        for (int e = 0; e < mesh.n_elements(); ++e)
        {
            const auto vertices = mesh.element_vertices(e);
            for (const auto &selection : selections)
            {
                if (selection->inside(e, vertices, barycenters.row(e)))
                {
                    geometry_ids[e] = selection->id(e, vertices, barycenters.row(e));
                    break;
                }
            }
        }
        mesh.set_geometry_ids(geometry_ids);
    }
 
    std::unique_ptr<Mesh> GeometryLoader::load_fem_entry(
        const json &j_mesh,
        const bool non_conforming) const
    {
        if (!is_param_valid(j_mesh, "mesh"))
            log_and_throw_error("Mesh {} is mising a \"mesh\" field!", j_mesh);
 
        if (j_mesh["extract"].get<std::string>() != "volume")
            log_and_throw_error("Only volumetric elements are implemented for FEM meshes!");
 
        std::unique_ptr<Mesh> mesh = MeshLoader(resources_).load_fem(j_mesh["mesh"].get<std::string>(), non_conforming);
        // --------------------------------------------------------------------
 
        // NOTE: Normaliziation is done before transformations are applied and/or any selection operators
        if (j_mesh["advanced"]["normalize_mesh"])
            mesh->normalize();
 
        // --------------------------------------------------------------------
 
        Selection::BBox bbox;
        mesh->bounding_box(bbox[0], bbox[1]);
 
        const std::string unit = j_mesh["unit"];
        double unit_scale = 1;
        if (!unit.empty())
            unit_scale = Units::convert(1, unit, units_.length());
 
        {
            MatrixNd A;
            VectorNd b;
            construct_affine_transformation(
                unit_scale,
                j_mesh["transformation"],
                (bbox[1] - bbox[0]).cwiseAbs().transpose(),
                A, b);
            mesh->apply_affine_transformation(A, b);
        }
 
        mesh->bounding_box(bbox[0], bbox[1]);
 
        // --------------------------------------------------------------------
        std::vector<std::shared_ptr<Selection>> surface_selections =
            is_param_valid(j_mesh, "surface_selection") ? Selection::build_selections(j_mesh["surface_selection"], bbox, resources_) : std::vector<std::shared_ptr<Selection>>();
 
        // --------------------------------------------------------------------
 
        const int n_refs = j_mesh["n_refs"];
        const double refinement_location = j_mesh["advanced"]["refinement_location"];
        // TODO: renable this
        // if (n_refs <= 0 && args["poly_bases"] == "MFSHarmonic" && mesh->has_poly())
        // {
        //  if (args["force_no_ref_for_harmonic"])
        //      logger().warn("Using harmonic bases without refinement");
        //  else
        //      n_refs = 1;
        // }
        if (n_refs > 0)
        {
            if (mesh->has_boundary_ids() && surface_selections.empty())
                log_and_throw_error("Unable to refine a mesh with stored surface selections; provide an explicit surface_selection to recompute them after refinement.");
 
            // Check if the stored volume selection is uniform.
            assert(mesh->n_elements() > 0);
            const int uniform_value = mesh->get_body_id(0);
            for (int i = 1; i < mesh->n_elements(); ++i)
                if (mesh->get_body_id(i) != uniform_value)
                    log_and_throw_error("Unable to apply stored nonuniform volume_selection because n_refs={} > 0!", n_refs);
 
            logger().info("Performing global h-refinement with {} refinements", n_refs);
            mesh->refine(n_refs, refinement_location);
            mesh->set_body_ids(std::vector<int>(mesh->n_elements(), uniform_value));
        }
 
        // --------------------------------------------------------------------
 
        if (j_mesh["advanced"]["min_component"].get<int>() != -1)
            log_and_throw_error("Option \"min_component\" in geometry not implement yet!");
        // TODO:
        // if (args["min_component"] > 0) {
        //  Eigen::SparseMatrix<int> adj;
        //  igl::facet_adjacency_matrix(boundary_triangles, adj);
        //  Eigen::MatrixXi C, counts;
        //  igl::connected_components(adj, C, counts);
        //  std::vector<int> valid;
        //  const int min_count = args["min_component"];
        //  for (int i = 0; i < counts.size(); ++i) {
        //      if (counts(i) >= min_count) {
        //          valid.push_back(i);
        //      }
        //  }
        //  tris.clear();
        //  for (int i = 0; i < C.size(); ++i) {
        //      for (int v : valid) {
        //          if (v == C(i)) {
        //              tris.emplace_back(boundary_triangles(i, 0), boundary_triangles(i, 1), boundary_triangles(i, 2));
        //              break;
        //          }
        //      }
        //  }
        //  boundary_triangles.resize(tris.size(), 3);
        //  for (int i = 0; i < tris.size(); ++i) {
        //      boundary_triangles.row(i) << std::get<0>(tris[i]), std::get<1>(tris[i]), std::get<2>(tris[i]);
        //  }
        // }
 
        // --------------------------------------------------------------------
 
        if (j_mesh["advanced"]["force_linear_geometry"].get<bool>())
            log_and_throw_error("Option \"force_linear_geometry\" in geometry not implement yet!");
        // TODO:
        // if (!iso_parametric()) {
        //  if (args["force_linear_geometry"] || mesh->orders().size() <= 0) {
        //      geom_disc_orders.resizeLike(disc_orders);
        //      geom_disc_orders.setConstant(1);
        //  } else {
        //      geom_disc_orders = mesh->orders();
        //  }
        // }
 
        // --------------------------------------------------------------------
 
        const std::vector<std::shared_ptr<Selection>> node_selections =
            is_param_valid(j_mesh, "point_selection") ? Selection::build_selections(j_mesh["point_selection"], bbox, resources_) : std::vector<std::shared_ptr<Selection>>();
 
        if (!node_selections.empty())
        {
            bool boundary_only = true;
            for (const auto &selection : node_selections)
            {
                if (!selection->boundary_only())
                {
                    boundary_only = false;
                    break;
                }
            }
 
            mesh->compute_node_ids([&](const size_t n_id, const RowVectorNd &p, bool is_boundary) {
                if (boundary_only && !is_boundary)
                    return -1;
 
                const std::vector<int> tmp = {int(n_id)};
                for (const auto &selection : node_selections)
                {
                    if (selection->boundary_only() && !is_boundary)
                        continue;
 
                    if (selection->inside(n_id, tmp, p))
                        return selection->id(n_id, tmp, p);
                }
                return std::numeric_limits<int>::max(); // default for no selected boundary
            });
        }
 
        if (!j_mesh["curve_selection"].is_null())
            log_and_throw_error("Geometry curve selections are not implemented!");
 
        // --------------------------------------------------------------------
 
        if (!surface_selections.empty())
        {
            bool boundary_only = true;
            for (const auto &selection : surface_selections)
            {
                if (!selection->boundary_only())
                {
                    boundary_only = false;
                    break;
                }
            }
 
            mesh->compute_boundary_ids([&](const size_t p_id, const std::vector<int> &vs, const RowVectorNd &p, bool is_boundary) {
                if (boundary_only && !is_boundary)
                    return -1;
 
                for (const auto &selection : surface_selections)
                {
                    if (selection->boundary_only() && !is_boundary)
                        continue;
 
                    if (selection->inside(p_id, vs, p))
                        return selection->id(p_id, vs, p);
                }
                return std::numeric_limits<int>::max(); // default for no selected boundary
            });
        }
 
        // --------------------------------------------------------------------
 
        // If the selection is of the form {"id_offset": ...}
        const json volume_selection = j_mesh["volume_selection"];
        if (volume_selection.is_object()
            && volume_selection.size() == 1
            && volume_selection.contains("id_offset"))
        {
            const int id_offset = volume_selection["id_offset"].get<int>();
            if (id_offset != 0)
            {
                const int n_body_ids = mesh->n_elements();
                std::vector<int> body_ids(n_body_ids);
                for (int i = 0; i < n_body_ids; ++i)
                    body_ids[i] = mesh->get_body_id(i) + id_offset;
                mesh->set_body_ids(body_ids);
            }
        }
        else
        {
            // Specified volume selection has priority over mesh's stored ids
            std::vector<std::shared_ptr<Selection>> volume_selections =
                Selection::build_selections(volume_selection, bbox, resources_);
 
            // Append the mesh's stored ids to the volume selection as a lowest priority selection
            if (mesh->has_body_ids())
                volume_selections.push_back(std::make_shared<SpecifiedSelection>(mesh->get_body_ids()));
 
            mesh->compute_body_ids([&](const size_t cell_id, const std::vector<int> &vs, const RowVectorNd &p) -> int {
                for (const auto &selection : volume_selections)
                {
                    // TODO: add vs to compute_body_ids
                    if (selection->inside(cell_id, vs, p))
                        return selection->id(cell_id, vs, p);
                }
                return 0;
            });
        }
 
        // --------------------------------------------------------------------
 
        if (is_param_valid(j_mesh, "geometry_selection"))
            apply_geometry_selection(*mesh, j_mesh["geometry_selection"]);
 
        // --------------------------------------------------------------------
 
        return mesh;
    }
 
    // ========================================================================
 
    std::unique_ptr<Mesh> GeometryLoader::load_fem(
        const json &geometry,
        const bool non_conforming) const
    {
        // --------------------------------------------------------------------
 
        if (geometry.empty())
            log_and_throw_error("Provided geometry is empty!");
 
        std::vector<json> geometries = utils::json_as_array(geometry);
 
        // --------------------------------------------------------------------
 
        std::unique_ptr<Mesh> mesh = nullptr;
 
        for (const json &geometry : geometries)
        {
            if (!geometry["enabled"].get<bool>() || geometry["is_obstacle"].get<bool>())
                continue;
 
            if (geometry["type"] != "mesh" && geometry["type"] != "mesh_array")
                log_and_throw_error("Invalid geometry type \"{}\" for FEM mesh!", geometry["type"]);
 
            const std::unique_ptr<Mesh> tmp_mesh = load_fem_entry(geometry, non_conforming);
 
            if (mesh == nullptr)
                mesh = tmp_mesh->copy();
            else
                mesh->append(tmp_mesh);
 
            if (geometry["type"] == "mesh_array")
            {
                Selection::BBox bbox;
                tmp_mesh->bounding_box(bbox[0], bbox[1]);
 
                const long dim = tmp_mesh->dimension();
                const bool is_offset_relative = geometry["array"]["relative"];
                const double offset = geometry["array"]["offset"];
                const VectorNd dimensions = (bbox[1] - bbox[0]);
                const VectorNi size = geometry["array"]["size"];
 
                for (int i = 0; i < size[0]; ++i)
                {
                    for (int j = 0; j < size[1]; ++j)
                    {
                        for (int k = 0; k < (size.size() > 2 ? size[2] : 1); ++k)
                        {
                            if (i == 0 && j == 0 && k == 0)
                                continue;
 
                            RowVectorNd translation = offset * Eigen::RowVector3d(i, j, k).head(dim);
                            if (is_offset_relative)
                                translation.array() *= dimensions.array();
 
                            const std::unique_ptr<Mesh> copy_mesh = tmp_mesh->copy();
                            copy_mesh->apply_affine_transformation(MatrixNd::Identity(dim, dim), translation);
                            mesh->append(copy_mesh);
                        }
                    }
                }
            }
        }
 
        // --------------------------------------------------------------------
 
        return mesh;
    }
 
    // ========================================================================
 
    SurfaceMesh GeometryLoader::load_surface(
        const json &j_mesh,
        const int dim) const
    {
        if (!is_param_valid(j_mesh, "mesh"))
            log_and_throw_error("Mesh obstacle {} is mising a \"mesh\" field!", j_mesh);
 
        SurfaceMesh result = MeshLoader(resources_).load_surface(j_mesh["mesh"].get<std::string>(), /*keep_volume=*/true);
        auto &vertices = result.vertices;
        auto &codim_vertices = result.points;
        auto &codim_edges = result.edges;
        auto &faces = result.faces;
        auto &tets = result.tets;
 
        const int prev_dim = vertices.cols();
        vertices.conservativeResize(vertices.rows(), dim);
        if (prev_dim < dim)
            vertices.rightCols(dim - prev_dim).setZero();
 
        // --------------------------------------------------------------------
 
        {
            const std::string unit = j_mesh["unit"];
            double unit_scale = 1;
            if (!unit.empty())
                unit_scale = Units::convert(1, unit, units_.length());
 
            const VectorNd mesh_dimensions = (vertices.colwise().maxCoeff() - vertices.colwise().minCoeff()).cwiseAbs();
            MatrixNd A;
            VectorNd b;
            construct_affine_transformation(unit_scale, j_mesh["transformation"], mesh_dimensions, A, b);
            vertices = vertices * A.transpose();
            vertices.rowwise() += b.transpose();
        }
 
        // Preserve volumetric (2D triangle) connectivity for optional volumetric paraview export.
        // In 3D a tet mesh's tets come from load_surface with keep_volume.
        if (dim == 2 && tets.rows() == 0 && faces.rows() > 0 && faces.cols() == 3)
            tets = faces;
 
        std::string extract = j_mesh["extract"];
        // Default: "volume" clashes with defaults for non obstacle, here assume volume is suface
        if (extract == "volume")
            extract = "surface";
 
        if (extract == "points")
        {
            // points -> vertices (drop edges and faces)
            codim_edges.resize(0, 0);
            faces.resize(0, 0);
            codim_vertices.resize(vertices.rows());
            for (int i = 0; i < codim_vertices.size(); ++i)
                codim_vertices[i] = i;
        }
        else if (extract == "edges" && faces.size() != 0)
        {
            // edges -> edges (drop faces)
            Eigen::MatrixXi edges;
            igl::edges(faces, edges);
            faces.resize(0, 0);
            codim_edges.conservativeResize(codim_edges.rows() + edges.rows(), 2);
            codim_edges.bottomRows(edges.rows()) = edges;
        }
        else if (extract == "surface" && dim == 2 && faces.size() != 0)
        {
            // surface (2D) -> boundary edges (drop faces and interior edges)
            Eigen::MatrixXi boundary_edges;
            igl::boundary_facets(faces, boundary_edges);
            codim_edges.conservativeResize(codim_edges.rows() + boundary_edges.rows(), 2);
            codim_edges.bottomRows(boundary_edges.rows()) = boundary_edges;
            faces.resize(0, 0); // Clear faces
        }
        // surface (3D) -> boundary faces
        // Surface extraction in 3D already has the desired faces from MeshLoader.
        else if (extract == "volume")
        {
            // volume -> undefined
            log_and_throw_error("Volumetric elements not supported for collision obstacles!");
        }
 
        if (j_mesh["n_refs"].get<int>() != 0)
        {
            if (faces.size() != 0)
                log_and_throw_error("Option \"n_refs\" for triangle obstacles not implement yet!");
 
            const int n_refs = j_mesh["n_refs"];
            const double refinement_location = j_mesh["advanced"]["refinement_location"];
            for (int i = 0; i < n_refs; i++)
            {
                const size_t n_vertices = vertices.rows();
                const size_t n_edges = codim_edges.rows();
                vertices.conservativeResize(n_vertices + n_edges, vertices.cols());
                codim_edges.conservativeResize(2 * n_edges, codim_edges.cols());
                for (size_t ei = 0; ei < n_edges; ei++)
                {
                    const int v0i = codim_edges(ei, 0);
                    const int v1i = codim_edges(ei, 1);
                    const int v2i = n_vertices + ei;
                    vertices.row(v2i) = (vertices.row(v1i) - vertices.row(v0i)) * refinement_location + vertices.row(v0i);
                    codim_edges.row(ei) << v0i, v2i;
                    codim_edges.row(n_edges + ei) << v2i, v1i;
                }
            }
        }
        return result;
    }
 
    // ========================================================================
 
    Obstacle GeometryLoader::load_obstacles(
        const json &geometry,
        const std::vector<json> &displacements,
        const std::vector<json> &dirichlets,
        const int dim,
        const bool) const
    {
        Obstacle obstacle;
 
        if (geometry.empty())
            return obstacle;
 
        std::vector<json> geometries = utils::json_as_array(geometry);
 
        for (const json &geometry : geometries)
        {
 
            if (!geometry["is_obstacle"].get<bool>())
                continue;
 
            if (!geometry["enabled"].get<bool>())
                continue;
 
            if (geometry["type"] == "mesh" || geometry["type"] == "mesh_array")
            {
                SurfaceMesh loaded = load_surface(geometry, dim);
                Eigen::MatrixXd &vertices = loaded.vertices;
                Eigen::VectorXi &codim_vertices = loaded.points;
                Eigen::MatrixXi &codim_edges = loaded.edges;
                Eigen::MatrixXi &faces = loaded.faces;
                Eigen::MatrixXi &tets = loaded.tets;
 
                if (geometry["type"] == "mesh_array")
                {
                    const Selection::BBox bbox{{vertices.colwise().minCoeff(), vertices.colwise().maxCoeff()}};
 
                    const bool is_offset_relative = geometry["array"]["relative"];
                    const double offset = geometry["array"]["offset"];
                    const VectorNd dimensions = (bbox[1] - bbox[0]);
                    const VectorNi size = geometry["array"]["size"];
 
                    const int N = size.head(dim).prod();
                    const int nV = vertices.rows(), nCV = codim_vertices.rows(), nCE = codim_edges.rows(), nF = faces.rows(), nT = tets.rows();
 
                    vertices.conservativeResize(N * nV, Eigen::NoChange);
                    codim_vertices.conservativeResize(N * nCV, Eigen::NoChange);
                    codim_edges.conservativeResize(N * nCE, Eigen::NoChange);
                    faces.conservativeResize(N * nF, Eigen::NoChange);
                    tets.conservativeResize(N * nT, tets.cols());
 
                    for (int i = 0; i < size[0]; ++i)
                    {
                        for (int j = 0; j < size[1]; ++j)
                        {
                            for (int k = 0; k < (size.size() > 2 ? size[2] : 1); ++k)
                            {
                                RowVectorNd translation = offset * Eigen::RowVector3d(i, j, k).head(vertices.cols());
                                if (is_offset_relative)
                                    translation.array() *= dimensions.array();
 
                                int n = i * size[1] + j;
                                if (size.size() > 2)
                                    n = n * size[2] + k;
                                if (n == 0)
                                    continue;
 
                                vertices.middleRows(n * nV, nV) = vertices.topRows(nV).rowwise() + translation;
                                if (nCV)
                                    codim_vertices.segment(n * nV, nV) = codim_vertices.head(nV).array() + n * nV;
                                if (nCE)
                                    codim_edges.middleRows(n * nCE, nCE) = codim_edges.topRows(nCE).array() + n * nV;
                                if (nF)
                                    faces.middleRows(n * nF, nF) = faces.topRows(nF).array() + n * nV;
                                if (nT)
                                    tets.middleRows(n * nT, nT) = tets.topRows(nT).array() + n * nV;
                            }
                        }
                    }
                }
 
                json displacement = "{\"value\":[0, 0, 0]}"_json;
                if (is_param_valid(geometry, "surface_selection"))
                {
                    if (!geometry["surface_selection"].is_number())
                        log_and_throw_error("Invalid surface_selection for obstacle, needs to be an integer!");
 
                    const int id = geometry["surface_selection"];
                    for (const json &disp : dirichlets)
                    {
                        if ((disp["id"].is_string() && disp["id"].get<std::string>() == "all")
                            || (disp["id"].is_number_integer() && disp["id"].get<int>() == id))
                        {
                            displacement = disp;
                            break;
                        }
                        else if (disp["id"].is_array())
                        {
                            for (const json &disp_id : disp["id"])
                            {
                                assert(disp_id.is_number_integer());
                                if (disp_id.get<int>() == id)
                                {
                                    displacement = disp;
                                    break;
                                }
                            }
                        }
                    }
                    for (const json &disp : displacements)
                    {
                        if ((disp["id"].is_string() && disp["id"].get<std::string>() == "all")
                            || (disp["id"].is_number_integer() && disp["id"].get<int>() == id))
                        {
                            displacement = disp;
                            break;
                        }
                        else if (disp["id"].is_array())
                        {
                            for (const json &disp_id : disp["id"])
                            {
                                assert(disp_id.is_number_integer());
                                if (disp_id.get<int>() == id)
                                {
                                    displacement = disp;
                                    break;
                                }
                            }
                        }
                    }
                }
 
                obstacle.append_mesh(
                    vertices, codim_vertices, codim_edges, faces, displacement, resources_, tets);
            }
            else if (geometry["type"] == "plane")
            {
                obstacle.append_plane(geometry["point"], geometry["normal"]);
            }
            else if (geometry["type"] == "ground")
            {
                VectorNd gravity = VectorNd::Zero(dim); // TODO: Expose as parameter
                gravity[1] = -9.81;
                const double height = geometry["height"];
                assert(gravity.norm() != 0);
                const VectorNd normal = -gravity.normalized();
                const VectorNd point = height * normal; // origin + height * normal
                obstacle.append_plane(point, normal);
            }
            else if (geometry["type"] == "mesh_sequence")
            {
                std::vector<std::string> mesh_files;
                if (geometry["mesh_sequence"].is_array())
                {
                    mesh_files = geometry["mesh_sequence"].get<std::vector<std::string>>();
                }
                else
                {
                    assert(geometry["mesh_sequence"].is_string());
                    const std::string sequence = geometry["mesh_sequence"].get<std::string>();
                    mesh_files = resources_.is_group(sequence) ? resources_.list(sequence) : resources_.glob(sequence);
                    // Sort the file names naturally
                    std::sort(mesh_files.begin(), mesh_files.end(), [](const std::string &p1, const std::string &p2) {
                        return strnatcmp(p1.c_str(), p2.c_str()) < 0;
                    });
                }
 
                std::vector<Eigen::MatrixXd> vertices(mesh_files.size());
                Eigen::VectorXi codim_vertices;
                Eigen::MatrixXi codim_edges;
                Eigen::MatrixXi faces;
                Eigen::MatrixXi tets;
 
                for (int i = 0; i < mesh_files.size(); ++i)
                {
                    json jmesh = geometry;
                    jmesh["mesh"] = mesh_files[i];
                    jmesh["n_refs"] = 0;
 
                    SurfaceMesh loaded = load_surface(jmesh, dim);
                    vertices[i] = std::move(loaded.vertices);
                    Eigen::VectorXi &tmp_codim_vertices = loaded.points;
                    Eigen::MatrixXi &tmp_codim_edges = loaded.edges;
                    Eigen::MatrixXi &tmp_faces = loaded.faces;
                    Eigen::MatrixXi &tmp_tets = loaded.tets;
                    if (i == 0)
                    {
                        codim_vertices = tmp_codim_vertices;
                        codim_edges = tmp_codim_edges;
                        faces = tmp_faces;
                        tets = tmp_tets;
                    }
                    else
                    {
                        assert((codim_vertices.array() == tmp_codim_vertices.array()).all());
                        assert((codim_edges.array() == tmp_codim_edges.array()).all());
                        assert((faces.array() == tmp_faces.array()).all());
                        assert((tets.array() == tmp_tets.array()).all());
                    }
                }
 
                obstacle.append_mesh_sequence(
                    vertices, codim_vertices, codim_edges, faces, geometry["fps"], tets);
            }
            else
            {
                log_and_throw_error("Invalid geometry type \"{}\" for obstacle!", geometry["type"]);
            }
        }
 
        obstacle.set_units(units_);
        return obstacle;
    }
 
    // ========================================================================
 
    void GeometryLoader::construct_affine_transformation(
        const double unit_scale,
        const json &transform,
        const VectorNd &mesh_dimensions,
        MatrixNd &A,
        VectorNd &b)
    {
        const int dim = mesh_dimensions.size();
 
        // -----
        // Scale
        // -----
 
        RowVectorNd scale;
        if (transform["dimensions"].is_array()) // default is nullptr
        {
            VectorNd modified_dimensions =
                (mesh_dimensions.array() == 0).select(1, mesh_dimensions);
 
            scale = transform["dimensions"];
            const int scale_size = scale.size();
            scale.conservativeResize(dim);
            if (scale_size < dim)
                scale.tail(dim - scale_size).setZero();
 
            scale.array() /= modified_dimensions.array();
        }
        else if (transform["scale"].is_number())
        {
            scale.setConstant(dim, transform["scale"].get<double>());
        }
        else
        {
            assert(transform["scale"].is_array());
            scale = transform["scale"];
            const int scale_size = scale.size();
            scale.conservativeResize(dim);
            if (scale_size < dim)
                scale.tail(dim - scale_size).setZero();
 
            if (scale_size == 0)
                scale.setOnes();
        }
 
        A = (unit_scale * scale).asDiagonal();
 
        // ------
        // Rotate
        // ------
 
        // Rotate around the models origin NOT the bodies center of mass.
        // We could expose this choice as a "rotate_around" field.
        MatrixNd R = MatrixNd::Identity(dim, dim);
        if (!transform["rotation"].is_null())
        {
            if (dim == 2)
            {
                if (transform["rotation"].is_number())
                    R = Eigen::Rotation2Dd(deg2rad(transform["rotation"].get<double>()))
                            .toRotationMatrix();
                else if (!transform["rotation"].is_array() || !transform["rotation"].empty())
                    log_and_throw_error("Invalid 2D rotation; 2D rotations can only be a angle in degrees.");
            }
            else if (dim == 3)
            {
                R = to_rotation_matrix(transform["rotation"], transform["rotation_mode"]);
            }
        }
 
        A = R * A; // Scale first, then rotate
 
        // ---------
        // Translate
        // ---------
 
        b = transform["translation"];
        const int translation_size = b.size();
        b.conservativeResize(dim);
        if (translation_size < dim)
            b.tail(dim - translation_size).setZero();
    }
 
} // namespace polyfem::mesh