WildRemesher.hpp

Go to the documentation of this file.

#pragma once
 
#include <polyfem/mesh/remesh/Remesher.hpp>
#include <polyfem/solver/SolveData.hpp>
 
#include <wmtk/TriMesh.h>
#include <wmtk/TetMesh.h>
#include <wmtk/ExecutionScheduler.hpp>
 
#include <type_traits>
 
namespace polyfem::mesh
{
    enum class CollapseEdgeTo
    {
        V0,
        V1,
        MIDPOINT,
        ILLEGAL
    };
 
    class TriOperationCache;
    class TetOperationCache;
 
    template <class WMTKMesh>
    class WildRemesher : public Remesher, public WMTKMesh
    {
        using This = WildRemesher<WMTKMesh>;
 
        // --------------------------------------------------------------------
        // typedefs
    public:
        // NOTE: This assumes triangle meshes are only used in 2D.
        static constexpr int DIM = [] {
            if constexpr (std::is_same_v<wmtk::TriMesh, WMTKMesh>)
                return 2;
            else
                return 3;
        }();
 
        static constexpr int VERTICES_PER_ELEMENT = [] {
            if constexpr (std::is_same_v<wmtk::TriMesh, WMTKMesh>)
                return 3;
            else
                return 4;
        }();
 
        static constexpr int EDGES_PER_ELEMENT = [] {
            if constexpr (std::is_same_v<wmtk::TriMesh, WMTKMesh>)
                return 3;
            else
                return 6;
        }();
 
        static constexpr int FACETS_PER_ELEMENT = [] {
            if constexpr (std::is_same_v<wmtk::TriMesh, WMTKMesh>)
                return 3;
            else
                return 4;
        }();
 
        using VectorNd = Eigen::Matrix<double, DIM, 1>;
 
        using Tuple = typename WMTKMesh::Tuple;
 
        static constexpr wmtk::ExecutionPolicy EXECUTION_POLICY = wmtk::ExecutionPolicy::kSeq;
 
        // --------------------------------------------------------------------
        // constructors
    public:
        WildRemesher(
            const State &state,
            const Eigen::MatrixXd &obstacle_displacements,
            const Eigen::MatrixXd &obstacle_vals,
            const double current_time,
            const double starting_energy);
 
        virtual ~WildRemesher() = 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) override;
 
    protected:
        void init_attributes_and_connectivity(
            const size_t num_vertices,
            const Eigen::MatrixXi &elements) override;
 
        // --------------------------------------------------------------------
        // main functions
    public:
        bool execute() override;
 
        virtual void split_edges() = 0;
        virtual void collapse_edges() = 0;
        virtual void smooth_vertices();
        virtual void swap_edges() { log_and_throw_error("WildRemesher::swap_edges not implemented!"); }
 
    protected:
        // Edge splitting
        virtual bool split_edge_before(const Tuple &t) override;
        virtual bool split_edge_after(const Tuple &t) override;
 
        // Edge collapse
        virtual bool collapse_edge_before(const Tuple &t) override;
        virtual bool collapse_edge_after(const Tuple &t) override;
 
        // Swap edge
        virtual bool swap_edge_before(const Tuple &t) override;
        virtual bool swap_edge_after(const Tuple &t) override;
 
        // Smooth_vertex
        virtual bool smooth_before(const Tuple &t) override;
        virtual bool smooth_after(const Tuple &t) override;
 
        bool invariants(const std::vector<Tuple> &new_tris) override;
 
        bool is_rest_inverted(const Tuple &loc) const;
        bool is_inverted(const Tuple &loc) const;
 
        // --------------------------------------------------------------------
        // getters
    public:
        int dim() const override { return DIM; }
 
