Mesh.hpp

Go to the documentation of this file.

#pragma once
 
#include <polyfem/Common.hpp>
#include <polyfem/mesh/mesh2D/Navigation.hpp>
#include <polyfem/utils/Types.hpp>
#include <polyfem/utils/HashUtils.hpp>
#include <polyfem/utils/Types.hpp>
 
#include <Eigen/Dense>
#include <geogram/mesh/mesh.h>
 
#include <memory>
 
namespace polyfem
{
    namespace mesh
    {
        class Mesh;
 
        struct MeshWithID
        {
            int id;
            std::unique_ptr<Mesh> mesh;
        };
 
        enum class ElementType
        {
            SIMPLEX = 0,                   
            REGULAR_INTERIOR_CUBE,         
            SIMPLE_SINGULAR_INTERIOR_CUBE, 
            MULTI_SINGULAR_INTERIOR_CUBE,  
            REGULAR_BOUNDARY_CUBE,         
            SIMPLE_SINGULAR_BOUNDARY_CUBE, 
            MULTI_SINGULAR_BOUNDARY_CUBE,  
            INTERFACE_CUBE,                
            INTERIOR_POLYTOPE,             
            BOUNDARY_POLYTOPE,             
            PRISM,                         
            PYRAMID,                       
            UNDEFINED,                     
        };
 
        class Mesh
        {
        protected:
            class EdgeNodes
            {
            public:
                int v1, v2;
                std::vector<int> nodes_ids;
                Eigen::MatrixXd nodes;
            };
 
            class FaceNodes
            {
            public:
                int v1, v2, v3;
                std::vector<int> nodes_ids;
                Eigen::MatrixXd nodes;
            };
 
            class CellNodes
            {
            public:
                int v1, v2, v3, v4;
                std::vector<int> nodes_ids;
                Eigen::MatrixXd nodes;
            };
 
        public:
            static std::unique_ptr<Mesh> create(const std::string &path, const bool non_conforming = false);
 
            static std::unique_ptr<Mesh> create(GEO::Mesh &M, const bool non_conforming = false);
 
            static std::unique_ptr<Mesh> create(const Eigen::MatrixXd &vertices, const Eigen::MatrixXi &cells, const bool non_conforming = false);
 
            static std::unique_ptr<Mesh> create(const int dim, const bool non_conforming = false);
 
            virtual std::unique_ptr<Mesh> copy() const = 0;
 
            std::vector<MeshWithID> split() const;
 
        protected:
            Mesh() = default;
 
        public:
            virtual ~Mesh() = default;
            Mesh(Mesh &&) = default;
            Mesh &operator=(Mesh &&) = default;
            Mesh(const Mesh &) = default;
            Mesh &operator=(const Mesh &) = default;
 
            virtual void refine(const int n_refinement, const double t) = 0;
 
            virtual bool is_volume() const = 0;
            int dimension() const { return (is_volume() ? 3 : 2); }
            virtual bool is_conforming() const = 0;
            int n_elements() const { return (is_volume() ? n_cells() : n_faces()); }
            int n_boundary_elements() const { return (is_volume() ? n_faces() : n_edges()); }
 
            virtual int n_cells() const = 0;
            virtual int n_faces() const = 0;
            virtual int n_edges() const = 0;
            virtual int n_vertices() const = 0;
 
            virtual int n_face_vertices(const int f_id) const = 0;
            virtual int n_cell_vertices(const int c_id) const = 0;
            virtual int edge_vertex(const int e_id, const int lv_id) const = 0;
            virtual int face_vertex(const int f_id, const int lv_id) const = 0;
            virtual int cell_vertex(const int f_id, const int lv_id) const = 0;
            int element_vertex(const int el_id, const int lv_id) const
            {
                return (is_volume() ? cell_vertex(el_id, lv_id) : face_vertex(el_id, lv_id));
            }
 
            std::vector<int> element_vertices(const int el_id) const
            {
                std::vector<int> res;
                for (int i = 0; i < n_cell_vertices(el_id); ++i)
                    res.push_back(element_vertex(el_id, i));
                std::sort(res.begin(), res.end());
                return res;
            }
 
            int boundary_element_vertex(const int primitive_id, const int lv_id) const
            {
                return (is_volume() ? face_vertex(primitive_id, lv_id) : edge_vertex(primitive_id, lv_id));
            }
 
