mesh.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2000-2022 Inria
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *
 *  * Redistributions of source code must retain the above copyright notice,
 *  this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright notice,
 *  this list of conditions and the following disclaimer in the documentation
 *  and/or other materials provided with the distribution.
 *  * Neither the name of the ALICE Project-Team nor the names of its
 *  contributors may be used to endorse or promote products derived from this
 *  software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 *  Contact: Bruno Levy
 *
 *     https://www.inria.fr/fr/bruno-levy
 *
 *     Inria,
 *     Domaine de Voluceau,
 *     78150 Le Chesnay - Rocquencourt
 *     FRANCE
 *
 */
 
#ifndef GEOGRAM_MESH_MESH
#define GEOGRAM_MESH_MESH
 
#include <geogram/basic/common.h>
#include <geogram/basic/range.h>
#include <geogram/basic/attributes.h>
#include <geogram/basic/vector_attribute.h>
#include <geogram/basic/geometry.h>
#include <tuple>
 
// GOMGEN (Graphite's pre-processor) does not understand
// modern C++ (but it is not a big drama since it is just
// used to generate the GUI).
#ifdef GOMGEN
#define MESH_NO_SYNTAXIC_SUGAR
#endif
 
namespace GEO {
 
    class Mesh;
 
    static constexpr index_t NO_VERTEX = NO_INDEX;
    static constexpr index_t NO_EDGE   = NO_INDEX;
    static constexpr index_t NO_FACET  = NO_INDEX;
    static constexpr index_t NO_CELL   = NO_INDEX;
    static constexpr index_t NO_CORNER = NO_INDEX;
 
    class GEOGRAM_API MeshSubElementsStore {
    public:
 
    MeshSubElementsStore(Mesh& mesh);
 
    virtual ~MeshSubElementsStore();
 
    index_t nb() const {
        return nb_;
    }
 
    AttributesManager& attributes() const {
        return const_cast<AttributesManager&>(attributes_);
    }
 
    index_as_iterator begin() const {
        return index_as_iterator(0);
    }
 
    index_as_iterator end() const {
        return index_as_iterator(nb());
    }
 
    protected:
 
    virtual void clear_store(
        bool keep_attributes, bool keep_memory = false
    );
 
    virtual void resize_store(index_t new_size);
 
 
    void reserve_store(index_t nb_to_reserve) {
        index_t nb = this->nb();
        resize_store(nb + nb_to_reserve);
        resize_store(nb);
    }
 
    index_t create_sub_elements(index_t nb) {
        index_t result = nb_;
        if(nb_ + nb > attributes_.size()) {
            index_t new_capacity=nb_ + nb;
            if(nb < 128) {
                new_capacity = std::max(index_t(16),attributes_.size());
                while(new_capacity < nb_ + nb) {
                    new_capacity *= 2;
                }
            }
            attributes_.reserve(new_capacity);
        }
        nb_ += nb;
        attributes_.resize(nb_);
        return result;
    }
 
    index_t create_sub_element() {
        index_t result = nb_;
        ++nb_;
        if(attributes_.capacity() < nb_) {
            index_t new_capacity =
                std::max(index_t(16),attributes_.capacity()*2);
            attributes_.reserve(new_capacity);
        }
        attributes_.resize(nb_);
        return result;
    }
 
    void adjust_store() {
        attributes_.resize(nb_);
    }
 
    void copy(
        const MeshSubElementsStore& rhs,
        bool copy_attributes = true
    ) {
        nb_ = rhs.nb();
        if(copy_attributes) {
            attributes_.copy(rhs.attributes_);
        } else {
            attributes_.clear(false,false);
            attributes_.resize(rhs.attributes_.size());
        }
    }
 
    protected:
    Mesh& mesh_;
    AttributesManager attributes_;
    index_t nb_;
    };
 
 
    /**************************************************************************/
 
    class GEOGRAM_API MeshElements {
    public:
    MeshElements();
    virtual ~MeshElements();
 
    virtual void delete_elements(
        vector<index_t>& to_delete,
        bool remove_isolated_vertices=true
    ) = 0;
 
    virtual void permute_elements(vector<index_t>& permutation) = 0;
 
    virtual void clear(
        bool keep_attributes=true, bool keep_memory=false
    ) = 0;
 
    virtual void pop() = 0;
 
    protected:
    static bool has_non_zero(const GEO::vector<index_t>& I) {
        for(index_t i = 0; i < I.size(); i++) {
            if(I[i] != 0) {
                return true;
            }
        }
        return false;
    }
    };
 
    /**************************************************************************/
 
    class MeshEdges;
    class MeshFacetCornersStore;
    class MeshCellCornersStore;
 
