Commit 551161cd authored by Ján Bobot's avatar Ján Bobot Committed by Jakub Klinkovský
Browse files

Meshes/Geometry: refactored polygon and polyhedron decomposition into... 

Meshes/Geometry: refactored polygon and polyhedron decomposition into templated EntityDecomposer struct with specializations
parent f1389ef8

src/TNL/Meshes/Geometry/EntityDecomposer.h

0 → 100644
+235 −0
Original line number Diff line number Diff line 
#pragma once



#include <TNL/Meshes/Mesh.h>

#include <TNL/Meshes/MeshEntity.h>

#include <TNL/Meshes/MeshBuilder.h>

#include <TNL/Meshes/Topologies/Triangle.h>

#include <TNL/Meshes/Topologies/Polygon.h>

#include <TNL/Meshes/Geometry/getEntityCenter.h>

#include <TNL/Meshes/Geometry/getEntityMeasure.h>

#include <functional>



namespace TNL {

namespace Meshes {



enum class EntityDecomposerVersion

{

   ConnectEdgesToCentroid, ConnectEdgesToPoint

};



template< typename MeshConfig,

          typename Topology,

          EntityDecomposerVersion EntityDecomposerVersion_ = EntityDecomposerVersion::ConnectEdgesToCentroid, 

          EntityDecomposerVersion SubentityDecomposerVersion = EntityDecomposerVersion::ConnectEdgesToCentroid >

struct EntityDecomposer;



// Polygon

template< typename MeshConfig, EntityDecomposerVersion SubentityDecomposerVersion > // SubentityDecomposerVersion is not used for polygons

struct EntityDecomposer< MeshConfig, Topologies::Polygon, EntityDecomposerVersion::ConnectEdgesToCentroid, SubentityDecomposerVersion >

{

   using Device = typename Devices::Host;

   using Topology = typename Topologies::Polygon;

   using MeshEntityType = MeshEntity< MeshConfig, Device, Topology >;

   using VertexMeshEntityType = typename MeshEntityType::template SubentityTraits< 0 >::SubentityType;

   using PointType = typename VertexMeshEntityType::PointType;

   using GlobalIndexType = typename MeshConfig::GlobalIndexType;

   using LocalIndexType = typename MeshConfig::LocalIndexType;

   using PointCreationFunctorType = std::function< GlobalIndexType ( const PointType& ) >;

   using DecomposedEntityFunctorType = std::function< void ( GlobalIndexType, GlobalIndexType, GlobalIndexType ) >;

   

   static std::pair< GlobalIndexType, GlobalIndexType > getExtraPointsAndEntitiesCount( const MeshEntityType & entity )

   {

      const auto pointsCount = entity.template getSubentitiesCount< 0 >();

      if( pointsCount == 3 ) // polygon is triangle

         return { 0, 1 }; // No extra points and no decomposition

      return { 1, pointsCount }; // 1 extra centroid point and decomposition creates pointsCount triangles

   }



   static void decompose( const MeshEntityType & entity,

                          PointCreationFunctorType pointCreationFunctor,

                          DecomposedEntityFunctorType decomposedEntityFunctor )

   {

      const auto verticesCount = entity.template getSubentitiesCount< 0 >();

      if( verticesCount == 3 ) { // polygon is only copied as it's already a triangle

         const auto v0 = entity.template getSubentityIndex< 0 >( 0 );

         const auto v1 = entity.template getSubentityIndex< 0 >( 1 );

         const auto v2 = entity.template getSubentityIndex< 0 >( 2 );

         decomposedEntityFunctor( v0, v1, v2 );

      }

      else { // polygon is decomposed as it has got more than 3 vertices

         const auto v0 = pointCreationFunctor( getEntityCenter( entity.getMesh(), entity ) );

         // decompose polygon into triangles by connecting each edge to the centroid

         for( LocalIndexType j = 0, k = 1; k < verticesCount; j++, k++ ) {

            const auto v1 = entity.template getSubentityIndex< 0 >( j );

            const auto v2 = entity.template getSubentityIndex< 0 >( k );

            decomposedEntityFunctor( v0, v1, v2 );

         }

         { // wrap around term

            const auto v1 = entity.template getSubentityIndex< 0 >( verticesCount - 1 );

            const auto v2 = entity.template getSubentityIndex< 0 >( 0 );

            decomposedEntityFunctor( v0, v1, v2 );

         }

      }

   }

};



template< typename MeshConfig, EntityDecomposerVersion SubentityDecomposerVersion > // SubentityDecomposerVersion is not used for polygons

struct EntityDecomposer< MeshConfig, Topologies::Polygon, EntityDecomposerVersion::ConnectEdgesToPoint, SubentityDecomposerVersion >

{

