mesh.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 <algorithm>
#include <memory>

#include <ringmesh/basic/factory.h>
#include <ringmesh/basic/nn_search.h>

#include <ringmesh/mesh/aabb.h>

namespace GEO
{
    class AttributesManager;
}

namespace RINGMesh
{
    FORWARD_DECLARATION_DIMENSION_CLASS( GeoModel );
    FORWARD_DECLARATION_DIMENSION_CLASS( MeshBaseBuilder );
    FORWARD_DECLARATION_DIMENSION_CLASS( PointSetMeshBuilder );
    FORWARD_DECLARATION_DIMENSION_CLASS( LineMeshBuilder );
    FORWARD_DECLARATION_DIMENSION_CLASS( SurfaceMeshBuilder );
    FORWARD_DECLARATION_DIMENSION_CLASS( VolumeMeshBuilder );
    FORWARD_DECLARATION_DIMENSION_CLASS( SurfaceMesh );
    struct EdgeLocalVertex;
    struct PolygonLocalEdge;
    struct CellLocalFacet;
} // namespace RINGMesh

namespace RINGMesh
{
    using MeshType = std::string;

    struct RINGMESH_API ElementLocalVertex
    {
        ElementLocalVertex() = default;
        ElementLocalVertex( index_t element_id, index_t local_vertex_id )
            : element_id_( element_id ), local_vertex_id_( local_vertex_id )
        {
        }
        ElementLocalVertex( EdgeLocalVertex edge_local_vertex );
        ElementLocalVertex( PolygonLocalEdge polygon_local_edge );
        ElementLocalVertex( CellLocalFacet cell_local_facet );
        index_t element_id_{ NO_ID };
        index_t local_vertex_id_{ NO_ID };
    };

    struct RINGMESH_API EdgeLocalVertex
    {
        EdgeLocalVertex() = default;
        EdgeLocalVertex( index_t edge_id, index_t local_vertex_id )
            : edge_id_( edge_id ), local_vertex_id_( local_vertex_id )
        {
        }
        EdgeLocalVertex( ElementLocalVertex edge_local_vertex )
            : edge_id_( std::move( edge_local_vertex.element_id_ ) ),
              local_vertex_id_(
                  std::move( edge_local_vertex.local_vertex_id_ ) )
        {
        }
        index_t edge_id_{ NO_ID };
        index_t local_vertex_id_{ NO_ID };
    };

    struct RINGMESH_API PolygonLocalEdge
    {
        PolygonLocalEdge() = default;
        PolygonLocalEdge( index_t polygon_id, index_t local_edge_id )
            : polygon_id_( polygon_id ), local_edge_id_( local_edge_id )
        {
        }
        PolygonLocalEdge( ElementLocalVertex polygon_local_vertex )
            : polygon_id_( std::move( polygon_local_vertex.element_id_ ) ),
              local_edge_id_(
                  std::move( polygon_local_vertex.local_vertex_id_ ) )
        {
        }
        index_t polygon_id_{ NO_ID };
        index_t local_edge_id_{ NO_ID };
    };

    struct RINGMESH_API CellLocalFacet
    {
        CellLocalFacet() = default;
        CellLocalFacet( index_t cell_id, index_t local_facet_id )
            : cell_id_( cell_id ), local_facet_id_( local_facet_id )
        {
        }
        index_t cell_id_{ NO_ID };
        index_t local_facet_id_{ NO_ID };
    };

    template < index_t DIMENSION >
    class MeshBase
    {
        ringmesh_disable_copy_and_move( MeshBase );
        ringmesh_template_assert_2d_or_3d( DIMENSION );
        friend class MeshBaseBuilder< DIMENSION >;

    public:
        virtual ~MeshBase() = default;

        virtual void save_mesh( const std::string& filename ) const = 0;

        virtual std::tuple< index_t, std::vector< index_t > >
            connected_components() const = 0;

        // TODO maybe reimplement the function with a RINGMesh::Mesh??
        virtual void print_mesh_bounded_attributes() const = 0;
        virtual const vecn< DIMENSION >& vertex( index_t v_id ) const = 0;
        /*
         * @brief Gets the number of vertices in the Mesh.
         */
        virtual index_t nb_vertices() const = 0;

        virtual GEO::AttributesManager& vertex_attribute_manager() const = 0;

