StateLoad.cpp

Go to the documentation of this file.

#include <polyfem/legacy/State.hpp>
 
#include <polyfem/assembler/Mass.hpp>
#include <polyfem/io/ResourceIO.hpp>
 
#include <polyfem/mesh/GeometryLoader.hpp>
#include <polyfem/mesh/mesh2D/CMesh2D.hpp>
#include <polyfem/mesh/mesh2D/NCMesh2D.hpp>
#include <polyfem/mesh/mesh3D/CMesh3D.hpp>
#include <polyfem/mesh/mesh3D/NCMesh3D.hpp>
#include <polyfem/mesh/MeshReader.hpp>
 
#include <polyfem/utils/Selection.hpp>
 
#include <polyfem/utils/JSONUtils.hpp>
 
#include <igl/Timer.h>
namespace polyfem::legacy
{
    using namespace basis;
    using namespace mesh;
    using namespace utils;
 
    void State::reset_mesh()
    {
        bases.clear();
        pressure_bases.clear();
        geom_bases_.clear();
        boundary_nodes.clear();
        local_boundary.clear();
        local_neumann_boundary.clear();
        polys.clear();
        poly_edge_to_data.clear();
        obstacle.clear();
 
        mass.resize(0, 0);
        rhs.resize(0, 0);
 
        n_bases = 0;
        n_pressure_bases = 0;
    }
 
    void State::load_mesh(GEO::Mesh &meshin, const std::function<int(const size_t, const std::vector<int> &, const RowVectorNd &, bool)> &boundary_marker, bool non_conforming, bool skip_boundary_sideset)
    {
        reset_mesh();
 
        igl::Timer timer;
        timer.start();
        logger().info("Loading mesh...");
        mesh = Mesh::create(mesh::MeshReader::from_geogram(meshin), non_conforming);
        if (!mesh)
        {
            logger().error("Unable to load the mesh");
            return;
        }
 
        RowVectorNd min, max;
        mesh->bounding_box(min, max);
 
        logger().info("mesh bb min [{}], max [{}]", min, max);
 
        if (!skip_boundary_sideset)
            mesh->compute_boundary_ids(boundary_marker);
 
        std::vector<std::shared_ptr<assembler::Assembler>> assemblers;
        assemblers.push_back(assembler);
        assemblers.push_back(mass_matrix_assembler);
        if (mixed_assembler != nullptr)
            // TODO: assemblers.push_back(mixed_assembler);
            mixed_assembler->set_size(mesh->dimension());
        if (pressure_assembler != nullptr)
            assemblers.push_back(pressure_assembler);
        set_materials(assemblers);
 
        timer.stop();
        logger().info(" took {}s", timer.getElapsedTime());
 
        timer.start();
        logger().info("Loading obstacles...");
        const polyfem::io::FileSystemIO resources(args["root_path"].get<std::string>());
        const mesh::GeometryLoader geometry_loader(units, resources);
        obstacle = geometry_loader.load_obstacles(
            args["geometry"],
            utils::json_as_array(args["boundary_conditions"]["obstacle_displacements"]),
            utils::json_as_array(args["boundary_conditions"]["dirichlet_boundary"]),
            mesh->dimension());
 
        timer.stop();
        logger().info(" took {}s", timer.getElapsedTime());
 
        out_geom.init_sampler(*mesh, args["output"]["paraview"]["vismesh_rel_area"]);
    }
 
