geomodel_mesh_entity.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2017, Association Scientifique pour la Geologie et ses
 * Applications (ASGA). 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 ASGA 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 ASGA 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.
 *
 *     http://www.ring-team.org
 *
 *     RING Project
 *     Ecole Nationale Superieure de Geologie - GeoRessources
 *     2 Rue du Doyen Marcel Roubault - TSA 70605
 *     54518 VANDOEUVRE-LES-NANCY
 *     FRANCE
 */

#pragma once

#include <ringmesh/basic/common.h>

#include <ringmesh/geomodel/entity_type.h>
#include <ringmesh/geomodel/geomodel_entity.h>
#include <ringmesh/geomodel/geomodel_indexing_types.h>

#include <ringmesh/mesh/mesh.h>

namespace RINGMesh
{
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelGeologicalEntity );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelMeshEntityConstAccess );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelBuilderTopologyBase );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelBuilderTopology );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelBuilderGeometryBase );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelBuilderGeometry );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelBuilderGeology );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelBuilderRemovalBase );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelBuilderRemoval );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModelBuilderInfo );
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModel );

    ALIAS_2D_AND_3D( GeoModel );
} // namespace RINGMesh

namespace RINGMesh
{
    template < index_t DIMENSION >
    class RINGMESH_API GeoModelMeshEntity : public GeoModelEntity< DIMENSION >
    {
        ringmesh_disable_copy_and_move( GeoModelMeshEntity );
        ringmesh_template_assert_2d_or_3d( DIMENSION );
        friend class GeoModelMeshEntityAccess< DIMENSION >;
        friend class GeoModelMeshEntityConstAccess< DIMENSION >;

    public:
        virtual ~GeoModelMeshEntity();

        virtual MeshEntityType type_name() const = 0;

        gmme_id gmme() const
        {
            return gmme_id( type_name(), this->index() );
        }
        MeshEntityType mesh_entity_type() const
        {
            return gmme().type();
        }
        bool is_valid() const final
        {
            return is_mesh_valid() && are_geomodel_vertex_indices_valid();
        }
        virtual bool is_connectivity_valid() const;

        bool is_parent_connectivity_valid() const;

        index_t nb_boundaries() const
        {
            return static_cast< index_t >( boundaries_.size() );
        }
        const gmme_id& boundary_gmme( index_t x ) const;
        const GeoModelMeshEntity< DIMENSION >& boundary( index_t x ) const;

        index_t nb_incident_entities() const
        {
            return static_cast< index_t >( incident_entities_.size() );
        }
        const gmme_id& incident_entity_gmme( index_t x ) const;
        const GeoModelMeshEntity< DIMENSION >& incident_entity(
            index_t x ) const;

        bool has_inside_border() const;
        bool is_inside_border( const GeoModelMeshEntity& rhs ) const;

        bool has_parent() const
        {
            return nb_parents() != 0;
        }

        bool has_parent( const GeologicalEntityType& parent_type ) const
        {
            return could_be_undefined_parent_gmge( parent_type ).is_defined();
        }

        index_t nb_parents() const
        {
            return static_cast< index_t >( parents_.size() );
        }

        const gmge_id& parent_gmge( index_t id ) const;

        const gmge_id parent_gmge(
            const GeologicalEntityType& parent_type ) const;
        const GeoModelGeologicalEntity< DIMENSION >& parent(
            index_t parent_index ) const;
        const GeoModelGeologicalEntity< DIMENSION >& parent(
            const GeologicalEntityType& parent_type ) const;

        void save( const std::string& filename ) const
        {
            mesh_->save_mesh( filename );
        }
        const NNSearch< DIMENSION >& vertex_nn_search() const
        {
            return mesh_->vertex_nn_search();
        }

        index_t nb_vertices() const
        {
            return mesh_->nb_vertices();
        }
        const vecn< DIMENSION >& vertex( index_t vertex_index ) const
        {
            return mesh_->vertex( vertex_index );
        }

        virtual index_t nb_mesh_elements() const = 0;
        virtual index_t nb_mesh_element_vertices(
            index_t mesh_element_index ) const = 0;
        virtual index_t mesh_element_vertex_index(
            const ElementLocalVertex& element_local_vertex ) const = 0;
        const vecn< DIMENSION >& mesh_element_vertex(
            const ElementLocalVertex& element_local_vertex ) const
        {
            return vertex( mesh_element_vertex_index( element_local_vertex ) );
        }