        const NNSearch< DIMENSION >& vertex_nn_search() const
        {
            if( !vertex_nn_search_ )
            {
                std::vector< vecn< DIMENSION > > vec_vertices( nb_vertices() );
                for( auto v : range( nb_vertices() ) )
                {
                    vec_vertices[v] = vertex( v );
                }
                vertex_nn_search_.reset(
                    new NNSearch< DIMENSION >( vec_vertices, true ) );
            }
            return *vertex_nn_search_.get();
        }

        virtual MeshType type_name() const = 0;

        virtual std::string default_extension() const = 0;

        virtual bool is_mesh_valid() const = 0;

    protected:
        MeshBase() = default;

    protected:
        mutable std::unique_ptr< NNSearch< DIMENSION > > vertex_nn_search_{};
    };
    ALIAS_2D_AND_3D( MeshBase );

    template < index_t DIMENSION >
    class PointSetMesh : public MeshBase< DIMENSION >
    {
        friend class PointSetMeshBuilder< DIMENSION >;

    public:
        static std::unique_ptr< PointSetMesh< DIMENSION > > create_mesh(
            const MeshType type = "" );
        std::tuple< index_t, std::vector< index_t > >
            connected_components() const final;
        bool is_mesh_valid() const override
        {
            return true;
        }

    protected:
        PointSetMesh() = default;
    };
    ALIAS_2D_AND_3D( PointSetMesh );

    template < index_t DIMENSION >
    using PointSetMeshFactory = Factory< MeshType, PointSetMesh< DIMENSION > >;
    ALIAS_2D_AND_3D( PointSetMeshFactory );

    template < index_t DIMENSION >
    class LineMesh : public MeshBase< DIMENSION >
    {
        friend class LineMeshBuilder< DIMENSION >;

    public:
        static std::unique_ptr< LineMesh< DIMENSION > > create_mesh(
            const MeshType type = "" );

        /*
         * @brief Gets the index of an edge vertex.
         * @param[in] edge_local_vertex index of the edge and of the
         * local index of the vertex, in {0,1}
         * @return the global index of vertex in \p edge_local_vertex.
         */
        virtual index_t edge_vertex(
            const ElementLocalVertex& edge_local_vertex ) const = 0;

        virtual index_t nb_edges() const = 0;

        double edge_length( index_t edge_id ) const
        {
            const vecn< DIMENSION >& e0 =
                this->vertex( edge_vertex( ElementLocalVertex( edge_id, 0 ) ) );
            const vecn< DIMENSION >& e1 =
                this->vertex( edge_vertex( ElementLocalVertex( edge_id, 1 ) ) );
            return ( e1 - e0 ).length();
        }

        vecn< DIMENSION > edge_barycenter( index_t edge_id ) const
        {
            const vecn< DIMENSION >& e0 =
                this->vertex( edge_vertex( ElementLocalVertex( edge_id, 0 ) ) );
            const vecn< DIMENSION >& e1 =
                this->vertex( edge_vertex( ElementLocalVertex( edge_id, 1 ) ) );
            return ( e1 + e0 ) / 2.;
        }

        const NNSearch< DIMENSION >& edge_nn_search() const
        {
            if( !edge_nn_search_ )
            {
                std::vector< vecn< DIMENSION > > edge_centers( nb_edges() );
                for( auto e : range( nb_edges() ) )
                {
                    edge_centers[e] = edge_barycenter( e );
                }
                edge_nn_search_.reset(
                    new NNSearch< DIMENSION >( edge_centers, true ) );
            }
            return *edge_nn_search_.get();
        }
        const LineAABBTree< DIMENSION >& edge_aabb() const
        {
            if( !edge_aabb_ )
            {
                edge_aabb_.reset( new LineAABBTree< DIMENSION >( *this ) );
            }
            return *edge_aabb_.get();
        }

        virtual GEO::AttributesManager& edge_attribute_manager() const = 0;