    class GEOGRAM_API MeshVertices :
        public MeshSubElementsStore, public MeshElements {
    public:
        MeshVertices(Mesh& mesh);
        ~MeshVertices() override;
 
        void remove_isolated();
 
        void delete_elements(
            vector<index_t>& to_delete, bool remove_isolated_vertices=true
        ) override;
 
        void permute_elements(vector<index_t>& permutation) override;
 
        index_t create_vertex() {
            return MeshSubElementsStore::create_sub_element();
        }
 
        index_t create_vertex(const double* coords) {
            // Sanity check:
            // It is not correct to call create_vertex(point_ptr(v)) since
            // create_vertex may realloc the points coordinates vector, thus
            // invalidate point_ptr(v).
            geo_debug_assert(
                nb() == 0 ||
                coords < point_ptr(0) ||
                coords >= point_ptr(0) + nb() * dimension()
            );
            index_t result = create_vertex();
            for(index_t c=0; c<dimension(); ++c) {
                point_ptr(result)[c] = coords[c];
            }
            return result;
        }
 
        template <index_t DIM> index_t create_vertex(
            const vecng<DIM,double>& p
        ) {
            geo_debug_assert(dimension() == DIM);
            return create_vertex(p.data());
        }
 
        index_t create_vertices(index_t nb) {
            return MeshSubElementsStore::create_sub_elements(nb);
        }
 
        void clear(
            bool keep_attributes=true, bool keep_memory=false
        ) override;
 
        void set_single_precision();
 
        void set_double_precision();
 
        bool single_precision() const {
            return point_fp32_.is_bound();
        }
 
        bool double_precision() const {
            return point_.is_bound();
        }
 
        index_t dimension() const {
            return
                single_precision() ?
                point_fp32_.dimension() :
                point_.dimension() ;
        }
 
        void set_dimension(index_t dim) {
            if(single_precision()) {
                point_fp32_.redim(dim);
            } else {
                point_.redim(dim);
            }
        }
 
        const double* point_ptr(index_t v) const {
            geo_debug_assert(v < nb());
            geo_debug_assert(!single_precision());
            return &point_[v*point_.dimension()];
        }
 
        double* point_ptr(index_t v) {
            geo_debug_assert(v < nb());
            geo_debug_assert(!single_precision());
            return &point_[v*point_.dimension()];
        }
 
 
        template<index_t DIM=3> vecng<DIM,double>& point(index_t v) {
            geo_debug_assert(v < nb());
            geo_debug_assert(!single_precision());
            geo_debug_assert(dimension() >= DIM);
            return *reinterpret_cast<vecng<DIM,double>*>(
                &point_[v*point_.dimension()]
            );
        }
 
        template <index_t DIM=3> const vecng<DIM,double>& point(
            index_t v
        ) const {
            geo_debug_assert(v < nb());
            geo_debug_assert(!single_precision());
            geo_debug_assert(dimension() >= DIM);
            return *reinterpret_cast<const vecng<DIM,double>*>(
                &point_[v*point_.dimension()]
            );
        }
 
        const float* single_precision_point_ptr(index_t v) const {
            geo_debug_assert(v < nb());
            geo_debug_assert(single_precision());
            return &point_fp32_[v*point_fp32_.dimension()];
        }
 
        float* single_precision_point_ptr(index_t v) {
            geo_debug_assert(v < nb());
            geo_debug_assert(single_precision());
            return &point_fp32_[v*point_fp32_.dimension()];
        }
 
        void assign_points(
            vector<double>& points, index_t dim, bool steal_arg
        );
 
        void assign_points(
            const double* points, index_t dim, index_t nb_pts
        );
 
        void pop() override;
 
#ifndef MESH_NO_SYNTAXIC_SUGAR
 
        template <index_t DIM = 3> auto points() const {
            typedef vecng<DIM,double> vecn;
            return transform_range_ref(
                index_range(0, nb()),
                [this](index_t v)->const vecn& {
                    // for MSVC that cannot chose among const/non-const versions
                    return *reinterpret_cast<const vecn*>(point_ptr(v));
                    // return point<DIM>(v); // MSVC does not understand this one
                }
            );
        }
 
        template <index_t DIM = 3> auto points() {
            typedef vecng<DIM,double> vecn;
            return transform_range_ref(
                index_range(0, nb()),
                [this](index_t v)->vecn& {
                    // for MSVC that cannot chose among const/non-const versions
                    return *reinterpret_cast<vecn*>(point_ptr(v));
                    // return point<DIM>(v); // MSVC does not understand this one
                }
            );
        }
 
#endif
 
    protected:
 
        void clear_store(
            bool keep_attributes, bool keep_memory = false
        ) override;
 
        void resize_store(index_t new_size) override;
 
        void bind_point_attribute(index_t dim, bool single_precision=false);
 
        void copy(const MeshVertices& rhs, bool copy_attributes=true) {
            index_t dim = rhs.dimension();
            if(point_fp32_.is_bound()) {
                point_fp32_.destroy();
            }
            if(point_.is_bound()) {
                point_.destroy();
            }
            MeshSubElementsStore::copy(rhs, copy_attributes);
            if(rhs.single_precision()) {
                point_fp32_.bind_if_is_defined(attributes(),"point_fp32");
                if(!point_fp32_.is_bound()) {
                    point_fp32_.create_vector_attribute(
                        attributes(), "point_fp32", dim
                    );
                }
            } else {
                point_.bind_if_is_defined(attributes(),"point");
                if(!point_.is_bound()) {
                    point_.create_vector_attribute(
                        attributes(), "point", dim
                    );
                }
            }
            // Even if we do not copy the attributes, we need at least
            // to copy the coordinates of the points !!
            if(!copy_attributes) {
                if(rhs.single_precision()) {
                    Memory::copy(
                        single_precision_point_ptr(0),
                        rhs.single_precision_point_ptr(0),
                        rhs.dimension()*rhs.nb()*sizeof(float)
                    );
                } else {
                    Memory::copy(
                        point_ptr(0),
                        rhs.point_ptr(0),
                        rhs.dimension()*rhs.nb()*sizeof(double)
                    );
                }
            }
        }
 
        MeshEdges& edges_;
        MeshFacetCornersStore& facet_corners_;
        MeshCellCornersStore& cell_corners_;
        Attribute<double> point_;
        Attribute<float> point_fp32_;
 
        friend class Mesh;
        friend class GeogramIOHandler;
    };
 
    /*************************************************************************/
 
    class GEOGRAM_API MeshEdges :
        public MeshSubElementsStore, public MeshElements {
    public:
        MeshEdges(Mesh& mesh);
        ~MeshEdges() override;
 
        index_t vertex(index_t e, index_t lv) const {
            geo_debug_assert(e < nb());
            geo_debug_assert(lv < 2);
            return edge_vertex_[2*e+lv];
        }
 
        void set_vertex(index_t e, index_t lv, index_t v) {
            geo_debug_assert(e < nb());
            geo_debug_assert(lv < 2);
            edge_vertex_[2*e+lv] = v;
        }
 
 
        index_t* vertex_index_ptr(index_t c) {
            geo_debug_assert(c < 2*nb());
            return &(edge_vertex_[c]);
        }
 
        const index_t* vertex_index_ptr(index_t c) const {
            geo_debug_assert(c < 2*nb());
            return &(edge_vertex_[c]);
        }
 
        index_t create_edge() {
            return create_sub_element();
        }
 
        index_t create_edges(index_t nb) {
            return create_sub_elements(nb);
        }
 
        index_t create_edge(index_t v1, index_t v2) {
            index_t result = create_edge();
            set_vertex(result,0,v1);
            set_vertex(result,1,v2);
            return result;
        }
 
        void delete_elements(
            vector<index_t>& to_delete, bool remove_isolated_vertices=true
        ) override;
 
        void permute_elements(vector<index_t>& permutation) override;
 
        void clear(
            bool keep_attributes=true, bool keep_memory=false
        ) override;
 
        void pop() override;
 
    protected:
        void clear_store(
            bool keep_attributes, bool keep_memory = false
        ) override;
 
        void resize_store(index_t new_size) override;
 
        index_t create_sub_element() {
            edge_vertex_.push_back(NO_VERTEX);
            edge_vertex_.push_back(NO_VERTEX);
            return MeshSubElementsStore::create_sub_element();
        }
 
        index_t create_sub_elements(index_t nb_in) {
            edge_vertex_.resize(2*(nb()+nb_in),NO_VERTEX);
            return MeshSubElementsStore::create_sub_elements(nb_in);
        }
 
        void copy(const MeshEdges& rhs, bool copy_attributes=true) {
            MeshSubElementsStore::copy(rhs, copy_attributes);
            edge_vertex_ = rhs.edge_vertex_;
        }
 
        vector<index_t> edge_vertex_;
        friend class Mesh;
        friend class GeogramIOHandler;
    };
 
    /**************************************************************************/
 
    class GEOGRAM_API MeshFacetsStore : public MeshSubElementsStore {
    public:
        MeshFacetsStore(Mesh& mesh);
 
        index_t corners_begin(index_t f) const {
            geo_debug_assert(f < nb());
            return (is_simplicial_ ? 3*f : facet_ptr_[f]);
        }
 
        index_t corners_end(index_t f) const {
            geo_debug_assert(f < nb());
            return (is_simplicial_ ? 3*(f+1): facet_ptr_[f+1]);
        }
 
        index_t nb_corners(index_t f) const {
            geo_debug_assert(f < nb());
            return (is_simplicial_ ? 3 : facet_ptr_[f+1] - facet_ptr_[f]);
        }
 
        index_t corner(index_t f, index_t lv) const {
            geo_debug_assert(f < nb());
            geo_debug_assert(lv < nb_corners(f));
            return corners_begin(f)+lv;
        }
 
        bool are_simplices() const {
            return is_simplicial_;
        }
 