   using Device = typename Devices::Host;

   using Topology = typename Topologies::Polygon;

   using MeshEntityType = MeshEntity< MeshConfig, Device, Topology >;

   using VertexMeshEntityType = typename MeshEntityType::template SubentityTraits< 0 >::SubentityType;

   using PointType = typename VertexMeshEntityType::PointType;

   using GlobalIndexType = typename MeshConfig::GlobalIndexType;

   using LocalIndexType = typename MeshConfig::LocalIndexType;

   using PointCreationFunctorType = std::function< GlobalIndexType ( const PointType& ) >;

   using DecomposedEntityFunctorType = std::function< void ( GlobalIndexType, GlobalIndexType, GlobalIndexType ) >;



   static std::pair< GlobalIndexType, GlobalIndexType > getExtraPointsAndEntitiesCount( const MeshEntityType & entity )

   {

      const auto pointsCount = entity.template getSubentitiesCount< 0 >();

      return { 0, pointsCount - 2 }; // no extra points and there is a triangle for every non-adjacent edge to the 0th point (pointsCount - 2)

   }



   static void decompose( const MeshEntityType & entity,

                          PointCreationFunctorType pointCreationFunctor,

                          DecomposedEntityFunctorType decomposedEntityFunctor )

   {

      // decompose polygon into triangles by connecting 0th point to each non-adjacent edge

      const auto verticesCount = entity.template getSubentitiesCount< 0 >();

      const auto v0 = entity.template getSubentityIndex< 0 >( 0 );

      for( LocalIndexType j = 1, k = 2; k < verticesCount; j++, k++ ) {

         const auto v1 = entity.template getSubentityIndex< 0 >( j );

         const auto v2 = entity.template getSubentityIndex< 0 >( k );

         decomposedEntityFunctor( v0, v1, v2 );

      }

   }

};



// Polyhedron

template< typename MeshConfig, EntityDecomposerVersion SubentityDecomposerVersion >

struct EntityDecomposer< MeshConfig, Topologies::Polyhedron, EntityDecomposerVersion::ConnectEdgesToCentroid, SubentityDecomposerVersion >

{

   using Device = typename Devices::Host;

   using CellTopology = typename Topologies::Polyhedron;

   using FaceTopology = typename Topologies::Polygon;

   using MeshEntityType = MeshEntity< MeshConfig, Device, CellTopology >;

   using VertexMeshEntityType = typename MeshEntityType::template SubentityTraits< 0 >::SubentityType;

   using PointType = typename VertexMeshEntityType::PointType;

   using GlobalIndexType = typename MeshConfig::GlobalIndexType;

   using LocalIndexType = typename MeshConfig::LocalIndexType;

   using PointCreationFunctorType = std::function< GlobalIndexType ( const PointType& ) >;

   using DecomposedEntityFunctorType = std::function< void ( GlobalIndexType, GlobalIndexType, GlobalIndexType, GlobalIndexType ) >;

   using SubentityDecomposer = EntityDecomposer< MeshConfig, FaceTopology, SubentityDecomposerVersion >;

   

   static std::pair< GlobalIndexType, GlobalIndexType > getExtraPointsAndEntitiesCount( const MeshEntityType & entity )