        bool is_mesh_valid() const override
        {
            bool valid{ true };

            if( this->nb_vertices() < 2 )
            {
                Logger::err( "LineMesh", "Mesh has less than 2 vertices " );
                valid = false;
            }

            if( nb_edges() == 0 )
            {
                Logger::err( "LineMesh", "Mesh has no edge" );
                valid = false;
            }

            // No isolated vertices
            std::vector< index_t > nb( this->nb_vertices(), 0 );
            for( auto p : range( nb_edges() ) )
            {
                for( auto v : range( 2 ) )
                {
                    nb[edge_vertex( { p, v } )]++;
                }
            }
            auto nb_isolated_vertices =
                static_cast< index_t >( std::count( nb.begin(), nb.end(), 0 ) );
            if( nb_isolated_vertices > 0 )
            {
                Logger::warn( "LineMesh", "Mesh has ", nb_isolated_vertices,
                    " isolated vertices " );
                valid = false;
            }

            return valid;
        }
        std::tuple< index_t, std::vector< index_t > >
            connected_components() const final;

    protected:
        LineMesh() = default;

    protected:
        mutable std::unique_ptr< NNSearch< DIMENSION > > edge_nn_search_{};
        mutable std::unique_ptr< LineAABBTree< DIMENSION > > edge_aabb_{};
    };
    ALIAS_2D_AND_3D( LineMesh );

    template < index_t DIMENSION >
    using LineMeshFactory = Factory< MeshType, LineMesh< DIMENSION > >;
    ALIAS_2D_AND_3D( LineMeshFactory );

    template < index_t DIMENSION >
    class SurfaceMeshBase : public MeshBase< DIMENSION >
    {
        friend class SurfaceMeshBuilder< DIMENSION >;

    public:
        static std::unique_ptr< SurfaceMesh< DIMENSION > > create_mesh(
            const MeshType type = "" );

        virtual index_t polygon_vertex(
            const ElementLocalVertex& polygon_local_vertex ) const = 0;

        virtual index_t nb_polygons() const = 0;
        virtual index_t nb_polygon_vertices( index_t polygon_id ) const = 0;
        ElementLocalVertex next_polygon_vertex(
            const ElementLocalVertex& polygon_local_vertex ) const
        {
            const index_t local_vertex_id =
                polygon_local_vertex.local_vertex_id_;
            ringmesh_assert(
                local_vertex_id
                < nb_polygon_vertices( polygon_local_vertex.element_id_ ) );
            if( local_vertex_id
                != nb_polygon_vertices( polygon_local_vertex.element_id_ ) - 1 )
            {
                return { polygon_local_vertex.element_id_,
                    local_vertex_id + 1 };
            }
            return { polygon_local_vertex.element_id_, 0 };
        }
        PolygonLocalEdge next_on_border(
            const PolygonLocalEdge& polygon_local_edge ) const;

        ElementLocalVertex prev_polygon_vertex(
            const ElementLocalVertex& polygon_local_vertex ) const
        {
            ringmesh_assert(
                polygon_local_vertex.local_vertex_id_
                < nb_polygon_vertices( polygon_local_vertex.element_id_ ) );
            if( polygon_local_vertex.local_vertex_id_ > 0 )
            {
                return { polygon_local_vertex.element_id_,
                    polygon_local_vertex.local_vertex_id_ - 1 };
            }
            return { polygon_local_vertex.element_id_,
                nb_polygon_vertices( polygon_local_vertex.element_id_ ) - 1 };
        }

        PolygonLocalEdge prev_on_border(
            const PolygonLocalEdge& polygon_local_edge ) const;

        index_t vertex_index_in_polygon(
            index_t polygon_index, index_t vertex_id ) const;

        index_t closest_vertex_in_polygon(
            index_t polygon_index, const vecn< DIMENSION >& query_point ) const;

        index_t polygon_from_vertex_ids( index_t in0, index_t in1 ) const;

        std::vector< index_t > polygons_around_vertex(
            index_t vertex_id, bool border_only, index_t first_polygon ) const;

        virtual index_t polygon_adjacent(
            const PolygonLocalEdge& polygon_local_edge ) const = 0;

        virtual GEO::AttributesManager& polygon_attribute_manager() const = 0;
        virtual bool polygons_are_simplicies() const = 0;
        bool is_triangle( index_t polygon_id ) const
        {
            return nb_polygon_vertices( polygon_id ) == 3;
        }

        PolygonType polygone_type( index_t polygon_id ) const
        {
            if( is_triangle( polygon_id ) )
            {
                return PolygonType::TRIANGLE;
            }
            if( nb_polygon_vertices( polygon_id ) == 4 )
            {
                return PolygonType::QUAD;
            }
            return PolygonType::UNDEFINED;
        }