        const index_t* corners_begin_ptr(index_t f) const {
            geo_debug_assert(!is_simplicial_);
            geo_debug_assert(f < nb());
            return &facet_ptr_[f];
        }
 
    protected:
        void clear_store(
            bool keep_attributes, bool keep_memory = false
        ) override;
 
        void resize_store(index_t new_size) override;
 
        index_t create_sub_element() {
            if(!is_simplicial_) {
                facet_ptr_.push_back(NO_CORNER);
            }
            return MeshSubElementsStore::create_sub_element();
        }
 
        index_t create_sub_elements(index_t nb) {
            if(!is_simplicial_) {
                for(index_t i=0; i<nb; ++i) {
                    facet_ptr_.push_back(NO_CORNER);
                }
            }
            return MeshSubElementsStore::create_sub_elements(nb);
        }
 
        void copy(const MeshFacetsStore& rhs, bool copy_attributes=true) {
            MeshSubElementsStore::copy(rhs,copy_attributes);
            is_simplicial_ = rhs.is_simplicial_;
            facet_ptr_ = rhs.facet_ptr_;
        }
 
    protected:
        bool is_simplicial_;
        vector<index_t> facet_ptr_;
        friend class Mesh;
        friend class GeogramIOHandler;
    };
 
    /*************************************************************************/
 
    class GEOGRAM_API MeshFacetCornersStore : public MeshSubElementsStore {
    public:
        MeshFacetCornersStore(Mesh& mesh);
 
        index_t vertex(index_t c) const {
            geo_assert(c < nb());
            return corner_vertex_[c];
        }
 
        index_t adjacent_facet(index_t c) const {
            geo_assert(c < nb());
            return corner_adjacent_facet_[c];
        }
 
        const index_t* adjacent_facet_ptr(index_t c) const {
            geo_assert(c < nb());
            return &corner_adjacent_facet_[c];
        }
 
 
        index_t* adjacent_facet_ptr(index_t c) {
            geo_assert(c < nb());
            return &corner_adjacent_facet_[c];
        }
 
        void set_vertex(index_t c, index_t v) {
            geo_debug_assert(c < nb());
            geo_debug_assert(v < vertices_.nb());
            corner_vertex_[c] = v;
        }
 
        void set_vertex_no_check(index_t c, index_t v) {
            geo_debug_assert(c < nb());
            corner_vertex_[c] = v;
        }
 
        void set_adjacent_facet(index_t c, index_t f) {
            geo_debug_assert(c < nb());
            geo_debug_assert(f == NO_FACET || f < facets_.nb());
            corner_adjacent_facet_[c] = f;
        }
 
        index_t* vertex_index_ptr(index_t c) {
            geo_debug_assert(c < nb());
            return &(corner_vertex_[c]);
        }
 
        const index_t* vertex_index_ptr(index_t c) const {
            geo_debug_assert(c < nb());
            return &(corner_vertex_[c]);
        }
 
        template <index_t DIM=3> vecng<DIM,double>& point(index_t c) {
            geo_debug_assert(c < nb());
            return vertices_.point<DIM>(vertex(c));
        }
 
        template <index_t DIM=3> const vecng<DIM,double>& point(
            index_t c
        ) const {
            geo_debug_assert(c < nb());
            return vertices_.point(vertex(c));
        }
 
    protected:
        void clear_store(
            bool keep_attributes, bool keep_memory = false
        ) override;
 
        void resize_store(index_t new_size) override;
 
        index_t create_sub_element(index_t v, index_t f = NO_FACET) {
            corner_vertex_.push_back(v);
            corner_adjacent_facet_.push_back(f);
            return MeshSubElementsStore::create_sub_element();
        }
 
        index_t create_sub_elements(index_t nb) {
            for(index_t i=0; i<nb; ++i) {
                corner_vertex_.push_back(NO_VERTEX);
            }
            for(index_t i=0; i<nb; ++i) {
                corner_adjacent_facet_.push_back(NO_FACET);
            }
            return MeshSubElementsStore::create_sub_elements(nb);
        }
 
        void copy(
            const MeshFacetCornersStore& rhs, bool copy_attributes=true
        ) {
            MeshSubElementsStore::copy(rhs, copy_attributes);
            corner_vertex_ = rhs.corner_vertex_;
            corner_adjacent_facet_ = rhs.corner_adjacent_facet_;
        }
 
    protected:
        MeshVertices& vertices_;
        MeshFacetsStore& facets_;
        vector<index_t> corner_vertex_;
        vector<index_t> corner_adjacent_facet_;
 
        friend class MeshFacets;
        friend class Mesh;
        friend class GeogramIOHandler;
    };
 
    /*************************************************************************/
 
    class GEOGRAM_API MeshFacets : public MeshFacetsStore, public MeshElements {
    public:
 
        MeshFacets(Mesh& mesh);
 
        index_t nb_vertices(index_t f) const {
            return nb_corners(f);
        }
 
        index_t vertex(index_t f, index_t lv) const {
            return facet_corners_.vertex(corner(f,lv));
        }
 
 
        template <index_t DIM=3> const vecng<DIM,double>& point(
            index_t f, index_t lv
        ) const {
            return vertices_.point<DIM>(vertex(f,lv));
        }
 
        template <index_t DIM=3> vecng<DIM,double>& point(
            index_t f, index_t lv
        ) {
            return vertices_.point<DIM>(vertex(f,lv));
        }
 
        void set_vertex(index_t f, index_t lv, index_t v) {
            facet_corners_.set_vertex(corner(f,lv),v);
        }
 
        index_t find_vertex(index_t f, index_t v) const {
            for(index_t lv=0; lv<nb_vertices(f); ++lv) {
                if(vertex(f,lv) == v) {
                    return lv;
                }
            }
            return NO_VERTEX;
        }
 
        index_t find_common_vertex(index_t f1, index_t f2) const {
            for(index_t lv=0; lv<nb_vertices(f1); ++lv) {
                index_t v = vertex(f1,lv);
                if(find_vertex(f2,v) != NO_VERTEX) {
                    return lv;
                }
            }
            return NO_VERTEX;
        }
 
        index_t adjacent(index_t f, index_t le) const {
            return facet_corners_.adjacent_facet(corner(f,le));
        }
 
        index_t find_adjacent(index_t f, index_t f2) const {
            for(index_t le=0; le<nb_vertices(f); ++le) {
                if(adjacent(f,le) == f2) {
                    return le;
                }
            }
            return NO_INDEX;
        }
 
        void set_adjacent(index_t f, index_t le, index_t f2) {
            facet_corners_.set_adjacent_facet(corner(f,le),f2);
        }
 
        index_t next_corner_around_facet(index_t f, index_t c) const {
            geo_debug_assert(f < nb());
            geo_debug_assert(c >= corners_begin(f) && c < corners_end(f));
            return c + 1 == corners_end(f) ? corners_begin(f) : c + 1;
        }
 
        index_t prev_corner_around_facet(index_t f, index_t c) const {
            geo_debug_assert(f < nb());
            geo_debug_assert(c >= corners_begin(f) && c < corners_end(f));
            return c == corners_begin(f) ? corners_end(f) - 1 : c - 1;
        }
 