        virtual double mesh_element_size(
            index_t mesh_element_index ) const = 0;
        virtual vecn< DIMENSION > mesh_element_barycenter(
            index_t mesh_element_index ) const = 0;
        vecn< DIMENSION > entity_barycenter() const
        {
            vecn< DIMENSION > result;
            for( auto v : range( nb_vertices() ) )
            {
                result += vertex( v );
            }
            ringmesh_assert( nb_vertices() > 0 );
            return result / static_cast< double >( nb_vertices() );
        }

        virtual double size() const
        {
            double size = 0.0;
            for( auto i : range( nb_mesh_elements() ) )
            {
                size += mesh_element_size( i );
            }
            return size;
        }

        GEO::AttributesManager& vertex_attribute_manager() const
        {
            return mesh_->vertex_attribute_manager();
        }

    protected:
        GeoModelMeshEntity( const GeoModel< DIMENSION >& geomodel, index_t id )
            : GeoModelEntity< DIMENSION >( geomodel, id )
        {
        }

        virtual void copy_mesh_entity(
            const GeoModelMeshEntity< DIMENSION >& from )
        {
            this->copy_name( from );
            this->id_ = from.index();
            boundaries_ = from.boundaries_;
            incident_entities_ = from.incident_entities_;
            parents_ = from.parents_;
        }

        virtual bool is_index_valid() const final;
        virtual bool is_mesh_valid() const
        {
            return mesh_ != nullptr;
        }

        void set_mesh( std::shared_ptr< MeshBase< DIMENSION > > mesh )
        {
            ringmesh_assert( mesh != nullptr );
            mesh_ = std::move( mesh );
        }

        bool is_boundary_connectivity_valid() const;
        bool is_incident_entity_connectivity_valid() const;
        bool are_geomodel_vertex_indices_valid() const;

        void unbind_vertex_mapping_attribute() const;
        void bind_vertex_mapping_attribute() const;

        virtual void change_mesh_data_structure( const MeshType& type ) = 0;

    private:
        gmge_id defined_parent_gmge(
            const GeologicalEntityType& parent_type ) const;

        gmge_id could_be_undefined_parent_gmge(
            const GeologicalEntityType& parent_type ) const;

    protected:
        std::vector< index_t > boundaries_{};

        std::vector< index_t > incident_entities_{};

        std::vector< index_t > parents_{};

    private:
        std::shared_ptr< MeshBase< DIMENSION > > mesh_{};
    };
    ALIAS_2D_AND_3D( GeoModelMeshEntity );

    template < index_t DIMENSION >
    class RINGMESH_API Corner final : public GeoModelMeshEntity< DIMENSION >
    {
        ringmesh_disable_copy_and_move( Corner );
        ringmesh_template_assert_2d_or_3d( DIMENSION );

    public:
        friend class GeoModelMeshEntityAccess< DIMENSION >;
        friend class GeoModelMeshEntityConstAccess< DIMENSION >;

        virtual ~Corner()
        {
            this->unbind_vertex_mapping_attribute();
        }

        static MeshEntityType type_name_static()
        {
            return MeshEntityType( "Corner" );
        }

        MeshEntityType type_name() const final
        {
            return type_name_static();
        }

        bool is_on_voi() const final;

        index_t nb_mesh_elements() const final
        {
            return 0;
        }

        index_t nb_mesh_element_vertices( index_t mesh_element = 0 ) const final
        {
            ringmesh_unused( mesh_element );
            index_t nb_vertices = point_set_mesh_->nb_vertices();
            ringmesh_assert( nb_vertices < 2 );
            return nb_vertices;
        }

        const Line< DIMENSION >& incident_entity( index_t x ) const;

        double mesh_element_size( index_t mesh_element = 0 ) const final
        {
            ringmesh_unused( mesh_element );
            return 0.0;
        }

        double size() const final
        {
            return 0.0;
        }

        vecn< DIMENSION > mesh_element_barycenter(
            index_t mesh_element = 0 ) const final
        {
            ringmesh_unused( mesh_element );
            return this->vertex( 0 );
        }