        bool is_edge_on_border(
            const PolygonLocalEdge& polygon_local_edge ) const
        {
            return polygon_adjacent( polygon_local_edge ) == NO_ID;
        }

        bool is_polygon_on_border( index_t polygon_index ) const
        {
            for( auto v : range( nb_polygon_vertices( polygon_index ) ) )
            {
                if( is_edge_on_border( PolygonLocalEdge( polygon_index, v ) ) )
                {
                    return true;
                }
            }
            return false;
        }

        double polygon_edge_length(
            const PolygonLocalEdge& polygon_local_edge ) const
        {
            const vecn< DIMENSION >& e0 =
                this->vertex( polygon_edge_vertex( polygon_local_edge, 0 ) );
            const vecn< DIMENSION >& e1 =
                this->vertex( polygon_edge_vertex( polygon_local_edge, 1 ) );
            return ( e1 - e0 ).length();
        }
        vecn< DIMENSION > polygon_edge_barycenter(
            const PolygonLocalEdge& polygon_local_edge ) const
        {
            const vecn< DIMENSION >& e0 =
                this->vertex( polygon_edge_vertex( polygon_local_edge, 0 ) );
            const vecn< DIMENSION >& e1 =
                this->vertex( polygon_edge_vertex( polygon_local_edge, 1 ) );
            return ( e1 + e0 ) / 2.;
        }
        index_t polygon_edge_vertex( const PolygonLocalEdge& polygon_local_edge,
            index_t vertex_id ) const
        {
            ringmesh_assert( vertex_id < 2 );
            if( vertex_id == 0 )
            {
                return polygon_vertex( polygon_local_edge );
            }
            return polygon_vertex( ElementLocalVertex(
                polygon_local_edge.polygon_id_,
                ( polygon_local_edge.local_edge_id_ + vertex_id )
                    % nb_polygon_vertices( polygon_local_edge.polygon_id_ ) ) );
        }

        vecn< DIMENSION > polygon_barycenter( index_t polygon_id ) const
        {
            vecn< DIMENSION > result;
            ringmesh_assert( nb_polygon_vertices( polygon_id ) >= 1 );
            for( auto v : range( nb_polygon_vertices( polygon_id ) ) )
            {
                result += this->vertex(
                    polygon_vertex( ElementLocalVertex( polygon_id, v ) ) );
            }
            return ( 1.0 / nb_polygon_vertices( polygon_id ) ) * result;
        }
        virtual double polygon_area( index_t polygon_id ) const = 0;

        const NNSearch< DIMENSION >& polygon_nn_search() const
        {
            if( !nn_search_ )
            {
                std::vector< vecn< DIMENSION > > polygon_centers(
                    nb_polygons() );
                for( auto p : range( nb_polygons() ) )
                {
                    polygon_centers[p] = polygon_barycenter( p );
                }
                nn_search_.reset(
                    new NNSearch< DIMENSION >( polygon_centers, true ) );
            }
            return *nn_search_.get();
        }
        const SurfaceAABBTree< DIMENSION >& polygon_aabb() const
        {
            if( !polygon_aabb_ )
            {
                polygon_aabb_.reset(
                    new SurfaceAABBTree< DIMENSION >( *this ) );
            }
            return *polygon_aabb_;
        }

        bool is_mesh_valid() const override
        {
            bool valid{ true };

            if( this->nb_vertices() < 3 )
            {
                Logger::warn( "SurfaceMesh has less than 3 vertices " );
                valid = false;
            }
            if( nb_polygons() == 0 )
            {
                Logger::warn( "SurfaceMesh has no polygon" );
                valid = false;
            }

            // No isolated vertices
            std::vector< index_t > nb( this->nb_vertices(), 0 );
            for( auto p : range( nb_polygons() ) )
            {
                for( auto v : range( nb_polygon_vertices( p ) ) )
                {
                    nb[polygon_vertex( { p, v } )]++;
                }
            }
            auto nb_isolated_vertices =
                static_cast< index_t >( std::count( nb.begin(), nb.end(), 0 ) );
            if( nb_isolated_vertices > 0 )
            {
                Logger::warn( "SurfaceMesh", "Mesh has ", nb_isolated_vertices,
                    " isolated vertices " );
                valid = false;
            }

            return valid;
        }
        std::tuple< index_t, std::vector< index_t > >
            connected_components() const final;