        index_t find_edge(index_t f, index_t v1, index_t v2) const {
            for(index_t c1 = corners_begin(f); c1 != corners_end(f); ++c1) {
                index_t c2 = next_corner_around_facet(f,c1);
                if(
                    facet_corners_.vertex(c1) == v1 &&
                    facet_corners_.vertex(c2) == v2
                ) {
                    return c1 - corners_begin(f);
                }
            }
            return NO_INDEX;
        }
 
        void delete_elements(
            vector<index_t>& to_delete,
            bool remove_isolated_vertices=true
        ) override;
 
        void permute_elements(vector<index_t>& permutation) override;
 
        void clear(
            bool keep_attributes=true, bool keep_memory=false
        ) override;
 
        index_t create_facets(
            index_t nb_facets, index_t nb_vertices_per_polygon
        ) {
            if(nb_vertices_per_polygon != 3) {
                is_not_simplicial();
            }
 
            index_t first_facet = nb();
            index_t co = facet_corners_.nb();
            facet_corners_.create_sub_elements(
                nb_facets*nb_vertices_per_polygon
            );
            index_t result = create_sub_elements(nb_facets);
 
            if(!is_simplicial_) {
                for(index_t f=first_facet; f<=first_facet+nb_facets; ++f) {
                    facet_ptr_[f] = co;
                    co += nb_vertices_per_polygon;
                }
                geo_debug_assert(facet_ptr_.size() == nb()+1);
                geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
            }
            return result;
        }
 
        void reserve(index_t nb_to_reserve) {
            facet_corners_.reserve_store(nb_to_reserve*3);
            this->reserve_store(nb_to_reserve);
        }
 
        index_t create_triangles(index_t nb_triangles) {
            return create_facets(nb_triangles, 3);
        }
 
        index_t create_quads(index_t nb_quads) {
            return create_facets(nb_quads, 4);
        }
 
        index_t create_triangle(index_t v1, index_t v2, index_t v3) {
            geo_debug_assert(v1 != v2);
            geo_debug_assert(v2 != v3);
            geo_debug_assert(v3 != v1);
            facet_corners_.create_sub_element(v1);
            facet_corners_.create_sub_element(v2);
            facet_corners_.create_sub_element(v3);
            index_t result = create_sub_element();
            if(!is_simplicial_) {
                facet_ptr_[result+1] = facet_corners_.nb();
                geo_debug_assert(facet_ptr_.size() == nb()+1);
                geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
            }
            return result;
        }
 
        index_t create_quad(index_t v1, index_t v2, index_t v3, index_t v4) {
            is_not_simplicial();
            facet_corners_.create_sub_element(v1);
            facet_corners_.create_sub_element(v2);
            facet_corners_.create_sub_element(v3);
            facet_corners_.create_sub_element(v4);
            index_t result = create_sub_element();
            facet_ptr_[result+1] = facet_corners_.nb();
            geo_debug_assert(facet_ptr_.size() == nb()+1);
            geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
            return result;
        }
 
        index_t create_polygon(index_t nb_vertices) {
            if(nb_vertices != 3) {
                is_not_simplicial();
            }
            for(index_t i=0; i<nb_vertices; ++i) {
                facet_corners_.create_sub_element(NO_VERTEX);
            }
            index_t result = create_sub_element();
            if(!is_simplicial_) {
                facet_ptr_[result+1] = facet_corners_.nb();
                geo_debug_assert(facet_ptr_.size() == nb()+1);
                geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
            }
            return result;
        }
 
        index_t create_polygon(index_t nb_vertices, const index_t* vertices) {
            if(nb_vertices != 3) {
                is_not_simplicial();
            }
            for(index_t i=0; i<nb_vertices; ++i) {
                facet_corners_.create_sub_element(vertices[i]);
            }
            index_t result = create_sub_element();
            if(!is_simplicial_) {
                facet_ptr_[result+1] = facet_corners_.nb();
                geo_debug_assert(facet_ptr_.size() == nb()+1);
                geo_debug_assert(facet_ptr_[nb()] == facet_corners_.nb());
            }
            return result;
        }
 
        index_t create_polygon(const vector<index_t>& vertices) {
            return create_polygon(vertices.size(), vertices.data());
        }
 
        void connect();
 
        void connect(index_t facets_begin, index_t facets_end);
 
        void triangulate();
 
        void flip(index_t f);
 
        void compute_borders();
 
        void assign_triangle_mesh(
            coord_index_t dim,
            vector<double>& vertices,
            vector<index_t>& triangles,
            bool steal_args
        );
 
        /*
         * \brief Copies a triangle mesh into this Mesh.
         * \details Facet adjacence are not computed.
         *   Facet and corner attributes are zeroed.
         * \param[in] triangles facet to vertex links
         * \param[in] steal_args if set, vertices and triangles
         * are 'stolen' from the arguments
         * (using vector::swap).
         */
        void assign_triangle_mesh(
            vector<index_t>& triangles,
            bool steal_args
        );
 
        void pop() override;
 
        index_range corners(index_t f) const {
            geo_debug_assert(f < nb());
            return index_range(
                index_as_iterator(corners_begin(f)),
                index_as_iterator(corners_end(f))
            );
        }
 
#ifndef MESH_NO_SYNTAXIC_SUGAR
 
        auto vertices(index_t f) const {
            geo_debug_assert(f < nb());
            return transform_range(
                corners(f), [this](index_t c)->index_t {
                    return facet_corners_.vertex(c);
                }
            );
        }
 
        auto adjacent(index_t f) const {
            geo_debug_assert(f < nb());
            return transform_range(
                corners(f), [this](index_t c)->index_t {
                    return facet_corners_.adjacent_facet(c);
                }
            );
        }
 
        template <index_t DIM=3> auto points(index_t f) const {
            geo_debug_assert(f < nb());
            typedef vecng<DIM,double> vecn;
            return transform_range_ref(
                corners(f), [this](index_t c)->const vecn& {
                    index_t v = facet_corners_.vertex(c);
                    return vertices_.point<DIM>(v);
                }
            );
        }
 
        template <index_t DIM=3> auto points(index_t f) {
            geo_debug_assert(f < nb());
            typedef vecng<DIM,double> vecn;
            return transform_range_ref(
                corners(f), [this](index_t c)->vecn& {
                    index_t v = facet_corners_.vertex(c);
                    return vertices_.point<DIM>(v);
                }
            );
        }
 
        auto triangles(index_t f) const {
            geo_debug_assert(f < nb());
            index_t v0 = facet_corners_.vertex(corners_begin(f));
            return transform_range(
                index_range(
                    index_as_iterator(corners_begin(f)+1),
                    index_as_iterator(corners_end(f)-1)
                ),
                [this,v0](index_t c)->std::tuple<index_t, index_t, index_t> {
                    return std::make_tuple(
                        v0,
                        facet_corners_.vertex(c),
                        facet_corners_.vertex(c+1)
                    );
                }
            );
        }
 
        template <index_t DIM=3> auto triangle_points(index_t f) const {
            geo_debug_assert(f < nb());
            typedef vecng<DIM,double> vecn;
            return transform_range(
                triangles(f),
                [this](std::tuple<index_t, index_t, index_t> T)
                ->std::tuple<const vecn&, const vecn&, const vecn&> {
                    return std::tuple<const vecn&, const vecn&, const vecn&>(
                        vertices_.point<DIM>(std::get<0>(T)),
                        vertices_.point<DIM>(std::get<1>(T)),
                        vertices_.point<DIM>(std::get<2>(T))
                    );
                }
            );
        }
 