   {

      const auto& mesh = entity.getMesh();

      GlobalIndexType extraPointsCount = 1, // there is one new centroid point

                      entitiesCount = 0;

      const auto facesCount = entity.template getSubentitiesCount< 2 >();

      for( LocalIndexType i = 0; i < facesCount; i++ ) {

         const auto face = mesh.template getEntity< 2 >( entity.template getSubentityIndex< 2 >( i ) );



         GlobalIndexType faceExtraPoints, faceEntitiesCount;

         std::tie( faceExtraPoints, faceEntitiesCount ) = SubentityDecomposer::getExtraPointsAndEntitiesCount( face );

         extraPointsCount += faceExtraPoints; // add extra points from decomposition of faces

         entitiesCount += faceEntitiesCount; // there is a new tetrahedron per triangle of a face

      }

      return { extraPointsCount, entitiesCount };

   }



   static void decompose( const MeshEntityType & entity,

                          PointCreationFunctorType pointCreationFunctor,

                          DecomposedEntityFunctorType decomposedEntityFunctor )

   {

      const auto & mesh = entity.getMesh();

      const auto v3 = pointCreationFunctor( getEntityCenter( mesh, entity ) );



      // Lambda for creating tetrahedron from decomposed triangles of faces

      auto decomposedSubentityFunctor = [&] ( GlobalIndexType v0, GlobalIndexType v1, GlobalIndexType v2 ) {

         decomposedEntityFunctor( v0, v1, v2, v3 );

      };



      const auto facesCount = entity.template getSubentitiesCount< 2 >();

      for( LocalIndexType i = 0; i < facesCount; i++ ) {

         const auto face = mesh.template getEntity< 2 >( entity.template getSubentityIndex< 2 >( i ) );

         SubentityDecomposer::decompose( face, pointCreationFunctor, decomposedSubentityFunctor );

      }

   }

};



template< typename MeshConfig, EntityDecomposerVersion SubentityDecomposerVersion >

struct EntityDecomposer< MeshConfig, Topologies::Polyhedron, EntityDecomposerVersion::ConnectEdgesToPoint, SubentityDecomposerVersion >

{

   // https://mathoverflow.net/a/7672

   using Device = typename Devices::Host;

   using CellTopology = typename Topologies::Polyhedron;

   using FaceTopology = typename Topologies::Polygon;

   using MeshEntityType = MeshEntity< MeshConfig, Device, CellTopology >;

   using MeshSubentityType = MeshEntity< MeshConfig, Device, FaceTopology >;

   using VertexMeshEntityType = typename MeshEntityType::template SubentityTraits< 0 >::SubentityType;

   using PointType = typename VertexMeshEntityType::PointType;

   using GlobalIndexType = typename MeshConfig::GlobalIndexType;

   using LocalIndexType = typename MeshConfig::LocalIndexType;

   using RealType = typename MeshConfig::RealType;

   using PointCreationFunctorType = std::function< GlobalIndexType ( const PointType& ) >;

   using DecomposedEntityFunctorType = std::function< void ( GlobalIndexType, GlobalIndexType, GlobalIndexType, GlobalIndexType ) >;

   using SubentityDecomposer = EntityDecomposer< MeshConfig, FaceTopology, SubentityDecomposerVersion >;

   

   static std::pair< GlobalIndexType, GlobalIndexType > getExtraPointsAndEntitiesCount( const MeshEntityType & entity )

   {

      const auto& mesh = entity.getMesh();

      const auto v3 = entity.template getSubentityIndex< 0 >( 0 );

      GlobalIndexType extraPointsCount = 0,

                      entitiesCount = 0;

      const auto facesCount = entity.template getSubentitiesCount< 2 >();

      for( LocalIndexType i = 0; i < facesCount; i++ ) {

         const auto face = mesh.template getEntity< 2 >( entity.template getSubentityIndex< 2 >( i ) );

         if( !faceContainsPoint( face, v3 ) ) { // include only faces, that don't contain point v3

            GlobalIndexType faceExtraPoints, faceEntitiesCount;

            std::tie( faceExtraPoints, faceEntitiesCount ) = SubentityDecomposer::getExtraPointsAndEntitiesCount( face );

            extraPointsCount += faceExtraPoints; // add extra points from decomposition of faces

            entitiesCount += faceEntitiesCount; // there is a new tetrahedron per triangle of a face

         }

      }

      return { extraPointsCount, entitiesCount };

   }