        const PointSetMesh< DIMENSION >& mesh() const
        {
            return *point_set_mesh_;
        }

    protected:
        Corner( const GeoModel< DIMENSION >& geomodel,
            index_t id,
            const MeshType& type )
            : GeoModelMeshEntity< DIMENSION >( geomodel, id )
        {
            update_mesh_storage_type(
                PointSetMesh< DIMENSION >::create_mesh( type ) );
        }

        index_t mesh_element_vertex_index(
            const ElementLocalVertex& element_local_vertex = ElementLocalVertex(
                0, 0 ) ) const final
        {
            ringmesh_unused( element_local_vertex );
            return 0;
        }

        bool is_mesh_valid() const final;

    private:
        void update_mesh_storage_type(
            std::unique_ptr< PointSetMesh< DIMENSION > > mesh )
        {
            point_set_mesh_ = std::move( mesh );
            GeoModelMeshEntity< DIMENSION >::set_mesh( point_set_mesh_ );
        }

        void change_mesh_data_structure( const MeshType& type ) final;

    private:
        std::shared_ptr< PointSetMesh< DIMENSION > > point_set_mesh_{};
    };
    ALIAS_2D_AND_3D( Corner );

    template < index_t DIMENSION >
    class RINGMESH_API Line final : public GeoModelMeshEntity< DIMENSION >
    {
        ringmesh_disable_copy_and_move( Line );
        ringmesh_template_assert_2d_or_3d( DIMENSION );

    public:
        friend class GeoModelMeshEntityAccess< DIMENSION >;

        virtual ~Line()
        {
            this->unbind_vertex_mapping_attribute();
        }

        static MeshEntityType type_name_static()
        {
            return MeshEntityType( "Line" );
        }

        MeshEntityType type_name() const final
        {
            return type_name_static();
        }

        bool is_on_voi() const final;

        const Corner< DIMENSION >& boundary( index_t x ) const;

        const Surface< DIMENSION >& incident_entity( index_t x ) const
        {
            return static_cast< const Surface< DIMENSION >& >(
                GeoModelMeshEntity< DIMENSION >::incident_entity( x ) );
        }

        bool is_connectivity_valid() const final;

        const LineAABBTree< DIMENSION >& edge_aabb() const
        {
            return line_mesh_->edge_aabb();
        }

        const NNSearch< DIMENSION >& edge_nn_search() const
        {
            return line_mesh_->edge_nn_search();
        }

        index_t nb_mesh_elements() const final
        {
            return line_mesh_->nb_edges();
        }

        index_t nb_mesh_element_vertices( index_t mesh_element = 0 ) const final
        {
            ringmesh_unused( mesh_element );
            return 2;
        }

        index_t mesh_element_vertex_index(
            const ElementLocalVertex& element_local_vertex ) const final
        {
            ringmesh_assert(
                element_local_vertex.element_id_ < nb_mesh_elements() );
            ringmesh_assert( element_local_vertex.local_vertex_id_ < 2 );
            return line_mesh_->edge_vertex( element_local_vertex );
        }

        bool is_closed() const
        {
            ringmesh_assert( this->nb_boundaries() == 2 );
            return ( this->boundary_gmme( 0 ).is_defined() )
                   && ( this->boundary_gmme( 0 ) == this->boundary_gmme( 1 ) );
        }

        double mesh_element_size( index_t edge_index ) const final
        {
            ringmesh_assert( edge_index < nb_mesh_elements() );
            return line_mesh_->edge_length( edge_index );
        }

        vecn< DIMENSION > mesh_element_barycenter(
            index_t edge_index ) const final
        {
            ringmesh_assert( edge_index < nb_mesh_elements() );
            return line_mesh_->edge_barycenter( edge_index );
        }

        bool is_first_corner_first_vertex() const;

        const LineMesh< DIMENSION >& mesh() const
        {
            return *line_mesh_;
        }

    protected:
        Line( const GeoModel< DIMENSION >& geomodel,
            index_t id,
            const MeshType& type )
            : GeoModelMeshEntity< DIMENSION >( geomodel, id )
        {
            update_mesh_storage_type(
                LineMesh< DIMENSION >::create_mesh( type ) );
        }

        bool is_mesh_valid() const final;