        template <index_t DIM=3> auto triangle_points(index_t f) {
            geo_debug_assert(f < nb());
            typedef vecng<DIM,double> vecn;
            return transform_range(
                triangles(f),
                [this](std::tuple<index_t, index_t, index_t> T)
                ->std::tuple<vecn&, vecn&, vecn&> {
                    return std::tuple<vecn&, vecn&, vecn&>(
                        vertices_.point(std::get<0>(T)),
                        vertices_.point(std::get<1>(T)),
                        vertices_.point(std::get<2>(T))
                    );
                }
            );
        }
 
#endif
 
    protected:
 
        void is_simplicial() {
            if(!is_simplicial_) {
                is_simplicial_ = true;
                facet_ptr_.resize(1);
                facet_ptr_[0] = 0;
            }
        }
 
        void is_not_simplicial() {
            if(is_simplicial_) {
                is_simplicial_ = false;
                facet_ptr_.resize(nb()+1);
                for(index_t f=0; f<facet_ptr_.size(); ++f) {
                    facet_ptr_[f] = 3*f;
                }
            }
        }
 
    protected:
        MeshVertices& vertices_;
        MeshFacetCornersStore& facet_corners_;
        friend class Mesh;
        friend class GeogramIOHandler;
        friend void GEOGRAM_API tessellate_facets(
            Mesh& M, index_t max_nb_vertices
        );
    };
 
    /*************************************************************************/
 
    enum MeshCellType {
        MESH_TET = 0,
        MESH_HEX = 1,
        MESH_PRISM = 2,
        MESH_PYRAMID = 3,
        MESH_CONNECTOR = 4,
        MESH_NB_CELL_TYPES = 5
    };
 
    struct CellDescriptor {
        index_t nb_vertices;
 
        index_t nb_facets;
 
        index_t nb_vertices_in_facet[6];
 
        index_t facet_vertex[6][4];
 
        index_t nb_edges;
 
        index_t edge_vertex[12][2];
 
        index_t edge_adjacent_facet[12][2];
    };
 
 
    namespace MeshCellDescriptors {
        GEOGRAM_API extern CellDescriptor*
        cell_type_to_cell_descriptor[GEO::MESH_NB_CELL_TYPES];
 
        GEOGRAM_API extern CellDescriptor tet_descriptor;
        GEOGRAM_API extern CellDescriptor hex_descriptor;
        GEOGRAM_API extern CellDescriptor prism_descriptor;
        GEOGRAM_API extern CellDescriptor pyramid_descriptor;
        GEOGRAM_API extern CellDescriptor connector_descriptor;
    }
 
    class GEOGRAM_API MeshCellsStore : public MeshSubElementsStore {
    public:
        MeshCellsStore(Mesh& mesh);
 
        bool are_simplices() const {
            return is_simplicial_;
        }
 
        MeshCellType type(index_t c) const {
            geo_debug_assert(c < nb());
            return is_simplicial_ ? MESH_TET : MeshCellType(cell_type_[c]);
        }
 
        const CellDescriptor& descriptor(index_t c) const;
 
        static const CellDescriptor& cell_type_to_cell_descriptor(
            MeshCellType t
        );
 
        index_t nb_corners(index_t c) const {
            geo_debug_assert(c < nb());
            return descriptor(c).nb_vertices;
        }
 
        index_t corners_begin(index_t c) const {
            geo_debug_assert(c < nb());
            return is_simplicial_ ? 4*c : cell_ptr_[c];
        }
 
        index_t corners_end(index_t c) const {
            geo_debug_assert(c < nb());
            return is_simplicial_ ? 4*(c+1) : cell_ptr_[c] + nb_corners(c);
        }
 
        index_t corner(index_t c, index_t lv) const {
            geo_debug_assert(c < nb());
            // There seems to be a linkage problem under MSVC for the
            // following assertion check...
#ifndef GEO_OS_WINDOWS
            geo_debug_assert(lv < nb_corners(c));
#endif
            return corners_begin(c) + lv;
        }
 
        index_t nb_facets(index_t c) const {
            geo_debug_assert(c < nb());
            return descriptor(c).nb_facets;
        }
 
        index_t facets_begin(index_t c) const {
            geo_debug_assert(c < nb());
            return is_simplicial_ ? 4*c : cell_ptr_[c];
        }
 
        index_t facets_end(index_t c) const {
            geo_debug_assert(c < nb());
            return is_simplicial_ ? 4*(c+1) : cell_ptr_[c] + nb_facets(c);
        }
 
        index_t facet(index_t c, index_t lf) const {
            geo_debug_assert(c < nb());
            geo_debug_assert(lf < nb_facets(c));
            return facets_begin(c) + lf;
        }
 
        index_t nb_edges(index_t c) const {
            return descriptor(c).nb_edges;
        }
 
    protected:
        void clear_store(
            bool keep_attributes, bool keep_memory = false
        ) override;
 
        void resize_store(index_t new_size) override;
 
        index_t create_sub_element(MeshCellType type) {
            if(!is_simplicial_) {
                cell_ptr_.push_back(NO_CORNER);
                cell_type_.push_back(Numeric::uint8(type));
            }
            return MeshSubElementsStore::create_sub_element();
        }
 
        index_t create_sub_elements(index_t nb, MeshCellType type) {
            if(!is_simplicial_) {
                for(index_t i=0; i<nb; ++i) {
                    cell_ptr_.push_back(NO_CORNER);
                    cell_type_.push_back(Numeric::uint8(type));
                }
            }
            return MeshSubElementsStore::create_sub_elements(nb);
        }
 
        void copy(
            const MeshCellsStore& rhs, bool copy_attributes=true
        ) {
            MeshSubElementsStore::copy(rhs, copy_attributes);
            is_simplicial_ = rhs.is_simplicial_;
            cell_type_ = rhs.cell_type_;
            cell_ptr_ = rhs.cell_ptr_;
        }
 
    protected:
        bool is_simplicial_;
        vector<Numeric::uint8> cell_type_;
        vector<index_t> cell_ptr_;
 
    protected:
        friend class Mesh;
        friend class GeogramIOHandler;
    };
 
    /*************************************************************************/
 
    class GEOGRAM_API MeshCellCornersStore : public MeshSubElementsStore {
    public:
        MeshCellCornersStore(Mesh& mesh);
 
        index_t vertex(index_t c) const {
            geo_assert(c < nb());
            return corner_vertex_[c];
        }
 
        void set_vertex(index_t c, index_t v) {
            geo_debug_assert(c < nb());
            geo_debug_assert(v < vertices_.nb());
            corner_vertex_[c] = v;
        }
 
        index_t* vertex_index_ptr(index_t c) {
            geo_debug_assert(c < nb());
            return &(corner_vertex_[c]);
        }
 
        const index_t* vertex_index_ptr(index_t c) const {
            geo_debug_assert(c < nb());
            return &(corner_vertex_[c]);
        }
 
        template <index_t DIM=3> vecng<DIM,double>& point(index_t c) {
            geo_debug_assert(c < nb());
            return vertices_.point<DIM>(vertex(c));
        }
 