   static void decompose( const MeshEntityType & entity,

                          PointCreationFunctorType pointCreationFunctor,

                          DecomposedEntityFunctorType decomposedEntityFunctor )

   {

      const auto & mesh = entity.getMesh();

      const auto v3 = entity.template getSubentityIndex< 0 >( 0 );



      // Lambda for creating tetrahedron by connecting decomposed triangles of faces to 0th point of polyhedron

      auto decomposedSubentityFunctor = [&] ( GlobalIndexType v0, GlobalIndexType v1, GlobalIndexType v2 ) {

         decomposedEntityFunctor( v0, v1, v2, v3 );

      };



      const auto facesCount = entity.template getSubentitiesCount< 2 >();

      for( LocalIndexType i = 0; i < facesCount; i++ ) {

         const auto face = mesh.template getEntity< 2 >( entity.template getSubentityIndex< 2 >( i ) );

         if( !faceContainsPoint( face, v3 ) ) { // include only faces, that don't contain point v3

            SubentityDecomposer::decompose( face, pointCreationFunctor, decomposedSubentityFunctor );

         }

      }

   }



private:

   static bool faceContainsPoint( const MeshSubentityType & face, const GlobalIndexType point )

   {

      const LocalIndexType pointsCount = face.template getSubentitiesCount< 0 >();

      for( LocalIndexType i = 0; i < pointsCount; i++ ) {

         const auto facePoint = face.template getSubentityIndex< 0 >( i );

         if( point == facePoint )

            return true;

      }

      return false;

   }

};



} // namespace Meshes

} // namespace TNL
 
 No newline at end of file

src/TNL/Meshes/Geometry/PolygonDecomposer.h

deleted100644 → 0
+0 −123
Original line number Diff line number Diff line 
#pragma once



#include <TNL/Meshes/Mesh.h>

#include <TNL/Meshes/MeshEntity.h>

#include <TNL/Meshes/MeshBuilder.h>

#include <TNL/Meshes/Topologies/Triangle.h>

#include <TNL/Meshes/Topologies/Polygon.h>

#include <TNL/Meshes/Geometry/getEntityCenter.h>



namespace TNL {

namespace Meshes {



enum class PolygonDecomposerVersion

{

   ConnectEdgesToCentroid, ConnectEdgesToPoint

};



template< typename MeshConfig, PolygonDecomposerVersion >

struct PolygonDecomposer;



template< typename MeshConfig >

struct PolygonDecomposer< MeshConfig, PolygonDecomposerVersion::ConnectEdgesToCentroid >

{

   using Device = typename Devices::Host;

   using Topology = typename Topologies::Polygon;

   using MeshEntityType = MeshEntity< MeshConfig, Device, Topology >;

   using GlobalIndexType = typename MeshConfig::GlobalIndexType;

   using LocalIndexType = typename MeshConfig::LocalIndexType;



   static GlobalIndexType getExtraPointsCount( const MeshEntityType & entity )

   {

      return 1;

   }



   static GlobalIndexType getEntitiesCount( const MeshEntityType & entity )

   {

      const auto pointsCount = entity.template getSubentitiesCount< 0 >();

      return ( pointsCount == 3 ) ? 1 : pointsCount; // polygon is decomposed only if it has more than 3 vertices

   }



   template< typename MeshType, typename EntitySeedGetterType >

   static void decompose( MeshBuilder< MeshType > & meshBuilder,

                          EntitySeedGetterType entitySeedGetter,

                          GlobalIndexType & pointsCount,

                          GlobalIndexType & entitiesCount,

                          const MeshEntityType & entity )