    private:
        void update_mesh_storage_type(
            std::unique_ptr< LineMesh< DIMENSION > > mesh )
        {
            line_mesh_ = std::move( mesh );
            GeoModelMeshEntity< DIMENSION >::set_mesh( line_mesh_ );
        }

        void change_mesh_data_structure( const MeshType& type ) final;

    private:
        std::shared_ptr< LineMesh< DIMENSION > > line_mesh_{};
    };
    ALIAS_2D_AND_3D( Line );

    template < index_t DIMENSION >
    class RINGMESH_API SurfaceBase : public GeoModelMeshEntity< DIMENSION >
    {
        ringmesh_disable_copy_and_move( SurfaceBase );
        ringmesh_template_assert_2d_or_3d( DIMENSION );

    public:
        friend class GeoModelMeshEntityAccess< DIMENSION >;

        virtual ~SurfaceBase()
        {
            this->unbind_vertex_mapping_attribute();
        }

        MeshEntityType type_name() const final
        {
            return type_name_static();
        }

        static MeshEntityType type_name_static()
        {
            return MeshEntityType( "Surface" );
        }

        const Line< DIMENSION >& boundary( index_t x ) const;

        const Region< DIMENSION >& incident_entity( index_t x ) const;

        bool is_simplicial() const
        {
            return surface_mesh_->polygons_are_simplicies();
        }

        const SurfaceAABBTree< DIMENSION >& polygon_aabb() const
        {
            return surface_mesh_->polygon_aabb();
        }

        const NNSearch< DIMENSION >& polygon_nn_search() const
        {
            return surface_mesh_->polygon_nn_search();
        }

        GEO::AttributesManager& polygon_attribute_manager() const
        {
            return surface_mesh_->polygon_attribute_manager();
        }

        index_t nb_mesh_elements() const final
        {
            return surface_mesh_->nb_polygons();
        }

        index_t nb_mesh_element_vertices( index_t polygon_index ) const final
        {
            ringmesh_assert( polygon_index < nb_mesh_elements() );
            return surface_mesh_->nb_polygon_vertices( polygon_index );
        }

        index_t mesh_element_vertex_index(
            const ElementLocalVertex& element_local_vertex ) const final
        {
            ringmesh_assert(
                element_local_vertex.element_id_ < nb_mesh_elements() );
            ringmesh_assert( element_local_vertex.local_vertex_id_
                             < nb_mesh_element_vertices(
                                   element_local_vertex.element_id_ ) );
            return surface_mesh_->polygon_vertex( element_local_vertex );
        }

        index_t polygon_adjacent_index(
            const PolygonLocalEdge& polygon_local_edge ) const
        {
            ringmesh_assert(
                polygon_local_edge.polygon_id_ < nb_mesh_elements() );
            ringmesh_assert(
                polygon_local_edge.local_edge_id_
                < nb_mesh_element_vertices( polygon_local_edge.polygon_id_ ) );
            return surface_mesh_->polygon_adjacent( polygon_local_edge );
        }

        vecn< DIMENSION > mesh_element_barycenter(
            index_t polygon_index ) const final
        {
            ringmesh_assert( polygon_index < nb_mesh_elements() );
            return surface_mesh_->polygon_barycenter( polygon_index );
        }

        double mesh_element_size( index_t polygon_index ) const final
        {
            ringmesh_assert( polygon_index < nb_mesh_elements() );
            return surface_mesh_->polygon_area( polygon_index );
        }

        const SurfaceMesh< DIMENSION >& mesh() const
        {
            return *surface_mesh_;
        }

    protected:
        SurfaceBase( const GeoModel< DIMENSION >& geomodel,
            index_t id,
            const MeshType type )
            : GeoModelMeshEntity< DIMENSION >( geomodel, id )
        {
            update_mesh_storage_type(
                SurfaceMesh< DIMENSION >::create_mesh( type ) );
        }

        bool is_mesh_valid() const final;

    private:
        void update_mesh_storage_type(
            std::unique_ptr< SurfaceMesh< DIMENSION > > mesh )
        {
            surface_mesh_ = std::move( mesh );
            GeoModelMeshEntity< DIMENSION >::set_mesh( surface_mesh_ );
        }

        void change_mesh_data_structure( const MeshType& type ) final;

