Remesher.hpp

Go to the documentation of this file.

#pragma once
 
#include <polyfem/legacy/State.hpp>
#include <polyfem/utils/Types.hpp>
#include <polyfem/utils/Timer.hpp>
 
#include <unordered_map>
#include <variant>
 
namespace polyfem::time_integrator
{
    class ImplicitTimeIntegrator;
} // namespace polyfem::time_integrator
 
#define POLYFEM_REMESHER_SCOPED_TIMER(name) polyfem::utils::Timer __polyfem_timer(Remesher::timings[name])
 
namespace polyfem::mesh
{
    using State = polyfem::legacy::State;
 
    class Remesher
    {
        // --------------------------------------------------------------------
        // typedefs
    public:
        template <typename T>
        using EdgeMap = std::unordered_map<
            std::array<size_t, 2>, T,
            polyfem::utils::HashUnorderedArray<size_t, 2>,
            polyfem::utils::EqualUnorderedArray<size_t, 2>>;
        template <typename T>
        using FaceMap = std::unordered_map<
            std::array<size_t, 3>, T,
            polyfem::utils::HashUnorderedArray<size_t, 3>,
            polyfem::utils::EqualUnorderedArray<size_t, 3>>;
        template <typename T>
        using TetMap = std::unordered_map<
            std::array<size_t, 4>, T,
            polyfem::utils::HashUnorderedArray<size_t, 4>,
            polyfem::utils::EqualUnorderedArray<size_t, 4>>;
        template <typename T>
        using BoundaryMap = std::variant<EdgeMap<T>, FaceMap<T>>;
 
        // --------------------------------------------------------------------
        // constructors
    public:
        Remesher(const State &state,
                 const Eigen::MatrixXd &obstacle_displacements,
                 const Eigen::MatrixXd &obstacle_quantities,
                 const double current_time,
                 const double starting_energy);
 
        virtual ~Remesher() = default;
 
        virtual void init(
            const Eigen::MatrixXd &rest_positions,
            const Eigen::MatrixXd &positions,
            const Eigen::MatrixXi &elements,
            const Eigen::MatrixXd &projection_quantities,
            const BoundaryMap<int> &boundary_to_id,
            const std::vector<int> &body_ids,
            const EdgeMap<double> &elastic_energy,
            const EdgeMap<double> &contact_energy);
 
    protected:
        virtual void init_attributes_and_connectivity(
            const size_t num_vertices,
            const Eigen::MatrixXi &elements) = 0;
 
        // --------------------------------------------------------------------
        // main functions
    public:
        virtual bool execute() = 0;
 
    protected:
        void project_quantities();
 
        void cache_before();
 
        // --------------------------------------------------------------------
        // getters
    public:
        virtual int dim() const = 0;
 
        virtual bool is_volume() const { return dim() == 3; }
 
        virtual Eigen::MatrixXd rest_positions() const = 0;
        virtual Eigen::MatrixXd displacements() const = 0;
        virtual Eigen::MatrixXd positions() const = 0;
        virtual Eigen::MatrixXi edges() const = 0;
        virtual Eigen::MatrixXi elements() const = 0;
        virtual Eigen::MatrixXi boundary_edges() const = 0;
        virtual Eigen::MatrixXi boundary_faces() const = 0;
        virtual Eigen::MatrixXd projection_quantities() const = 0;
        virtual BoundaryMap<int> boundary_ids() const = 0;
        virtual std::vector<int> body_ids() const = 0;
        virtual std::vector<int> boundary_nodes(const Eigen::VectorXi &vertex_to_basis) const = 0;
 
        virtual int n_quantities() const = 0;
 
        const Obstacle &obstacle() const { return state.obstacle; }
        const Eigen::MatrixXd &obstacle_displacements() const { return m_obstacle_displacements; }
        const Eigen::MatrixXd &obstacle_quantities() const { return m_obstacle_quantities; }
 
        // --------------------------------------------------------------------
        // setters
    public:
        virtual void set_rest_positions(const Eigen::MatrixXd &positions) = 0;
        virtual void set_positions(const Eigen::MatrixXd &positions) = 0;
        virtual void set_projection_quantities(const Eigen::MatrixXd &projection_quantities) = 0;
        virtual void set_fixed(const std::vector<bool> &fixed) = 0;
        virtual void set_boundary_ids(const BoundaryMap<int> &boundary_to_id) = 0;
        virtual void set_body_ids(const std::vector<int> &body_ids) = 0;
 
        // --------------------------------------------------------------------
        // utilities
    public:
        void write_mesh(const std::string &path) const;
 
        static Eigen::MatrixXd combine_time_integrator_quantities(
            const std::shared_ptr<time_integrator::ImplicitTimeIntegrator> &time_integrator);
 
        static void split_time_integrator_quantities(
            const Eigen::MatrixXd &quantities,
            const int dim,
            Eigen::MatrixXd &x_prevs,
            Eigen::MatrixXd &v_prevs,
            Eigen::MatrixXd &a_prevs);
 
        void init_assembler(const std::vector<int> &body_ids) const;
 
        static int build_bases(
            const Mesh &mesh,
            const std::string &assembler_formulation,
            std::vector<polyfem::basis::ElementBases> &bases,
            std::vector<LocalBoundary> &local_boundary,
            Eigen::VectorXi &vertex_to_basis);
 
        const State &state;
 
        // --------------------------------------------------------------------
        // members
    public:
        int max_op_attempts = 1;
 
    protected:
        // TODO: Drop this and only use a local EdgeOperationCache
        struct GlobalProjectionCache
        {
            Eigen::MatrixXd rest_positions;
            Eigen::MatrixXi elements;
            Eigen::MatrixXd projection_quantities;
        };
 
        GlobalProjectionCache global_projection_cache;
 
        const json args;
        const Eigen::MatrixXd m_obstacle_displacements;
        Eigen::MatrixXd m_obstacle_quantities;
        const double current_time;
        const double starting_energy;
 
        // --------------------------------------------------------------------
        // statistics
    public:
        static void log_timings();
 
        static std::unordered_map<std::string, utils::Timing> timings;
        static double total_time;  // = 0;
        static size_t num_solves;  // = 0;
        static size_t total_ndofs; // = 0;
    };
 
} // namespace polyfem::mesh