   {

      const auto verticesCount = entity.template getSubentitiesCount< 0 >();

      if( verticesCount == 3 ) { // polygon is not decomposed as it's already a triangle

         auto & entitySeed = entitySeedGetter( entitiesCount++ );

         entitySeed.setCornersCount( 3 );

         entitySeed.setCornerId( 0, entity.template getSubentityIndex< 0 >( 0 ) );

         entitySeed.setCornerId( 1, entity.template getSubentityIndex< 0 >( 1 ) );

         entitySeed.setCornerId( 2, entity.template getSubentityIndex< 0 >( 2 ) );

      }

      else { // polygon is decomposed as it has got more than 3 vertices

         // add centroid of entity to points

         const auto entityCenter = getEntityCenter( entity.getMesh(), entity );

         const auto entityCenterIdx = pointsCount++;

         meshBuilder.setPoint( entityCenterIdx, entityCenter );

         // decompose polygon into triangles by connecting each edge to the centroid

         for( LocalIndexType j = 0, k = 1; k < verticesCount; j++, k++ ) {

            auto & entitySeed = entitySeedGetter( entitiesCount++ );

            entitySeed.setCornersCount( 3 );

            entitySeed.setCornerId( 0, entity.template getSubentityIndex< 0 >( j ) );

            entitySeed.setCornerId( 1, entity.template getSubentityIndex< 0 >( k ) );

            entitySeed.setCornerId( 2, entityCenterIdx );

         }

         { // wrap around term

            auto & entitySeed = entitySeedGetter( entitiesCount++ );

            entitySeed.setCornersCount( 3 );

            entitySeed.setCornerId( 0, entity.template getSubentityIndex< 0 >( verticesCount - 1 ) );

            entitySeed.setCornerId( 1, entity.template getSubentityIndex< 0 >( 0 ) );

            entitySeed.setCornerId( 2, entityCenterIdx );

         }

      }

   }

};



template< typename MeshConfig >

struct PolygonDecomposer< MeshConfig, PolygonDecomposerVersion::ConnectEdgesToPoint >

{

   using Device = typename Devices::Host;

   using Topology = typename Topologies::Polygon;

   using MeshEntityType = MeshEntity< MeshConfig, Device, Topology >;

   using GlobalIndexType = typename MeshConfig::GlobalIndexType;

   using LocalIndexType = typename MeshConfig::LocalIndexType;



   static GlobalIndexType getExtraPointsCount( const MeshEntityType & entity )

   {

      return 0; // No extra points are added

   }



   static GlobalIndexType getEntitiesCount( const MeshEntityType & entity )

   {

      const auto pointsCount = entity.template getSubentitiesCount< 0 >();

      return pointsCount - 2; // there is a new triangle for every non-adjacent edge to the 0th point

   }



   template< typename MeshType, typename EntitySeedGetterType >

   static void decompose( MeshBuilder< MeshType > & meshBuilder,

                          EntitySeedGetterType entitySeedGetter,

                          GlobalIndexType & pointsCount,

                          GlobalIndexType & entitiesCount,

                          const MeshEntityType & entity )

   {

      // decompose polygon into triangles by connecting 0th point to each non-adjacent edge

      const auto verticesCount = entity.template getSubentitiesCount< 0 >();

      const auto v0 = entity.template getSubentityIndex< 0 >( 0 );

      for( LocalIndexType j = 1, k = 2; k < verticesCount; j++, k++ ) {

         auto & entitySeed = entitySeedGetter( entitiesCount++ );

         const auto v1 = entity.template getSubentityIndex< 0 >( j );

         const auto v2 = entity.template getSubentityIndex< 0 >( k );

         entitySeed.setCornersCount( 3 );

         entitySeed.setCornerId( 0, v0 );

         entitySeed.setCornerId( 1, v1 );

         entitySeed.setCornerId( 2, v2 );

      }

   }

};



} // namespace Meshes

} // namespace TNL
 
 No newline at end of file

src/TNL/Meshes/Geometry/getDecomposedMesh.h