    private:
        std::shared_ptr< SurfaceMesh< DIMENSION > > surface_mesh_{};
    };

    template < index_t DIMENSION >
    class RINGMESH_API Surface final : public SurfaceBase< DIMENSION >
    {
    };

    template <>
    class RINGMESH_API Surface< 2 > final : public SurfaceBase< 2 >
    {
        friend class GeoModelMeshEntityAccess< 2 >;

    private:
        Surface( const GeoModel2D& geomodel, index_t id, const MeshType type )
            : SurfaceBase< 2 >( geomodel, id, type )
        {
        }

    public:
        bool is_on_voi() const final;
        bool side( index_t i ) const
        {
            return sides_[i];
        }

        bool is_meshed() const
        {
            return mesh().nb_polygons() > 0;
        }

    private:
        std::vector< bool > sides_{};
    };

    template <>
    class RINGMESH_API Surface< 3 > final : public SurfaceBase< 3 >
    {
        friend class GeoModelMeshEntityAccess< 3 >;

    private:
        Surface( const GeoModel3D& geomodel, index_t id, const MeshType type )
            : SurfaceBase< 3 >( geomodel, id, type )
        {
        }

    public:
        bool is_on_voi() const final;
        const Region< 3 >& incident_entity( index_t x ) const;
    };
    ALIAS_2D_AND_3D( Surface );

    template < index_t DIMENSION >
    class RINGMESH_API Region final : public GeoModelMeshEntity< DIMENSION >
    {
        ringmesh_disable_copy_and_move( Region );
        ringmesh_template_assert_3d( DIMENSION );

    public:
        friend class GeoModelMeshEntityAccess< DIMENSION >;

        virtual ~Region()
        {
            this->unbind_vertex_mapping_attribute();
        }

        static MeshEntityType type_name_static()
        {
            return MeshEntityType( "Region" );
        }

        MeshEntityType type_name() const final
        {
            return type_name_static();
        }

        bool is_on_voi() const final;

        const Surface< DIMENSION >& boundary( index_t x ) const;

        bool is_connectivity_valid() const final;

        bool is_meshed() const
        {
            return volume_mesh_->nb_cells() > 0;
        }

        bool is_simplicial() const
        {
            return volume_mesh_->cells_are_simplicies();
        }

        const VolumeAABBTree< DIMENSION >& cell_aabb() const
        {
            return volume_mesh_->cell_aabb();
        }

        const NNSearch< DIMENSION >& cell_nn_search() const
        {
            return volume_mesh_->cell_nn_search();
        }

        GEO::AttributesManager& cell_attribute_manager() const
        {
            return volume_mesh_->cell_attribute_manager();
        }

        index_t nb_mesh_elements() const final
        {
            return volume_mesh_->nb_cells();
        }