    protected:
        SurfaceMeshBase() = default;

    protected:
        mutable std::unique_ptr< NNSearch< DIMENSION > > nn_search_{};
        mutable std::unique_ptr< SurfaceAABBTree< DIMENSION > > polygon_aabb_{};
    };
    ALIAS_2D_AND_3D( SurfaceMeshBase );

    template < index_t DIMENSION >
    class SurfaceMesh : public SurfaceMeshBase< DIMENSION >
    {
    };

    template < index_t DIMENSION >
    using SurfaceMeshFactory = Factory< MeshType, SurfaceMesh< DIMENSION > >;
    ALIAS_2D_AND_3D( SurfaceMeshFactory );

    template <>
    class RINGMESH_API SurfaceMesh< 3 > : public SurfaceMeshBase< 3 >
    {
    public:
        double polygon_area( index_t polygon_id ) const;

        vec3 polygon_normal( index_t polygon_id ) const
        {
            const vec3& p1 = this->vertex(
                this->polygon_vertex( ElementLocalVertex( polygon_id, 0 ) ) );
            const vec3& p2 = this->vertex(
                this->polygon_vertex( ElementLocalVertex( polygon_id, 1 ) ) );
            const vec3& p3 = this->vertex(
                this->polygon_vertex( ElementLocalVertex( polygon_id, 2 ) ) );
            vec3 norm = cross( p2 - p1, p3 - p1 );
            return normalize( norm );
        }

        vec3 normal_at_vertex( index_t vertex_id, index_t p0 = NO_ID ) const
        {
            ringmesh_assert( vertex_id < nb_vertices() );
            index_t p = 0;
            while( p0 == NO_ID && p < nb_polygons() )
            {
                for( auto lv : range( nb_polygon_vertices( p ) ) )
                {
                    if( polygon_vertex( ElementLocalVertex( p, lv ) )
                        == vertex_id )
                    {
                        p0 = p;
                        break;
                    }
                }
                p++;
            }

            std::vector< index_t > polygon_ids =
                polygons_around_vertex( vertex_id, false, p0 );
            vec3 norm;
            for( auto polygon_id : polygon_ids )
            {
                norm += polygon_normal( polygon_id );
            }
            return normalize( norm );
        }
    };

    template <>
    class RINGMESH_API SurfaceMesh< 2 > : public SurfaceMeshBase< 2 >
    {
    public:
        double polygon_area( index_t polygon_id ) const override
        {
            double result = 0.0;
            if( nb_polygon_vertices( polygon_id ) == 0 )
            {
                return result;
            }
            const vec2& p1 =
                vertex( polygon_vertex( ElementLocalVertex( polygon_id, 0 ) ) );
            for( auto i : range( 1, nb_polygon_vertices( polygon_id ) - 1 ) )
            {
                const vec2& p2 = vertex(
                    polygon_vertex( ElementLocalVertex( polygon_id, i ) ) );
                const vec2& p3 = vertex(
                    polygon_vertex( ElementLocalVertex( polygon_id, i + 1 ) ) );
                result += GEO::Geom::triangle_signed_area( p1, p2, p3 );
            }
            return std::fabs( result );
        }
    };

    ALIAS_2D_AND_3D( SurfaceMesh );

    template < index_t DIMENSION >
    class VolumeMesh : public MeshBase< DIMENSION >
    {
        ringmesh_template_assert_3d( DIMENSION );
        friend class VolumeMeshBuilder< DIMENSION >;

    public:
        static std::unique_ptr< VolumeMesh< DIMENSION > > create_mesh(
            const MeshType type = "" );

        virtual index_t cell_vertex(
            const ElementLocalVertex& cell_local_vertex ) const = 0;

        virtual index_t cell_edge_vertex(
            index_t cell_id, index_t edge_id, index_t vertex_id ) const = 0;

        virtual index_t cell_facet_vertex(
            const CellLocalFacet& cell_local_facet,
            index_t vertex_id ) const = 0;

        virtual index_t cell_facet(
            const CellLocalFacet& cell_local_facet ) const = 0;

        double cell_edge_length( index_t cell_id, index_t edge_id ) const
        {
            const vecn< DIMENSION >& e0 =
                this->vertex( cell_edge_vertex( cell_id, edge_id, 0 ) );
            const vecn< DIMENSION >& e1 =
                this->vertex( cell_edge_vertex( cell_id, edge_id, 1 ) );
            return ( e1 - e0 ).length();
        }