            virtual bool is_boundary_vertex(const int vertex_global_id) const = 0;
            virtual bool is_boundary_edge(const int edge_global_id) const = 0;
            virtual bool is_boundary_face(const int face_global_id) const = 0;
            virtual bool is_boundary_element(const int element_global_id) const = 0;
 
            virtual bool save(const std::string &path) const = 0;
 
        private:
            virtual bool build_from_matrices(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F) = 0;
 
        public:
            virtual void attach_higher_order_nodes(const Eigen::MatrixXd &V, const std::vector<std::vector<int>> &nodes) = 0;
            inline const Eigen::MatrixXi &orders() const { return orders_; }
            inline bool is_rational() const { return is_rational_; }
            inline void set_is_rational(const bool in_is_rational) { is_rational_ = in_is_rational; }
 
            virtual void normalize() = 0;
 
            virtual void compute_elements_tag() = 0;
            virtual void update_elements_tag() { assert(false); }
 
            virtual double edge_length(const int gid) const
            {
                assert(false);
                return 0;
            }
            virtual double quad_area(const int gid) const
            {
                assert(false);
                return 0;
            }
            virtual double tri_area(const int gid) const
            {
                assert(false);
                return 0;
            }
 
            virtual RowVectorNd point(const int global_index) const = 0;
            virtual void set_point(const int global_index, const RowVectorNd &p) = 0;
 
            virtual RowVectorNd edge_barycenter(const int e) const = 0;
            virtual RowVectorNd face_barycenter(const int f) const = 0;
            virtual RowVectorNd cell_barycenter(const int c) const = 0;
 
            void edge_barycenters(Eigen::MatrixXd &barycenters) const;
            void face_barycenters(Eigen::MatrixXd &barycenters) const;
            void cell_barycenters(Eigen::MatrixXd &barycenters) const;
            virtual void compute_element_barycenters(Eigen::MatrixXd &barycenters) const = 0;
 
            virtual void elements_boxes(std::vector<std::array<Eigen::Vector3d, 2>> &boxes) const = 0;
            virtual void barycentric_coords(const RowVectorNd &p, const int el_id, Eigen::MatrixXd &coord) const = 0;
 
            virtual void bounding_box(RowVectorNd &min, RowVectorNd &max) const = 0;
 
            bool is_spline_compatible(const int el_id) const;
            bool is_cube(const int el_id) const;
            bool is_polytope(const int el_id) const;
            bool is_simplex(const int el_id) const;
            bool is_prism(const int el_id) const;
 
            bool is_pyramid(const int el_id) const;
 
            const std::vector<ElementType> &elements_tag() const { return elements_tag_; }
            void set_tag(const int el, const ElementType type) { elements_tag_[el] = type; }
 
            void compute_node_ids(const std::function<int(const size_t, const RowVectorNd &, bool)> &marker);
 
            virtual void load_boundary_ids(const std::string &path);
 
            virtual void compute_boundary_ids(const std::function<int(const size_t, const std::vector<int> &, const RowVectorNd &, bool)> &marker) = 0;
 
            virtual void compute_body_ids(const std::function<int(const size_t, const std::vector<int> &, const RowVectorNd &)> &marker) = 0;
 
            virtual void set_boundary_ids(const std::vector<int> &boundary_ids)
            {
                assert(boundary_ids.size() == n_boundary_elements());
                boundary_ids_ = boundary_ids;
            }
            virtual void set_body_ids(const std::vector<int> &body_ids)
            {
                assert(body_ids.size() == n_elements());
                body_ids_ = body_ids;
            }
 
            virtual int get_default_boundary_id(const int primitive) const
            {
                if (is_volume() ? is_boundary_face(primitive) : is_boundary_edge(primitive))
                    return std::numeric_limits<int>::max(); // default for no selected boundary
                else
                    return -1; // default for no boundary
            }
 
            virtual int get_boundary_id(const int primitive) const
            {
                return has_boundary_ids() ? (boundary_ids_.at(primitive)) : get_default_boundary_id(primitive);
            }
 
            virtual int get_node_id(const int node_id) const
            {
                if (has_node_ids())
                    return node_ids_.at(node_id);
                else
                    return -1; // default for no boundary
            }
 
            void update_nodes(const Eigen::VectorXi &in_node_to_node);
 
            virtual int get_body_id(const int primitive) const
            {
                if (has_body_ids())
                    return body_ids_.at(primitive);
                else
                    return 0;
            }
            virtual const std::vector<int> &get_body_ids() const
            {
                return body_ids_;
            }
 