        template <index_t DIM=3> const vecng<DIM,double>& point(
            index_t c
        ) const {
            geo_debug_assert(c < nb());
            return vertices_.point(vertex(c));
        }
 
    protected:
        void clear_store(
            bool keep_attributes, bool keep_memory = false
        ) override;
 
        void resize_store(index_t new_size) override;
 
        index_t create_sub_element(index_t v) {
            corner_vertex_.push_back(v);
            return MeshSubElementsStore::create_sub_element();
        }
 
        index_t create_sub_elements(index_t nb) {
            for(index_t i=0; i<nb; ++i) {
                corner_vertex_.push_back(NO_VERTEX);
            }
            return MeshSubElementsStore::create_sub_elements(nb);
        }
 
        void copy(
            const MeshCellCornersStore& rhs, bool copy_attributes=true
        ) {
            MeshSubElementsStore::copy(rhs, copy_attributes);
            corner_vertex_ = rhs.corner_vertex_;
        }
 
    protected:
        MeshVertices& vertices_;
        vector<index_t> corner_vertex_;
 
        friend class MeshCells;
        friend class Mesh;
        friend class GeogramIOHandler;
    };
 
    /*************************************************************************/
 
    class GEOGRAM_API MeshCellFacetsStore : public MeshSubElementsStore {
    public:
        MeshCellFacetsStore(Mesh& mesh);
 
        index_t adjacent_cell(index_t f) const {
            geo_assert(f < nb());
            return adjacent_cell_[f];
        }
 
        void set_adjacent_cell(index_t f, index_t c) {
            geo_debug_assert(f < nb());
            geo_debug_assert(c == NO_CELL || c < cells_.nb());
            adjacent_cell_[f] = c;
        }
 
        const index_t* adjacent_cell_ptr(index_t f) const {
            geo_assert(f < nb());
            return &adjacent_cell_[f];
        }
 
        index_t* adjacent_cell_ptr(index_t f) {
            geo_assert(f < nb());
            return &adjacent_cell_[f];
        }
 
    protected:
        void clear_store(
            bool keep_attributes, bool keep_memory = false
        ) override;
 
        void resize_store(index_t new_size) override;
 
        index_t create_sub_element(index_t c = NO_CELL) {
            adjacent_cell_.push_back(c);
            return MeshSubElementsStore::create_sub_element();
        }
 
        index_t create_sub_elements(index_t nb) {
            for(index_t i=0; i<nb; ++i) {
                adjacent_cell_.push_back(NO_CELL);
            }
            return MeshSubElementsStore::create_sub_elements(nb);
        }
 
        void copy(
            const MeshCellFacetsStore& rhs, bool copy_attributes=true
        ) {
            MeshSubElementsStore::copy(rhs, copy_attributes);
            adjacent_cell_ = rhs.adjacent_cell_;
        }
 
    protected:
        MeshVertices& vertices_;
        MeshCellsStore& cells_;
        vector<index_t> adjacent_cell_;
 
        friend class MeshCells;
        friend class Mesh;
        friend class GeogramIOHandler;
    };
 
    /*************************************************************************/
 
    class GEOGRAM_API MeshCells : public MeshCellsStore, public MeshElements {
    public:
        MeshCells(Mesh& mesh);
 
        index_t nb_vertices(index_t c) const {
            return nb_corners(c);
        }
 
        index_t vertex(index_t c, index_t lv) const {
            return cell_corners_.vertex(corner(c,lv));
        }
 
        void set_vertex(index_t c, index_t lv, index_t v) {
            cell_corners_.set_vertex(corner(c,lv),v);
        }
 
        template <index_t DIM=3> const vecng<DIM,double>& point(
            index_t c, index_t lv
        ) const {
            return vertices_.point<DIM>(vertex(c,lv));
        }
 
        template <index_t DIM=3> vecng<DIM,double>& point(
            index_t c, index_t lv
        ) {
            return vertices_.point<DIM>(vertex(c,lv));
        }
 
        index_t adjacent(index_t c, index_t lf) const {
            return cell_facets_.adjacent_cell(facet(c,lf));
        }
 
        void set_adjacent(index_t c, index_t lf, index_t c2) {
            cell_facets_.set_adjacent_cell(facet(c,lf),c2);
        }
 
        index_t facet_nb_vertices(index_t c, index_t lf) const {
            geo_debug_assert(lf < nb_facets(c));
            return descriptor(c).nb_vertices_in_facet[lf];
        }
 
        index_t facet_vertex(index_t c, index_t lf, index_t lv) const {
            geo_debug_assert(lv < facet_nb_vertices(c, lf));
            return cell_corners_.vertex(
                corner(c, descriptor(c).facet_vertex[lf][lv])
            );
        }
        index_t facet_corner(index_t c, index_t lf, index_t lc) const {
            geo_debug_assert(lc < facet_nb_vertices(c, lf));
            return corner(c, descriptor(c).facet_vertex[lf][lc]);
        }
 
        index_t edge_vertex(index_t c, index_t le, index_t lv) const {
            geo_debug_assert(le < nb_edges(c));
            geo_debug_assert(lv < 2);
            return cell_corners_.vertex(
                corner(c,descriptor(c).edge_vertex[le][lv])
            );
        }
 
        index_t edge_adjacent_facet(index_t c, index_t le, index_t lf) const {
            geo_debug_assert(le < nb_edges(c));
            geo_debug_assert(lf < 2);
            return descriptor(c).edge_adjacent_facet[le][lf];
        }
 
        index_range corners(index_t c) const {
            geo_debug_assert(c < nb());
            return index_range(
                index_as_iterator(corners_begin(c)),
                index_as_iterator(corners_end(c))
            );
        }
 
        index_range facets(index_t c) const {
            geo_debug_assert(c < nb());
            return index_range(
                index_as_iterator(facets_begin(c)),
                index_as_iterator(facets_end(c))
            );
        }
 
#ifndef MESH_NO_SYNTAXIC_SUGAR
 
        template <index_t DIM=3> auto points(index_t cell) const {
            geo_debug_assert(cell < nb());
            typedef vecng<DIM,double> vecn;
            return transform_range_ref(
                corners(cell), [this](index_t c)->const vecn& {
                    index_t v = cell_corners_.vertex(c);
                    return vertices_.point<DIM>(v);
                }
            );
        }
 
        template <index_t DIM=3> auto points(index_t cell) {
            geo_debug_assert(cell < nb());
            typedef vecng<DIM,double> vecn;
            return transform_range_ref(
                corners(cell), [this](index_t c)->vecn& {
                    index_t v = cell_corners_.vertex(c);
                    return vertices_.point<DIM>(v);
                }
            );
        }
 
        auto adjacent(index_t c) const {
            geo_debug_assert(c < nb());
            return transform_range(
                facets(c), [this](index_t f)->index_t {
                    return cell_facets_.adjacent_cell(f);
                }
            );
        }
 
#endif
 
        void clear(
            bool keep_attributes=true, bool keep_memory=false
        ) override;
 
        void delete_elements(
            vector<index_t>& to_delete,
            bool remove_isolated_vertices=true
        ) override;
 
        void permute_elements(vector<index_t>& permutation) override;
 