+63 −422

File changed.

Preview size limit exceeded, changes collapsed.

src/TNL/Meshes/Geometry/getPlanarMesh.h

+64 −33

File changed.

Preview size limit exceeded, changes collapsed.

src/UnitTests/Meshes/MeshGeometryTest.h

+30 −32
Original line number Diff line number Diff line
@@ -553,18 +553,18 @@ TEST( MeshGeometryTest, PolygonDecompositionTest ) 


   // Write decomposed mesh using 1st version

   {

      auto triangleMesh = getDecomposedMesh< PolygonDecomposerVersion::ConnectEdgesToCentroid >( mesh );

      auto triangleMesh = getDecomposedMesh< EntityDecomposerVersion::ConnectEdgesToCentroid >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( triangleMesh ) >;

      std::ofstream file( "polygon_decompositionTest_v1.vtk" );

      std::ofstream file( "polygon_decompositionTest_c.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( triangleMesh );

   }



   // Write decomposed mesh using 2nd version

   {

      auto triangleMesh = getDecomposedMesh< PolygonDecomposerVersion::ConnectEdgesToPoint >( mesh );

      auto triangleMesh = getDecomposedMesh< EntityDecomposerVersion::ConnectEdgesToPoint >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( triangleMesh ) >;

      std::ofstream file( "polygon_decompositionTest_v2.vtk" );

      std::ofstream file( "polygon_decompositionTest_p.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( triangleMesh );

   }
@@ -812,45 +812,40 @@ TEST( MeshGeometryTest, PolyhedronDecompositionTest ) 


   // Write decomposed mesh using 1st version

   {

      auto tetrahedronMesh = getDecomposedMesh< GetTetrahedronMeshVersion::V1 >( mesh );

      auto tetrahedronMesh = getDecomposedMesh< EntityDecomposerVersion::ConnectEdgesToCentroid,

                                                EntityDecomposerVersion::ConnectEdgesToCentroid >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( tetrahedronMesh ) >;

      std::ofstream file( "polyhedron_decompositionTest_v1.vtk" );

      std::ofstream file( "polyhedron_decompositionTest_cc.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( tetrahedronMesh );

   }



   // Write decomposed mesh using 2nd version

   {

      auto tetrahedronMesh = getDecomposedMesh< GetTetrahedronMeshVersion::V2 >( mesh );

      auto tetrahedronMesh = getDecomposedMesh< EntityDecomposerVersion::ConnectEdgesToCentroid,

                                                EntityDecomposerVersion::ConnectEdgesToPoint >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( tetrahedronMesh ) >;

      std::ofstream file( "polyhedron_decompositionTest_v2.vtk" );

      std::ofstream file( "polyhedron_decompositionTest_cp.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( tetrahedronMesh );

   }



   // Write decomposed mesh using 3rd version

   {

      auto tetrahedronMesh = getDecomposedMesh< GetTetrahedronMeshVersion::V3 >( mesh );

      auto tetrahedronMesh = getDecomposedMesh< EntityDecomposerVersion::ConnectEdgesToPoint,

                                                EntityDecomposerVersion::ConnectEdgesToCentroid >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( tetrahedronMesh ) >;

      std::ofstream file( "polyhedron_decompositionTest_v3.vtk" );

      std::ofstream file( "polyhedron_decompositionTest_pc.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( tetrahedronMesh );

   }



   // Write decomposed mesh using 4th version

   {

      auto tetrahedronMesh = getDecomposedMesh< GetTetrahedronMeshVersion::V4 >( mesh );

      auto tetrahedronMesh = getDecomposedMesh< EntityDecomposerVersion::ConnectEdgesToPoint,

                                                EntityDecomposerVersion::ConnectEdgesToPoint >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( tetrahedronMesh ) >;

      std::ofstream file( "polyhedron_decompositionTest_v4.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( tetrahedronMesh );

   }



   // Write decomposed mesh using 5th version

   {

      auto tetrahedronMesh = getDecomposedMesh< GetTetrahedronMeshVersion::V5 >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( tetrahedronMesh ) >;

      std::ofstream file( "polyhedron_decompositionTest_v5.vtk" );

      std::ofstream file( "polyhedron_decompositionTest_pp.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( tetrahedronMesh );

   }
@@ -916,18 +911,18 @@ TEST( MeshGeometryTest, Polygon3DGetPlanarMeshTest ) 


   // Write decomposed mesh using 1st version

   {

      auto planarMesh = getPlanarMesh< PolygonDecomposerVersion::ConnectEdgesToCentroid >( mesh );

      auto planarMesh = getPlanarMesh< EntityDecomposerVersion::ConnectEdgesToCentroid >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( planarMesh ) >;

      std::ofstream file( "polygon_planarTest_v1.vtk" );

      std::ofstream file( "polygon_planarTest_c.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( planarMesh );

   }



   // Write decomposed mesh using 2nd version

   {

      auto planarMesh = getPlanarMesh< PolygonDecomposerVersion::ConnectEdgesToPoint >( mesh );

      auto planarMesh = getPlanarMesh< EntityDecomposerVersion::ConnectEdgesToPoint >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( planarMesh ) >;

      std::ofstream file( "polygon_planarTest_v2.vtk" );

      std::ofstream file( "polygon_planarTest_p.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( planarMesh );

   }
@@ -1175,7 +1170,8 @@ TEST( MeshGeometryTest, PolyhedronGetPlanarMeshTest ) 


   // Write original decomposed mesh

   {

      auto tetrahedronMesh = getDecomposedMesh< GetTetrahedronMeshVersion::V1 >( mesh );

      auto tetrahedronMesh = getDecomposedMesh< EntityDecomposerVersion::ConnectEdgesToCentroid,

                                                EntityDecomposerVersion::ConnectEdgesToCentroid >( mesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( tetrahedronMesh ) >;

      std::ofstream file( "polyhedron_planarTest_orig.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );
@@ -1184,20 +1180,22 @@ TEST( MeshGeometryTest, PolyhedronGetPlanarMeshTest ) 


   // Write planar decomposed mesh using 1st version

   {

      auto planarMesh = getPlanarMesh< PolygonDecomposerVersion::ConnectEdgesToCentroid >( mesh );

      auto tetrahedronMesh = getDecomposedMesh< GetTetrahedronMeshVersion::V2 >( planarMesh );

      auto planarMesh = getPlanarMesh< EntityDecomposerVersion::ConnectEdgesToCentroid >( mesh );

      auto tetrahedronMesh = getDecomposedMesh< EntityDecomposerVersion::ConnectEdgesToPoint,

                                                EntityDecomposerVersion::ConnectEdgesToPoint >( planarMesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( tetrahedronMesh ) >;

      std::ofstream file( "polyhedron_planarTest_v1.vtk" );

      std::ofstream file( "polyhedron_planarTest_c.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( tetrahedronMesh );

   }



   // Write planar decomposed mesh using 2nd version

   {

      auto planarMesh = getPlanarMesh< PolygonDecomposerVersion::ConnectEdgesToPoint >( mesh );

      auto tetrahedronMesh = getDecomposedMesh< GetTetrahedronMeshVersion::V1 >( planarMesh );

      auto planarMesh = getPlanarMesh< EntityDecomposerVersion::ConnectEdgesToPoint >( mesh );

      auto tetrahedronMesh = getDecomposedMesh< EntityDecomposerVersion::ConnectEdgesToCentroid,

                                                EntityDecomposerVersion::ConnectEdgesToCentroid >( planarMesh );

      using VTKWriter = Meshes::Writers::VTKWriter< decltype( tetrahedronMesh ) >;

      std::ofstream file( "polyhedron_planarTest_v2.vtk" );

      std::ofstream file( "polyhedron_planarTest_p.vtk" );

      VTKWriter writer( file, VTK::FileFormat::ascii );

      writer.template writeEntities( tetrahedronMesh );

   }