        Eigen::MatrixXd rest_positions() const override;
        Eigen::MatrixXd displacements() const override;
        Eigen::MatrixXd positions() const override;
        Eigen::MatrixXi edges() const override;
        Eigen::MatrixXi elements() const override;
        Eigen::MatrixXi boundary_edges() const override;
        Eigen::MatrixXi boundary_faces() const override;
        Eigen::MatrixXd projection_quantities() const override;
        BoundaryMap<int> boundary_ids() const override;
        std::vector<int> body_ids() const override;
        std::vector<int> boundary_nodes(const Eigen::VectorXi &vertex_to_basis) const override;
 
        int n_quantities() const override { return m_n_quantities; };
 
        std::vector<Tuple> get_facets() const;
        std::vector<Tuple> get_elements() const;
 
        // --------------------------------------------------------------------
        // setters
    public:
        void set_rest_positions(const Eigen::MatrixXd &positions) override;
        void set_positions(const Eigen::MatrixXd &positions) override;
        void set_projection_quantities(const Eigen::MatrixXd &projection_quantities) override;
        void set_fixed(const std::vector<bool> &fixed) override;
        void set_boundary_ids(const BoundaryMap<int> &boundary_to_id) override;
        void set_body_ids(const std::vector<int> &body_ids) override;
 
        // --------------------------------------------------------------------
        // utilities
    public:
        double rest_edge_length(const Tuple &e) const;
        double deformed_edge_length(const Tuple &e) const;
 
        VectorNd rest_edge_center(const Tuple &e) const;
        VectorNd deformed_edge_center(const Tuple &e) const;
 
        Eigen::VectorXd edge_adjacent_element_volumes(const Tuple &e) const;
 
        double element_volume(const Tuple &e) const;
 
        bool is_boundary_vertex(const Tuple &v) const;
        bool is_body_boundary_vertex(const Tuple &v) const;
        bool is_boundary_edge(const Tuple &e) const;
        bool is_body_boundary_edge(const Tuple &e) const;
        bool is_boundary_facet(const Tuple &t) const;
        bool is_boundary_op() const;
 
        std::vector<Tuple> boundary_facets(std::vector<int> *boundary_ids = nullptr) const;
 
        std::array<Tuple, DIM> facet_vertices(const Tuple &t) const;
        std::array<size_t, DIM> facet_vids(const Tuple &t) const;
 
        std::array<Tuple, VERTICES_PER_ELEMENT> element_vertices(const Tuple &t) const;
        std::array<size_t, VERTICES_PER_ELEMENT> element_vids(const Tuple &t) const;
 
        void element_aabb(const Tuple &t, polyfem::VectorNd &el_min, polyfem::VectorNd &el_max) const;
 
        std::array<size_t, VERTICES_PER_ELEMENT> orient_preserve_element_reorder(
            const std::array<size_t, VERTICES_PER_ELEMENT> &conn, const size_t v0) const;
 
        std::vector<Tuple> get_one_ring_elements_for_vertex(const Tuple &t) const;
        std::vector<Tuple> get_one_ring_boundary_edges_for_vertex(const Tuple &v) const;
        std::array<Tuple, 2> get_boundary_faces_for_edge(const Tuple &e) const;
        std::vector<Tuple> get_one_ring_boundary_faces_for_vertex(const Tuple &v) const;
        std::vector<Tuple> get_edges_for_elements(const std::vector<Tuple> &elements) const;
 
        size_t facet_id(const Tuple &t) const;
 
        size_t element_id(const Tuple &t) const;
 
        Tuple tuple_from_element(size_t elem_id) const;
        Tuple tuple_from_facet(size_t elem_id, int local_facet_id) const;
 
        std::vector<Tuple> get_incident_elements_for_edge(const Tuple &t) const;
 
        void extend_local_patch(std::vector<Tuple> &patch) const;
 
        Tuple opposite_vertex_on_face(const Tuple &e) const
        {
            return e.switch_edge(*this).switch_vertex(*this);
        }
 
        CollapseEdgeTo collapse_boundary_edge_to(const Tuple &e) const;
 