        index_t create_cells(index_t nb_cells, MeshCellType type) {
 
            if(nb_cells == 0) {
                return NO_CELL;
            }
 
 
            if(type != MESH_TET) {
                is_not_simplicial();
            }
 
            const CellDescriptor& desc = cell_type_to_cell_descriptor(type);
 
            //   Note: there is padding, the same number of corners and
            // faces is created for each cell, so that a single cell
            // pointer is used for both.
 
            index_t cell_size = std::max(desc.nb_vertices, desc.nb_facets);
            index_t first_cell = nb();
            index_t co = cell_corners_.nb();
 
            cell_corners_.create_sub_elements(
                nb_cells*cell_size
            );
 
            cell_facets_.create_sub_elements(
                nb_cells*cell_size
            );
 
            index_t result = create_sub_elements(nb_cells, type);
 
            if(!is_simplicial_) {
                for(index_t c=first_cell; c<=first_cell+nb_cells; ++c) {
                    cell_ptr_[c] = co;
                    co += cell_size;
                }
 
                geo_debug_assert(cell_ptr_.size() == nb()+1);
                geo_debug_assert(cell_ptr_[nb()] == cell_corners_.nb());
                geo_debug_assert(cell_ptr_[nb()] == cell_facets_.nb());
            }
 
            return result;
        }
 
        index_t create_tets(index_t nb_tets) {
            return create_cells(nb_tets, MESH_TET);
        }
 
        index_t create_hexes(index_t nb_hexes) {
            return create_cells(nb_hexes, MESH_HEX);
        }
 
        index_t create_prisms(index_t nb_prisms) {
            return create_cells(nb_prisms, MESH_PRISM);
        }
 
        index_t create_pyramids(index_t nb_pyramids) {
            return create_cells(nb_pyramids, MESH_PYRAMID);
        }
 
        index_t create_tet(
            index_t v1, index_t v2, index_t v3, index_t v4,
            index_t adj1 = NO_CELL,
            index_t adj2 = NO_CELL,
            index_t adj3 = NO_CELL,
            index_t adj4 = NO_CELL
        ) {
            cell_corners_.create_sub_element(v1);
            cell_corners_.create_sub_element(v2);
            cell_corners_.create_sub_element(v3);
            cell_corners_.create_sub_element(v4);
            cell_facets_.create_sub_element(adj1);
            cell_facets_.create_sub_element(adj2);
            cell_facets_.create_sub_element(adj3);
            cell_facets_.create_sub_element(adj4);
            index_t result = create_sub_element(MESH_TET);
            if(!is_simplicial_) {
                cell_ptr_[nb()] = cell_corners_.nb();
            }
            geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
            return result;
        }
 
        index_t create_hex(
            index_t v1, index_t v2, index_t v3, index_t v4,
            index_t v5, index_t v6, index_t v7, index_t v8,
            index_t adj1 = NO_CELL,
            index_t adj2 = NO_CELL,
            index_t adj3 = NO_CELL,
            index_t adj4 = NO_CELL,
            index_t adj5 = NO_CELL,
            index_t adj6 = NO_CELL
        ) {
            is_not_simplicial();
            cell_corners_.create_sub_element(v1);
            cell_corners_.create_sub_element(v2);
            cell_corners_.create_sub_element(v3);
            cell_corners_.create_sub_element(v4);
            cell_corners_.create_sub_element(v5);
            cell_corners_.create_sub_element(v6);
            cell_corners_.create_sub_element(v7);
            cell_corners_.create_sub_element(v8);
            cell_facets_.create_sub_element(adj1);
            cell_facets_.create_sub_element(adj2);
            cell_facets_.create_sub_element(adj3);
            cell_facets_.create_sub_element(adj4);
            cell_facets_.create_sub_element(adj5);
            cell_facets_.create_sub_element(adj6);
            cell_facets_.create_sub_element(NO_CELL); // padding
            cell_facets_.create_sub_element(NO_CELL); // padding
            index_t result = create_sub_element(MESH_HEX);
            cell_ptr_[nb()] = cell_corners_.nb();
            geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
            return result;
        }
 
        index_t create_prism(
            index_t v1, index_t v2,
            index_t v3, index_t v4,
            index_t v5, index_t v6,
            index_t adj1 = NO_CELL,
            index_t adj2 = NO_CELL,
            index_t adj3 = NO_CELL,
            index_t adj4 = NO_CELL,
            index_t adj5 = NO_CELL
        ) {
            is_not_simplicial();
            cell_corners_.create_sub_element(v1);
            cell_corners_.create_sub_element(v2);
            cell_corners_.create_sub_element(v3);
            cell_corners_.create_sub_element(v4);
            cell_corners_.create_sub_element(v5);
            cell_corners_.create_sub_element(v6);
            cell_facets_.create_sub_element(adj1);
            cell_facets_.create_sub_element(adj2);
            cell_facets_.create_sub_element(adj3);
            cell_facets_.create_sub_element(adj4);
            cell_facets_.create_sub_element(adj5);
            cell_facets_.create_sub_element(NO_CELL); // padding
            index_t result = create_sub_element(MESH_PRISM);
            cell_ptr_[nb()] = cell_corners_.nb();
            geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
            return result;
        }
 
        index_t create_pyramid(
            index_t v1, index_t v2, index_t v3, index_t v4, index_t v5,
            index_t adj1 = NO_CELL,
            index_t adj2 = NO_CELL,
            index_t adj3 = NO_CELL,
            index_t adj4 = NO_CELL,
            index_t adj5 = NO_CELL
        ) {
            is_not_simplicial();
            cell_corners_.create_sub_element(v1);
            cell_corners_.create_sub_element(v2);
            cell_corners_.create_sub_element(v3);
            cell_corners_.create_sub_element(v4);
            cell_corners_.create_sub_element(v5);
            cell_facets_.create_sub_element(adj1);
            cell_facets_.create_sub_element(adj2);
            cell_facets_.create_sub_element(adj3);
            cell_facets_.create_sub_element(adj4);
            cell_facets_.create_sub_element(adj5);
            index_t result = create_sub_element(MESH_PYRAMID);
            cell_ptr_[nb()] = cell_corners_.nb();
            geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
            return result;
        }
 
        index_t create_connector(
            index_t v1, index_t v2, index_t v3, index_t v4,
            index_t adj1 = NO_CELL,
            index_t adj2 = NO_CELL,
            index_t adj3 = NO_CELL
        ) {
            is_not_simplicial();
            cell_corners_.create_sub_element(v1);
            cell_corners_.create_sub_element(v2);
            cell_corners_.create_sub_element(v3);
            cell_corners_.create_sub_element(v4);
            cell_facets_.create_sub_element(adj1);
            cell_facets_.create_sub_element(adj2);
            cell_facets_.create_sub_element(adj3);
            cell_facets_.create_sub_element(NO_CELL); // padding
            index_t result = create_sub_element(MESH_CONNECTOR);
            cell_ptr_[nb()] = cell_corners_.nb();
            geo_debug_assert(cell_facets_.nb() == cell_corners_.nb());
            return result;
        }
 
        void connect(
            bool remove_trivial_slivers = true, bool verbose_if_OK=false
        );
 
        void compute_borders();
 
        void compute_borders(Attribute<index_t>& facet_cell);
 
        void assign_tet_mesh(
            coord_index_t dim,
            vector<double>& vertices,
            vector<index_t>& tets,
            bool steal_args
        );
 