            void set_geometry_ids(const std::vector<int> &geometry_ids)
            {
                assert(geometry_ids.size() == n_elements());
                geometry_ids_ = geometry_ids;
            }
 
            int get_geometry_id(const int element) const
            {
                return has_geometry_ids() ? geometry_ids_.at(element) : 0;
            }
 
            const std::vector<int> &get_geometry_ids() const { return geometry_ids_; }
            bool has_geometry_ids() const { return !geometry_ids_.empty(); }
            bool has_node_ids() const { return !node_ids_.empty(); }
            bool has_boundary_ids() const { return !boundary_ids_.empty(); }
            virtual bool has_body_ids() const { return !body_ids_.empty(); }
 
            virtual void get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1) const = 0;
            virtual void get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1, const std::vector<bool> &valid_elements) const = 0;
 
            // /// @brief generate a triangular representation of every face
            // ///
            // /// @param[out] tris triangles connectivity
            // /// @param[out] pts triangles vertices
            // /// @param[out] ranges connection to original faces
            // virtual void triangulate_faces(Eigen::MatrixXi &tris, Eigen::MatrixXd &pts, std::vector<int> &ranges) const = 0;
 
            const std::vector<double> &cell_weights(const int cell_index) const { return cell_weights_[cell_index]; }
            void set_cell_weights(const std::vector<std::vector<double>> &in_cell_weights) { cell_weights_ = in_cell_weights; }
 
            virtual void prepare_mesh() {};
 
            bool has_poly() const
            {
                for (int i = 0; i < n_elements(); ++i)
                {
                    if (is_polytope(i))
                        return true;
                }
 
                return false;
            }
 
            bool has_prism() const
            {
                for (int i = 0; i < n_elements(); ++i)
                {
                    if (is_prism(i))
                        return true;
                }
 
                return false;
            }
 
            bool has_pyramids() const
            {
                for (int i = 0; i < n_elements(); ++i)
                {
                    if (is_pyramid(i))
                        return true;
                }
 
                return false;
            }
 
            bool is_simplicial() const
            {
                for (int i = 0; i < n_elements(); ++i)
                {
                    if (!is_simplex(i))
                        return false;
                }
 
                return true;
            }
 
            inline bool is_linear() const { return orders_.size() == 0 || orders_.maxCoeff() == 1; }
 
            std::vector<std::pair<int, int>> edges() const;
            std::vector<std::vector<int>> faces() const;
 
            std::unordered_map<std::pair<int, int>, size_t, polyfem::utils::HashPair> edges_to_ids() const;
            std::unordered_map<std::vector<int>, size_t, polyfem::utils::HashVector> faces_to_ids() const;
 
            inline const Eigen::VectorXi &in_ordered_vertices() const { return in_ordered_vertices_; }
            inline const Eigen::MatrixXi &in_ordered_edges() const { return in_ordered_edges_; }
            inline const Eigen::MatrixXi &in_ordered_faces() const { return in_ordered_faces_; }
 
            virtual void append(const Mesh &mesh);
 
            void append(const std::unique_ptr<Mesh> &mesh)
            {
                if (mesh != nullptr)
                    append(*mesh);
            }
 
            void apply_affine_transformation(const MatrixNd &A, const VectorNd &b);
 
        protected:
            virtual void remove_elements(const std::vector<bool> &keep) = 0;
            void filter_element_data(const std::vector<bool> &keep);
 
            virtual bool load(const std::string &path) = 0;
            virtual bool load(const GEO::Mesh &M) = 0;
 
            std::vector<ElementType> elements_tag_;
            std::vector<int> node_ids_;
            std::vector<int> boundary_ids_;
            std::vector<int> body_ids_;
            std::vector<int> geometry_ids_;
            Eigen::MatrixXi orders_;
            bool is_rational_ = false;
 
            std::vector<EdgeNodes> edge_nodes_;
            std::vector<FaceNodes> face_nodes_;
            std::vector<CellNodes> cell_nodes_;
            std::vector<std::vector<double>> cell_weights_;
 
            Eigen::VectorXi in_ordered_vertices_;
            Eigen::MatrixXi in_ordered_edges_;
            Eigen::MatrixXi in_ordered_faces_;
        };
    } // namespace mesh
} // namespace polyfem