    protected:
        using Operations = std::vector<std::pair<std::string, Tuple>>;
        virtual Operations renew_neighbor_tuples(
            const std::string &op, const std::vector<Tuple> &tris) const { return {}; }
 
        void cache_split_edge(const Tuple &e);
 
        void cache_collapse_edge(const Tuple &e, const CollapseEdgeTo collapse_to);
 
        void cache_swap_edge(const Tuple &e);
 
        // NOTE: Nothing to cache for vertex smoothing
 
        void map_edge_split_edge_attributes(const Tuple &t);
        void map_edge_split_boundary_attributes(const Tuple &t);
        void map_edge_split_element_attributes(const Tuple &t);
 
        void map_edge_collapse_vertex_attributes(const Tuple &t);
        void map_edge_collapse_boundary_attributes(const Tuple &t);
        void map_edge_collapse_edge_attributes(const Tuple &t);
 
        void map_edge_swap_edge_attributes(const Tuple &t);
        void map_edge_swap_element_attributes(const Tuple &t);
 
        // NOTE: Nothing to map for vertex smoothing
 
        // --------------------------------------------------------------------
        // members
    public:
        struct VertexAttributes
        {
            using VectorNd = WildRemesher<WMTKMesh>::VectorNd;
 
            VectorNd rest_position;
            VectorNd position;
 
            Eigen::MatrixXd projection_quantities;
 
            // Can the point be smoothed or moved around by operations?
            bool fixed = false;
            size_t partition_id = 0; // Vertices marked as fixed cannot be modified by any local operation
 
            VectorNd displacement() const { return position - rest_position; }
 
            VectorNd prev_position(const int i) const
            {
                assert(0 <= i && i < projection_quantities.cols() / 3);
                return rest_position + projection_quantities.col(i);
            }
 
            VectorNd position_i(const int i) const
            {
                assert(i >= 0);
                if (i == 0)
                    return rest_position;
                else if (i - 1 < projection_quantities.cols() / 3)
                    return prev_position(i - 1);
                else
                    return position;
            }
 
            static VertexAttributes edge_collapse(
                const VertexAttributes &v0,
                const VertexAttributes &v1,
                const CollapseEdgeTo collapse_to);
        };
 
        struct EdgeAttributes
        {
            int op_depth = 0;
            int op_attempts = 0;
            // clang-format off
            enum class EnergyRank { BOTTOM, MIDDLE, TOP };
            // clang-format on
            EnergyRank energy_rank = EnergyRank::MIDDLE;
        };
 
        struct BoundaryAttributes : public EdgeAttributes
        {
            int boundary_id = -1;
        };
 
        struct ElementAttributes
        {
            int body_id = 0;
        };
 
        void write_edge_ranks_mesh(
            const EdgeMap<typename EdgeAttributes::EnergyRank> &elastic_ranks,
            const EdgeMap<typename EdgeAttributes::EnergyRank> &contact_ranks) const;
 
        EdgeAttributes &edge_attr(const size_t e_id);
 
        const EdgeAttributes &edge_attr(const size_t e_id) const;
 
        wmtk::AttributeCollection<VertexAttributes> vertex_attrs;
        wmtk::AttributeCollection<BoundaryAttributes> boundary_attrs;
        wmtk::AttributeCollection<ElementAttributes> element_attrs;
 
    protected:
        wmtk::ExecutePass<WildRemesher, EXECUTION_POLICY> executor;
        int m_n_quantities;
        double total_volume;
 
        // TODO: make this thread local
        typename std::conditional<
            std::is_same<WMTKMesh, wmtk::TriMesh>::value,
            std::shared_ptr<TriOperationCache>,
            std::shared_ptr<TetOperationCache>>::type op_cache;
 
    private:
        wmtk::AttributeCollection<EdgeAttributes> edge_attrs; // not used for tri mesh
    };
 
    using WildTriRemesher = WildRemesher<wmtk::TriMesh>;
    using WildTetRemesher = WildRemesher<wmtk::TetMesh>;
 
} // namespace polyfem::mesh