        void assign_tet_mesh(
            vector<index_t>& tets,
            bool steal_args
        );
 
        void pop() override;
 
        index_t tet_adjacent(index_t t, index_t lf) const {
            geo_debug_assert(is_simplicial_);
            geo_debug_assert(t < nb());
            geo_debug_assert(lf < 4);
            return cell_facets_.adjacent_cell_[4*t+lf];
        }
 
        index_t find_tet_adjacent(index_t t, index_t t2) const {
            geo_debug_assert(is_simplicial_);
            geo_debug_assert(t < nb());
            geo_debug_assert(t2 < nb());
            for(index_t lf=0; lf<4; ++lf) {
                if(cell_facets_.adjacent_cell_[4*t+lf] == t2) {
                    return lf;
                }
            }
            return NO_FACET;
        }
 
        index_t tet_vertex(index_t t, index_t lv) const {
            geo_debug_assert(is_simplicial_);
            geo_debug_assert(t < nb());
            geo_debug_assert(lv < 4);
            return cell_corners_.corner_vertex_[4*t+lv];
        }
 
        index_t find_tet_vertex(index_t t, index_t v) const {
            geo_debug_assert(is_simplicial_);
            geo_debug_assert(t < nb());
            geo_debug_assert(v < vertices_.nb());
            for(index_t lv=0; lv<4; ++lv) {
                if(cell_corners_.corner_vertex_[4*t+lv] == v) {
                    return lv;
                }
            }
            return NO_VERTEX;
        }
 
        index_t tet_facet_vertex(
            index_t t, index_t lf, index_t lv
        ) const {
            geo_debug_assert(is_simplicial_);
            geo_debug_assert(t < nb());
            geo_debug_assert(lf < 4);
            geo_debug_assert(lv < 3);
            return cell_corners_.vertex(
                4 * t + local_tet_facet_vertex_index(lf,lv)
            );
        }
 
        index_t find_tet_facet(
            index_t t, index_t v1, index_t v2, index_t v3
        ) const {
            geo_debug_assert(is_simplicial_);
            for(index_t lf = 0; lf < 4; ++lf) {
                index_t w1 = tet_facet_vertex(t, lf, 0);
                index_t w2 = tet_facet_vertex(t, lf, 1);
                index_t w3 = tet_facet_vertex(t, lf, 2);
                if(
                    (v1 == w1 && v2 == w2 && v3 == w3) ||
                    (v1 == w2 && v2 == w3 && v3 == w1) ||
                    (v1 == w3 && v2 == w1 && v3 == w2)
                ) {
                    return lf;
                }
            }
            return NO_FACET;
        }
 
        static index_t local_tet_facet_vertex_index(index_t lf, index_t lv) {
            geo_debug_assert(lf < 4);
            geo_debug_assert(lv < 3);
            return MeshCellDescriptors::tet_descriptor.facet_vertex[lf][lv];
        }
 
    protected:
 
        void is_not_simplicial() {
            if(is_simplicial_) {
                is_simplicial_ = false;
                cell_ptr_.resize(nb()+1);
                cell_type_.assign(nb(), MESH_TET);
                for(index_t c=0; c<cell_ptr_.size(); ++c) {
                    cell_ptr_[c] = 4*c;
                }
            }
        }
 
        bool facets_match(
            index_t c1, index_t f1, index_t c2, index_t f2
        ) const;
 
        index_t find_cell_vertex(index_t c, index_t v) const {
            geo_debug_assert(c < nb());
            geo_debug_assert(v < vertices_.nb());
            for(index_t lv=0; lv<nb_vertices(c); ++lv) {
                if(vertex(c,lv) == v) {
                    return lv;
                }
            }
            return NO_VERTEX;
        }
 
        index_t find_cell_facet(
            index_t c1, index_t c2, index_t f2
        ) const {
            for(index_t f1=0; f1<nb_facets(c1); ++f1) {
                if(facets_match(c1,f1,c2,f2)) {
                    return f1;
                }
            }
            return NO_FACET;
        }
 
        bool triangular_facet_matches_quad_facet(
            index_t c1, index_t lf1,
            index_t c2, index_t lf2
        ) const;
 
 
        bool triangular_facets_have_common_edge(
            index_t c1, index_t f1,
            index_t c2, index_t f2,
            index_t& e1, index_t& e2
        ) const;
 
        bool create_connector(
            index_t c1, index_t lf1,
            const std::vector< std::pair<index_t, index_t> >& matches
        );
 
        void connect_tets();
 
    protected:
        MeshVertices& vertices_;
        MeshCellCornersStore& cell_corners_;
        MeshCellFacetsStore& cell_facets_;
        friend class Mesh;
        friend class GeogramIOHandler;
    };
 
    /*************************************************************************/
 
    enum MeshElementsFlags {
        MESH_NONE = 0,
        MESH_VERTICES = 1,
        MESH_FACETS = 2,
        MESH_EDGES  = 4,
        MESH_CELLS  = 8,
        MESH_ALL_ELEMENTS = 15,
        MESH_FACET_CORNERS = 16,
        MESH_CELL_CORNERS = 32,
        MESH_CELL_FACETS = 64,
        MESH_ALL_SUBELEMENTS = 65
    };
 
    /*************************************************************************/
 
    class GEOGRAM_API Mesh {
    public:
    MeshVertices          vertices;
    MeshEdges             edges;
    MeshFacets            facets;
    MeshFacetCornersStore facet_corners;
    MeshCells             cells;
    MeshCellCornersStore  cell_corners;
    MeshCellFacetsStore   cell_facets;
 
    Mesh(index_t dimension=3, bool single_precision=false);
 
    virtual ~Mesh();
 
    void clear(bool keep_attributes=true, bool keep_memory=false);
 
    void show_stats(const std::string& tag = "Mesh") const;
 
 
    void assert_is_valid();
 
 
    void copy(
        const Mesh& rhs,
        bool copy_attributes=true,
        MeshElementsFlags what=MESH_ALL_ELEMENTS
    );
 
 
    std::string get_attributes() const;
 
    std::string get_scalar_attributes() const;
 
    std::string get_vector_attributes(index_t max_dim = 0) const;
 
    index_t nb_subelements_types() const;
 
    MeshSubElementsStore& get_subelements_by_index(index_t i);
 
    const MeshSubElementsStore& get_subelements_by_index(index_t i) const;
 
 
    MeshSubElementsStore& get_subelements_by_type(MeshElementsFlags what);
 
    const MeshSubElementsStore& get_subelements_by_type(
        MeshElementsFlags what
    ) const;
 
    static std::string subelements_type_to_name(MeshElementsFlags what);
 
    static MeshElementsFlags name_to_subelements_type(
        const std::string& name
    );
 
    static bool parse_attribute_name(
        const std::string& full_attribute_name,
        MeshElementsFlags& where,
        std::string& attribute_name,
        index_t& component
    );
 
    protected:
    void display_attributes(
        const std::string& tag, const std::string& subelement_name,
        const MeshSubElementsStore& subelements
    ) const;
 
    Mesh(const Mesh& rhs) = delete;
 
    const Mesh& operator=(const Mesh& rhs) = delete;
    };
 
    /*************************************************************************/
}
 
#endif