        vecn< DIMENSION > cell_edge_barycenter(
            index_t cell_id, index_t edge_id ) const
        {
            const vecn< DIMENSION >& e0 =
                this->vertex( cell_edge_vertex( cell_id, edge_id, 0 ) );
            const vecn< DIMENSION >& e1 =
                this->vertex( cell_edge_vertex( cell_id, edge_id, 1 ) );
            return ( e1 + e0 ) / 2.;
        }

        virtual index_t nb_cell_facets( index_t cell_id ) const = 0;
        virtual index_t nb_cell_facets() const = 0;

        virtual index_t nb_cell_edges( index_t cell_id ) const = 0;

        virtual index_t nb_cell_facet_vertices(
            const CellLocalFacet& cell_local_facet ) const = 0;

        virtual index_t nb_cell_vertices( index_t cell_id ) const = 0;

        virtual index_t nb_cells() const = 0;

        virtual index_t cell_begin( index_t cell_id ) const = 0;

        virtual index_t cell_end( index_t cell_id ) const = 0;

        virtual index_t cell_adjacent(
            const CellLocalFacet& cell_local_facet ) const = 0;

        virtual GEO::AttributesManager& cell_attribute_manager() const = 0;

        virtual GEO::AttributesManager&
            cell_facet_attribute_manager() const = 0;

        virtual CellType cell_type( index_t cell_id ) const = 0;

        virtual bool cells_are_simplicies() const = 0;

        vecn< DIMENSION > cell_facet_barycenter(
            const CellLocalFacet& cell_local_facet ) const
        {
            vecn< DIMENSION > result;
            index_t nb_vertices = nb_cell_facet_vertices( cell_local_facet );
            for( auto v : range( nb_vertices ) )
            {
                result +=
                    this->vertex( cell_facet_vertex( cell_local_facet, v ) );
            }
            ringmesh_assert( nb_vertices > 0 );

            return result / static_cast< double >( nb_vertices );
        }
        vecn< DIMENSION > cell_barycenter( index_t cell_id ) const
        {
            vecn< DIMENSION > result;
            ringmesh_assert( nb_cell_vertices( cell_id ) >= 1 );
            for( auto v : range( nb_cell_vertices( cell_id ) ) )
            {
                result += this->vertex(
                    cell_vertex( ElementLocalVertex( cell_id, v ) ) );
            }
            return ( 1.0 / nb_cell_vertices( cell_id ) ) * result;
        }
        vecn< DIMENSION > cell_facet_normal(
            const CellLocalFacet& cell_local_facet ) const
        {
            ringmesh_assert( cell_local_facet.cell_id_ < nb_cells() );
            ringmesh_assert( cell_local_facet.local_facet_id_
                             < nb_cell_facets( cell_local_facet.cell_id_ ) );

            const vecn< DIMENSION >& p1 =
                this->vertex( cell_facet_vertex( cell_local_facet, 0 ) );
            const vecn< DIMENSION >& p2 =
                this->vertex( cell_facet_vertex( cell_local_facet, 1 ) );
            const vecn< DIMENSION >& p3 =
                this->vertex( cell_facet_vertex( cell_local_facet, 2 ) );

            return cross( p2 - p1, p3 - p1 );
        }

        virtual double cell_volume( index_t cell_id ) const = 0;

        std::vector< index_t > cells_around_vertex(
            index_t vertex_id, index_t cell_hint ) const;

        index_t find_cell_corner( index_t cell_id, index_t vertex_id ) const
        {
            for( auto v : range( nb_cell_vertices( cell_id ) ) )
            {
                if( cell_vertex( ElementLocalVertex( cell_id, v ) )
                    == vertex_id )
                {
                    return v;
                }
            }
            return NO_ID;
        }

        bool find_cell_from_colocated_vertex_within_distance_if_any(
            const vecn< DIMENSION >& vertex_vec,
            double distance,
            index_t& cell_id,
            index_t& cell_vertex_id ) const;