    void State::load_mesh(bool non_conforming,
                          const std::vector<std::string> &names,
                          const std::vector<Eigen::MatrixXi> &cells,
                          const std::vector<Eigen::MatrixXd> &vertices)
    {
        assert(names.size() == cells.size());
        assert(vertices.size() == cells.size());
        assert(names.empty());
        assert(cells.empty());
        assert(vertices.empty());
 
        reset_mesh();
 
        igl::Timer timer;
        timer.start();
 
        logger().info("Loading mesh ...");
        const polyfem::io::FileSystemIO resources(args["root_path"].get<std::string>());
        const mesh::GeometryLoader geometry_loader(units, resources);
        if (mesh == nullptr)
        {
            assert(is_param_valid(args, "geometry"));
            mesh = geometry_loader.load_fem(args["geometry"], non_conforming);
        }
 
        if (mesh == nullptr)
        {
            log_and_throw_error("unable to load the mesh!");
        }
 
        // if(!flipped_elements.empty())
        // {
        //  mesh->compute_elements_tag();
        //  for(auto el_id : flipped_elements)
        //      mesh->set_tag(el_id, ElementType::INTERIOR_POLYTOPE);
        // }
 
        RowVectorNd min, max;
        mesh->bounding_box(min, max);
 
        if (args["solver"]["advanced"]["characteristic_length"] > 0)
        {
            characteristic_length = args["solver"]["advanced"]["characteristic_length"];
        }
        else
        {
            characteristic_length = (max - min).norm();
        }
 
        logger().info("mesh bb min [{}], max [{}]", min, max);
 
        std::vector<std::shared_ptr<assembler::Assembler>> assemblers;
        assemblers.push_back(assembler);
        assemblers.push_back(mass_matrix_assembler);
        if (mixed_assembler != nullptr)
            // TODO: assemblers.push_back(mixed_assembler);
            mixed_assembler->set_size(mesh->dimension());
        if (pressure_assembler != nullptr)
            assemblers.push_back(pressure_assembler);
        set_materials(assemblers);
        pure_mass_matrix_assembler->set_size(mass_matrix_assembler->size());
 
        timer.stop();
        logger().info(" took {}s", timer.getElapsedTime());
 
        out_geom.init_sampler(*mesh, args["output"]["paraview"]["vismesh_rel_area"]);
 
        timer.start();
        logger().info("Loading obstacles...");
        obstacle = geometry_loader.load_obstacles(
            args["geometry"],
            utils::json_as_array(args["boundary_conditions"]["obstacle_displacements"]),
            utils::json_as_array(args["boundary_conditions"]["dirichlet_boundary"]),
            mesh->dimension(), non_conforming);
        timer.stop();
        logger().info(" took {}s", timer.getElapsedTime());
    }
 
    void State::build_mesh_matrices(Eigen::MatrixXd &V, Eigen::MatrixXi &F)
    {
        assert(bases.size() == mesh->n_elements());
        const size_t n_vertices = n_bases - obstacle.n_vertices();
        const int dim = mesh->dimension();
 
        V.resize(n_vertices, dim);
        F.resize(bases.size(), dim + 1); // TODO: this only works for triangles and tetrahedra
 
        for (int i = 0; i < bases.size(); i++)
        {
            const basis::ElementBases &element = bases[i];
            for (int j = 0; j < element.bases.size(); j++)
            {
                const basis::Basis &basis = element.bases[j];
                assert(basis.global().size() == 1);
                V.row(basis.global()[0].index) = basis.global()[0].node;
                if (j < F.cols()) // Only grab the corners of the triangles/tetrahedra
                    F(i, j) = basis.global()[0].index;
            }
        }
    }
 
    std::unordered_map<int, std::array<bool, 3>>
    State::boundary_conditions_ids(const std::string &bc_type) const
    {
        assert(args["boundary_conditions"].contains(bc_type));
        const std::vector<json> json_bcs = json_as_array(args["boundary_conditions"][bc_type]);
        std::unordered_map<int, std::array<bool, 3>> bcs;
        for (const json &bc : json_bcs)
        {
            assert(bc["dimension"].size() >= mesh->dimension());
            std::array<bool, 3> dimension{{true, true, true}};
            for (int d = 0; d < bc["dimension"].size(); ++d)
                dimension[d] = bc["dimension"][d];
 
            assert(bc.contains("id") && bc["id"].is_number_integer());
            bcs[bc["id"].get<int>()] = dimension;
        }
 
        return bcs;
    }
} // namespace polyfem::legacy