        index_t nb_mesh_element_vertices( index_t cell_index ) const final
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                return volume_mesh_->nb_cell_vertices( cell_index );
            }
            ringmesh_assert_not_reached;
            return NO_ID;
        }

        index_t mesh_element_vertex_index(
            const ElementLocalVertex& element_local_vertex ) const final
        {
            if( is_meshed() )
            {
                ringmesh_assert(
                    element_local_vertex.element_id_ < nb_mesh_elements() );
                ringmesh_assert( element_local_vertex.local_vertex_id_
                                 < nb_mesh_element_vertices(
                                       element_local_vertex.element_id_ ) );
                return volume_mesh_->cell_vertex( element_local_vertex );
            }
            ringmesh_assert_not_reached;
            return NO_ID;
        }

        CellType cell_type( index_t cell_index ) const
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                return volume_mesh_->cell_type( cell_index );
            }
            ringmesh_assert_not_reached;
            return CellType::UNDEFINED;
        }

        index_t nb_cell_edges( index_t cell_index ) const
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                return volume_mesh_->nb_cell_edges( cell_index );
            }
            ringmesh_assert_not_reached;
            return NO_ID;
        }

        index_t nb_cell_facets( index_t cell_index ) const
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                return volume_mesh_->nb_cell_facets( cell_index );
            }
            ringmesh_assert_not_reached;
            return NO_ID;
        }

        index_t nb_cell_facet_vertices(
            index_t cell_index, index_t facet_index ) const
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                ringmesh_assert( facet_index < nb_cell_facets( cell_index ) );
                return volume_mesh_->nb_cell_facet_vertices(
                    CellLocalFacet( cell_index, facet_index ) );
            }
            ringmesh_assert_not_reached;
            return NO_ID;
        }

        index_t cell_edge_vertex_index(
            index_t cell_index, index_t edge_index, index_t vertex_index ) const
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                ringmesh_assert( edge_index < nb_cell_edges( cell_index ) );
                ringmesh_assert(
                    vertex_index < nb_mesh_element_vertices( cell_index ) );
                return volume_mesh_->cell_edge_vertex(
                    cell_index, edge_index, vertex_index );
            }
            ringmesh_assert_not_reached;
            return NO_ID;
        }

        index_t cell_facet_vertex_index( index_t cell_index,
            index_t facet_index,
            index_t vertex_index ) const
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                ringmesh_assert( facet_index < nb_cell_facets( cell_index ) );
                ringmesh_assert(
                    vertex_index < nb_mesh_element_vertices( cell_index ) );
                return volume_mesh_->cell_facet_vertex(
                    CellLocalFacet( cell_index, facet_index ), vertex_index );
            }
            ringmesh_assert_not_reached;
            return NO_ID;
        }

        index_t cell_adjacent_index(
            index_t cell_index, index_t facet_index ) const
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                ringmesh_assert( facet_index < nb_cell_facets( cell_index ) );
                return volume_mesh_->cell_adjacent(
                    CellLocalFacet( cell_index, facet_index ) );
            }
            ringmesh_assert_not_reached;
            return NO_ID;
        }

        ElementLocalVertex find_cell_from_colocated_vertex_if_any(
            const vecn< DIMENSION >& vertex_vec ) const;

        double mesh_element_size( index_t cell_index ) const final
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                return volume_mesh_->cell_volume( cell_index );
            }
            ringmesh_assert_not_reached;
            return 0;
        }
        double size() const final
        {
            double result = 0.;
            for( auto i : range( this->nb_boundaries() ) )
            {
                const Surface< DIMENSION >& surface = boundary( i );
                for( auto t : range( surface.nb_mesh_elements() ) )
                {
                    const vecn< DIMENSION >& p0 = surface.mesh_element_vertex(
                        ElementLocalVertex( t, 0 ) );
                    for( auto v :
                        range( 1, surface.nb_mesh_element_vertices( t ) - 1 ) )
                    {
                        double cur_volume =
                            ( dot( p0,
                                cross( surface.mesh_element_vertex(
                                           ElementLocalVertex( t, v ) ),
                                    surface.mesh_element_vertex(
                                        ElementLocalVertex( t, v + 1 ) ) ) ) )
                            / 6.;
                        side( i ) ? result -= cur_volume : result += cur_volume;
                    }
                }
            }
            return std::fabs( result );
        }
        vecn< DIMENSION > mesh_element_barycenter(
            index_t cell_index ) const final
        {
            if( is_meshed() )
            {
                ringmesh_assert( cell_index < nb_mesh_elements() );
                return volume_mesh_->cell_barycenter( cell_index );
            }
            ringmesh_assert_not_reached;
            return vecn< DIMENSION >();
        }

        std::vector< index_t > cells_around_vertex(
            index_t vertex_id, index_t cell_hint ) const
        {
            return volume_mesh_->cells_around_vertex( vertex_id, cell_hint );
        }

        bool side( index_t i ) const
        {
            return sides_[i];
        }
        const VolumeMesh< DIMENSION >& mesh() const
        {
            return *volume_mesh_;
        }

    private:
        Region( const GeoModel< DIMENSION >& geomodel,
            index_t id,
            const MeshType type )
            : GeoModelMeshEntity< DIMENSION >( geomodel, id )
        {
            update_mesh_storage_type(
                VolumeMesh< DIMENSION >::create_mesh( type ) );
        }

        bool is_mesh_valid() const final;

        void update_mesh_storage_type(
            std::unique_ptr< VolumeMesh< DIMENSION > > mesh )
        {
            volume_mesh_ = std::move( mesh );
            GeoModelMeshEntity< DIMENSION >::set_mesh( volume_mesh_ );
        }

        void change_mesh_data_structure( const MeshType& type ) final;

        void copy_mesh_entity(
            const GeoModelMeshEntity< DIMENSION >& from ) final
        {
            const auto& region_from =
                dynamic_cast< const Region< DIMENSION >& >( from );
            GeoModelMeshEntity< DIMENSION >::copy_mesh_entity( from );
            sides_ = region_from.sides_;
        }

    protected:
        std::vector< bool > sides_{};

    private:
        std::shared_ptr< VolumeMesh< DIMENSION > > volume_mesh_{};
    };

    ALIAS_3D( Region );

    template < index_t DIMENSION >
    class RINGMESH_API GeoModelMeshEntityConstAccess
    {
        ringmesh_disable_copy_and_move( GeoModelMeshEntityConstAccess );
        ringmesh_template_assert_2d_or_3d( DIMENSION );
        friend class GeoModelBuilderGeometryBase< DIMENSION >;
        friend class GeoModelBuilderGeometry< DIMENSION >;
        friend class GeoModelBuilderTopologyBase< DIMENSION >;
        friend class GeoModelBuilderTopology< DIMENSION >;

    private:
        explicit GeoModelMeshEntityConstAccess(
            const GeoModelMeshEntity< DIMENSION >& gme )
            : gmme_( gme )
        {
        }

        ~GeoModelMeshEntityConstAccess() = default;

        const std::shared_ptr< MeshBase< DIMENSION > >& mesh() const
        {
            return gmme_.mesh_;
        }

        const std::vector< index_t >& incident_entity_relation_ids() const
        {
            return gmme_.incident_entities_;
        }

        const std::vector< index_t >& boundary_relation_ids() const
        {
            return gmme_.boundaries_;
        }

    private:
        const GeoModelMeshEntity< DIMENSION >& gmme_;
    };

    template < index_t DIMENSION >
    class RINGMESH_API GeoModelMeshEntityAccess
    {
        ringmesh_disable_copy_and_move( GeoModelMeshEntityAccess );
        ringmesh_template_assert_2d_or_3d( DIMENSION );
        friend class GeoModelBuilderTopologyBase< DIMENSION >;
        friend class GeoModelBuilderTopology< DIMENSION >;
        friend class GeoModelBuilderGeometryBase< DIMENSION >;
        friend class GeoModelBuilderGeometry< DIMENSION >;
        friend class GeoModelBuilderGeology< DIMENSION >;
        friend class GeoModelBuilderInfo< DIMENSION >;
        friend class GeoModelBuilderRemovalBase< DIMENSION >;
        friend class GeoModelBuilderRemoval< DIMENSION >;

    private:
        explicit GeoModelMeshEntityAccess(
            GeoModelMeshEntity< DIMENSION >& gme )
            : gmme_( gme )
        {
        }

        ~GeoModelMeshEntityAccess() = default;

        std::string& modifiable_name()
        {
            return gmme_.name_;
        }

        index_t& modifiable_index()
        {
            return gmme_.id_;
        }

        std::vector< index_t >& modifiable_boundaries()
        {
            return gmme_.boundaries_;
        }

        std::vector< index_t >& modifiable_incident_entities()
        {
            return gmme_.incident_entities_;
        }

        std::vector< bool >& modifiable_sides();

        std::vector< index_t >& modifiable_parents()
        {
            return gmme_.parents_;
        }

        std::shared_ptr< MeshBase< DIMENSION > >& modifiable_mesh()
        {
            return gmme_.mesh_;
        }

        void copy( const GeoModelMeshEntity< DIMENSION >& from )
        {
            gmme_.copy_mesh_entity( from );
        }

        void change_mesh_data_structure( const MeshType& type );

        template < template < index_t > class ENTITY >
        static std::unique_ptr< ENTITY< DIMENSION > > create_entity(
            const GeoModel< DIMENSION >& geomodel,
            index_t id,
            const MeshType& type )
        {
            return std::unique_ptr< ENTITY< DIMENSION > >(
                new ENTITY< DIMENSION >( geomodel, id, type ) );
        }

    private:
        GeoModelMeshEntity< DIMENSION >& gmme_;
    };
    ALIAS_2D_AND_3D( GeoModelMeshEntityAccess );
} // namespace RINGMesh