        const NNSearch< DIMENSION >& cell_facet_nn_search() const
        {
            if( !cell_facet_nn_search_ )
            {
                std::vector< vecn< DIMENSION > > cell_facet_centers(
                    nb_cell_facets() );
                index_t cf = 0;
                for( auto c : range( nb_cells() ) )
                {
                    for( auto f : range( nb_cell_facets( c ) ) )
                    {
                        cell_facet_centers[cf] =
                            cell_facet_barycenter( CellLocalFacet( c, f ) );
                        ++cf;
                    }
                }
                cell_facet_nn_search_.reset(
                    new NNSearch< DIMENSION >( cell_facet_centers, true ) );
            }
            return *cell_facet_nn_search_.get();
        }
        const NNSearch< DIMENSION >& cell_nn_search() const
        {
            if( !cell_nn_search_ )
            {
                std::vector< vecn< DIMENSION > > cell_centers( nb_cells() );
                for( auto c : range( nb_cells() ) )
                {
                    cell_centers[c] = cell_barycenter( c );
                }
                cell_nn_search_.reset(
                    new NNSearch< DIMENSION >( cell_centers, true ) );
            }
            return *cell_nn_search_.get();
        }
        const VolumeAABBTree< DIMENSION >& cell_aabb() const
        {
            if( !cell_aabb_ )
            {
                cell_aabb_.reset( new VolumeAABBTree< DIMENSION >( *this ) );
            }
            return *cell_aabb_.get();
        }

        bool is_mesh_valid() const override
        {
            bool valid{ true };

            if( this->nb_vertices() < 4 )
            {
                Logger::warn( "VolumeMesh has less than 4 vertices " );
                valid = false;
            }
            if( nb_cells() == 0 )
            {
                Logger::warn( "VolumeMesh has no cell" );
                valid = false;
            }

            // No isolated vertices
            std::vector< index_t > nb( this->nb_vertices(), 0 );
            for( auto c : range( nb_cells() ) )
            {
                for( auto v : range( nb_cell_vertices( c ) ) )
                {
                    nb[cell_vertex( { c, v } )]++;
                }
            }
            auto nb_isolated_vertices =
                static_cast< index_t >( std::count( nb.begin(), nb.end(), 0 ) );
            if( nb_isolated_vertices > 0 )
            {
                Logger::warn( "VolumeMesh", "Mesh has ", nb_isolated_vertices,
                    " isolated vertices " );
                valid = false;
            }

            return valid;
        }
        std::tuple< index_t, std::vector< index_t > >
            connected_components() const final;

    protected:
        VolumeMesh() = default;

    protected:
        mutable std::unique_ptr< NNSearch< DIMENSION > >
            cell_facet_nn_search_{};
        mutable std::unique_ptr< NNSearch< DIMENSION > > cell_nn_search_{};
        mutable std::unique_ptr< VolumeAABBTree< DIMENSION > > cell_aabb_{};
    };

    using VolumeMesh3D = VolumeMesh< 3 >;

    template < index_t DIMENSION >
    using VolumeMeshFactory = Factory< MeshType, VolumeMesh< DIMENSION > >;
    using VolumeMeshFactory3D = VolumeMeshFactory< 3 >;

    template < index_t DIMENSION >
    class MeshSetBase
    {
        ringmesh_disable_copy_and_move( MeshSetBase );
        ringmesh_template_assert_2d_or_3d( DIMENSION );

    public:
        void create_point_set_mesh( MeshType type );
        void create_line_mesh( MeshType type );
        void create_well_mesh( MeshType type );
        void create_surface_mesh( MeshType type );

    protected:
        MeshSetBase();
        virtual ~MeshSetBase() = default;

    public:
        std::unique_ptr< PointSetMesh< DIMENSION > > point_set_mesh{};
        std::unique_ptr< LineMesh< DIMENSION > > well_mesh{};
        std::unique_ptr< LineMesh< DIMENSION > > line_mesh{};
        std::unique_ptr< SurfaceMesh< DIMENSION > > surface_mesh{};
    };

    template < index_t DIMENSION >
    class MeshSet : public MeshSetBase< DIMENSION >
    {
    public:
        MeshSet() = default;
    };

    template <>
    class RINGMESH_API MeshSet< 3 > : public MeshSetBase< 3 >
    {
    public:
        MeshSet();

        void create_volume_mesh( MeshType type );

    public:
        std::unique_ptr< VolumeMesh3D > volume_mesh{};
    };
} // namespace RINGMesh