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

Refactored getDecomposedMesh function to use ParallelFor

parent c9b50715

src/Benchmarks/Mesh/MeshBenchmarks.h

+2 −2
Original line number Diff line number Diff line
@@ -341,7 +341,7 @@ struct MeshBenchmarks 
   static void benchmark_decomposition( Benchmark & benchmark, const Config::ParameterContainer & parameters, const M & mesh_src )

   {

      auto benchmark_func = [&] () {

         const auto decomposedMesh = getDecomposedMesh< DecomposerVersion >( mesh_src );

         auto meshBuilder = decomposeMesh< DecomposerVersion >( mesh_src );

      };



      benchmark.time< Devices::Host >( "CPU",
@@ -356,7 +356,7 @@ struct MeshBenchmarks 
   static void benchmark_decomposition( Benchmark & benchmark, const Config::ParameterContainer & parameters, const M & mesh_src )

   {

      auto benchmark_func = [&] () {

         const auto decomposedMesh = getDecomposedMesh< DecomposerVersion, SubDecomposerVersion >( mesh_src );

         auto meshBuilder = decomposeMesh< DecomposerVersion, SubDecomposerVersion >( mesh_src );

      };



      benchmark.time< Devices::Host >( "CPU",

src/TNL/Assert.h

+7 −0
Original line number Diff line number Diff line
@@ -217,6 +217,13 @@ fatalFailure() 
#endif

}



template< typename T >

::std::stringstream& operator<<( ::std::stringstream& ss, const std::pair< T, T >& pair )

{

   ss << '(' << pair.first << ',' <<  pair.second << ')';

   return ss;

}



template< typename T >

struct Formatter

{

src/TNL/Meshes/Geometry/EntityDecomposer.h

+38 −38
Original line number Diff line number Diff line
@@ -34,8 +34,6 @@ struct EntityDecomposer< MeshConfig, Topologies::Polygon, EntityDecomposerVersio 
   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 )

   {
@@ -45,29 +43,31 @@ struct EntityDecomposer< MeshConfig, Topologies::Polygon, EntityDecomposerVersio 
      return { 1, pointsCount }; // 1 extra centroid point and decomposition creates pointsCount triangles

   }



   template< typename AddPointFunctor,

             typename AddCellFunctor >

   static void decompose( const MeshEntityType& entity,

                          PointCreationFunctorType pointCreationFunctor,

                          DecomposedEntityFunctorType decomposedEntityFunctor )

                          AddPointFunctor&& addPoint,

                          AddCellFunctor&& addCell )

   {

      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 );

         addCell( v0, v1, v2 );

      }

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

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

         const auto v0 = addPoint( 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 );

            addCell( 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 );

            addCell( v0, v1, v2 );

         }

      }

   }
@@ -83,8 +83,6 @@ struct EntityDecomposer< MeshConfig, Topologies::Polygon, EntityDecomposerVersio 
   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 )

   {
@@ -92,9 +90,11 @@ struct EntityDecomposer< MeshConfig, Topologies::Polygon, EntityDecomposerVersio 
      return { 0, pointsCount - 2 }; // no extra points and there is a triangle for every non-adjacent edge to the 0th point (pointsCount - 2)

   }



   template< typename AddPointFunctor,

             typename AddCellFunctor >

   static void decompose( const MeshEntityType& entity,

                          PointCreationFunctorType pointCreationFunctor,

                          DecomposedEntityFunctorType decomposedEntityFunctor )

                          AddPointFunctor&& addPoint,

                          AddCellFunctor&& addCell )

   {

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

      const auto verticesCount = entity.template getSubentitiesCount< 0 >();
@@ -102,7 +102,7 @@ struct EntityDecomposer< MeshConfig, Topologies::Polygon, EntityDecomposerVersio 
      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 );

         addCell( v0, v1, v2 );

      }

   }

};
@@ -119,8 +119,6 @@ struct EntityDecomposer< MeshConfig, Topologies::Polyhedron, EntityDecomposerVer 
   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 )
@@ -140,22 +138,24 @@ struct EntityDecomposer< MeshConfig, Topologies::Polyhedron, EntityDecomposerVer 
      return { extraPointsCount, entitiesCount };

   }



   template< typename AddPointFunctor,

             typename AddCellFunctor >

   static void decompose( const MeshEntityType & entity,

                          PointCreationFunctorType pointCreationFunctor,

                          DecomposedEntityFunctorType decomposedEntityFunctor )

                          AddPointFunctor&& addPoint,

                          AddCellFunctor&& addCell )

   {

      const auto& mesh = entity.getMesh();

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

      const auto v3 = addPoint( 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 );

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

         addCell( 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 );

         SubentityDecomposer::decompose( face, std::forward< AddPointFunctor >( addPoint ), addFace );

      }

   }

};
@@ -174,8 +174,6 @@ struct EntityDecomposer< MeshConfig, Topologies::Polyhedron, EntityDecomposerVer 
   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 )
@@ -197,23 +195,25 @@ struct EntityDecomposer< MeshConfig, Topologies::Polyhedron, EntityDecomposerVer 
      return { extraPointsCount, entitiesCount };

   }



   template< typename AddPointFunctor,

             typename AddCellFunctor >

   static void decompose( const MeshEntityType& entity,

                          PointCreationFunctorType pointCreationFunctor,

                          DecomposedEntityFunctorType decomposedEntityFunctor )

                          AddPointFunctor&& addPoint,

                          AddCellFunctor&& addCell )

   {

      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 );

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

         addCell( 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 );

            SubentityDecomposer::decompose( face, std::forward< AddPointFunctor >( addPoint ), addFace );

         }

      }

   }

src/TNL/Meshes/Geometry/getDecomposedMesh.h

+144 −85
Original line number Diff line number Diff line
@@ -7,9 +7,9 @@ 
#include <TNL/Meshes/Topologies/Tetrahedron.h>

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

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

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

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

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

#include <TNL/Algorithms/ParallelFor.h>

#include <TNL/Algorithms/scan.h>



namespace TNL {

namespace Meshes {
@@ -21,67 +21,95 @@ struct TriangleConfig: public ParentConfig 
   using CellTopology = Topologies::Triangle;

};



template< EntityDecomposerVersion EntityDecomposerVersion,

template< EntityDecomposerVersion DecomposerVersion,

          typename MeshConfig,

          std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Polygon >::value, bool > = true >

auto

getDecomposedMesh( const Mesh< MeshConfig, Devices::Host > & inMesh )

auto // returns MeshBuilder

decomposeMesh( const Mesh< MeshConfig, Devices::Host >& inMesh )

{

   using namespace TNL;

   using namespace TNL::Containers;

   using namespace TNL::Algorithms;



   using TriangleMeshConfig = TriangleConfig< MeshConfig >;

   using TriangleMesh = Mesh< TriangleMeshConfig, Devices::Host >;

   using MeshBuilder = MeshBuilder< TriangleMesh >;

   using PointType = typename TriangleMesh::PointType;

   using GlobalIndexType = typename TriangleMesh::GlobalIndexType;

   using LocalIndexType = typename TriangleMesh::LocalIndexType;

   using EntityDecomposer = EntityDecomposer< MeshConfig, Topologies::Polygon, EntityDecomposerVersion >;

   using PointType = typename TriangleMesh::PointType;

   using EntityDecomposer = EntityDecomposer< MeshConfig, Topologies::Polygon, DecomposerVersion >;

   

   TriangleMesh outMesh;

   MeshBuilder meshBuilder;



   const GlobalIndexType inPointsCount = inMesh.template getEntitiesCount< 0 >();

   const GlobalIndexType inCellsCount = inMesh.template getEntitiesCount< 2 >();



   // Find the number of points and cells in the outMesh

   GlobalIndexType outPointsCount = inPointsCount;

   GlobalIndexType outCellsCount = 0;

   for( GlobalIndexType i = 0; i < inCellsCount; i++ ) {

      const auto cell = inMesh.template getEntity< 2 >( i );

      GlobalIndexType extraPointsCount, entitiesCount;

      std::tie( extraPointsCount, entitiesCount ) = EntityDecomposer::getExtraPointsAndEntitiesCount( cell );

      outPointsCount += extraPointsCount;

      outCellsCount += entitiesCount;

   }

   const GlobalIndexType inCellsCount = inMesh.template getEntitiesCount< TriangleMesh::getMeshDimension() >();



   // Find the number of output points and cells as well as

   // starting indeces at which every cell will start writing new decomposed points and cells

   using IndexPair = std::pair< GlobalIndexType, GlobalIndexType >;

   Array< IndexPair, Devices::Host > indeces( inCellsCount );

   auto setCounts = [&] ( GlobalIndexType i ) {

      const auto cell = inMesh.template getEntity< TriangleMesh::getMeshDimension() >( i );

      indeces[ i ] = EntityDecomposer::getExtraPointsAndEntitiesCount( cell );

   };

   ParallelFor< Devices::Host >::exec( 0, inCellsCount, setCounts );

   const auto lastCounts = indeces[ indeces.getSize() - 1 ];

   auto reduction = [] ( const IndexPair& a, const IndexPair& b ) -> IndexPair {

      return { a.first + b.first, a.second + b.second };

   };

   inplaceExclusiveScan( indeces, 0, indeces.getSize(), reduction, std::make_pair( 0, 0 ) );

   const auto lastIndexPair = indeces[ indeces.getSize() - 1 ];

   const GlobalIndexType outPointsCount = inPointsCount + lastIndexPair.first + lastCounts.first;

   const GlobalIndexType outCellsCount = lastIndexPair.second + lastCounts.second;

   meshBuilder.setPointsCount( outPointsCount );

   meshBuilder.setCellsCount( outCellsCount );



   // Copy the points from inMesh to outMesh

   GlobalIndexType setPointsCount = 0;

   for( ; setPointsCount < inPointsCount; setPointsCount++ ) {

      meshBuilder.setPoint( setPointsCount, inMesh.getPoint( setPointsCount ) );

   }

   auto copyPoint = [&] ( GlobalIndexType i ) mutable {

      meshBuilder.setPoint( i, inMesh.getPoint( i ) );

   };

   ParallelFor< Devices::Host >::exec( 0, inPointsCount, copyPoint );



   // Lambda for creating new points

   auto createPointFunc = [&] ( const PointType & point ) {

      const auto pointIdx = setPointsCount++;

   // Decompose each cell

   auto decomposeCell = [&] ( GlobalIndexType i ) mutable {

      const auto cell = inMesh.template getEntity< TriangleMesh::getMeshDimension() >( i );

      const auto indexPair = indeces[ i ];



      // Lambda for adding new points

      GlobalIndexType setPointIndex = inPointsCount + indexPair.first;

      auto addPoint = [&] ( const PointType& point ) {

         const auto pointIdx = setPointIndex++;

         meshBuilder.setPoint( pointIdx, point );

         return pointIdx;

      };



   // Lambda for setting decomposed triangle in meshBuilder

   GlobalIndexType setCellsCount = 0;

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

      auto entitySeed = meshBuilder.getCellSeed( setCellsCount++ );

      // Lambda for adding new cells

      GlobalIndexType setCellIndex = indexPair.second;

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

         auto entitySeed = meshBuilder.getCellSeed( setCellIndex++ );

         entitySeed.setCornerId( 0, v0 );

         entitySeed.setCornerId( 1, v1 );

         entitySeed.setCornerId( 2, v2 );

      };



   // Decompose each cell

   for( GlobalIndexType i = 0; i < inCellsCount; i++ ) {

      const auto cell = inMesh.template getEntity< 2 >( i );

      EntityDecomposer::decompose( cell, createPointFunc, setDecomposedCellFunc );

      EntityDecomposer::decompose( cell, addPoint, addCell );

   };

   ParallelFor< Devices::Host >::exec( 0, inCellsCount, decomposeCell );



   return meshBuilder;

}



template< EntityDecomposerVersion DecomposerVersion,

          typename MeshConfig,

          std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Polygon >::value, bool > = true >

auto // returns Mesh

getDecomposedMesh( const Mesh< MeshConfig, Devices::Host >& inMesh )

{

   using TriangleMeshConfig = TriangleConfig< MeshConfig >;

   using TriangleMesh = Mesh< TriangleMeshConfig, Devices::Host >;

   

   TriangleMesh outMesh;

   auto meshBuilder = decomposeMesh< DecomposerVersion >( inMesh );

   meshBuilder.build( outMesh );

   return outMesh;

}
@@ -93,68 +121,99 @@ struct TetrahedronConfig: public ParentConfig 
   using CellTopology = Topologies::Tetrahedron;

};



template< EntityDecomposerVersion EntityDecomposerVersion_,

          EntityDecomposerVersion SubentityDecomposerVersion,

template< EntityDecomposerVersion DecomposerVersion,

          EntityDecomposerVersion SubdecomposerVersion,

          typename MeshConfig,

          std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Polyhedron >::value, bool > = true >

auto

getDecomposedMesh( const Mesh< MeshConfig, Devices::Host > & inMesh )

auto // returns MeshBuilder

decomposeMesh( const Mesh< MeshConfig, Devices::Host > & inMesh )

{

   using namespace TNL;

   using namespace TNL::Containers;

   using namespace TNL::Algorithms;



   using TetrahedronMeshConfig = TetrahedronConfig< MeshConfig >;

   using TetrahedronMesh = Mesh< TetrahedronMeshConfig, Devices::Host >;

   using MeshBuilder = MeshBuilder< TetrahedronMesh >;

   using GlobalIndexType = typename TetrahedronMesh::GlobalIndexType;

   using LocalIndexType = typename TetrahedronMesh::LocalIndexType;

   using PointType = typename TetrahedronMesh::PointType;

   using EntityDecomposer = EntityDecomposer< MeshConfig, Topologies::Polyhedron, EntityDecomposerVersion_, SubentityDecomposerVersion >;

   using EntityDecomposer = EntityDecomposer< MeshConfig, Topologies::Polyhedron, DecomposerVersion, SubdecomposerVersion >;

   

   TetrahedronMesh outMesh;

   MeshBuilder< TetrahedronMesh > meshBuilder;

   MeshBuilder meshBuilder;



   const GlobalIndexType inPointsCount = inMesh.template getEntitiesCount< 0 >();

   const GlobalIndexType inCellsCount = inMesh.template getEntitiesCount< 3 >();



   // Find the number of points and cells in the outMesh

   GlobalIndexType outPointsCount = inPointsCount;

   GlobalIndexType outCellsCount = 0;

   for( GlobalIndexType i = 0; i < inCellsCount; i++ ) {

      const auto cell = inMesh.template getEntity< 3 >( i );

      GlobalIndexType extraPointsCount, entitiesCount;

      std::tie( extraPointsCount, entitiesCount ) = EntityDecomposer::getExtraPointsAndEntitiesCount( cell );

      outPointsCount += extraPointsCount;

      outCellsCount += entitiesCount;

   }

   const GlobalIndexType inCellsCount = inMesh.template getEntitiesCount< TetrahedronMesh::getMeshDimension() >();



   using IndexPair = std::pair< GlobalIndexType, GlobalIndexType >;

   Array< IndexPair, Devices::Host > indeces( inCellsCount );



   // Find the number of output points and cells as well as

   // starting indeces at which every cell will start writing new decomposed points and cells

   auto setCounts = [&] ( GlobalIndexType i ) {

      const auto cell = inMesh.template getEntity< TetrahedronMesh::getMeshDimension() >( i );

      indeces[ i ] = EntityDecomposer::getExtraPointsAndEntitiesCount( cell );

   };

   ParallelFor< Devices::Host >::exec( 0, inCellsCount, setCounts );

   const auto lastCounts = indeces[ indeces.getSize() - 1 ];

   auto reduction = [] ( const IndexPair& a, const IndexPair& b ) -> IndexPair {

      return { a.first + b.first, a.second + b.second };

   };

   inplaceExclusiveScan( indeces, 0, indeces.getSize(), reduction, std::make_pair( 0, 0 ) );

   const auto lastIndexPair = indeces[ indeces.getSize() - 1 ];

   const GlobalIndexType outPointsCount = inPointsCount + lastIndexPair.first + lastCounts.first;

   const GlobalIndexType outCellsCount = lastIndexPair.second + lastCounts.second;

   meshBuilder.setPointsCount( outPointsCount );

   meshBuilder.setCellsCount( outCellsCount );



   // Copy the points from inMesh to outMesh

   GlobalIndexType setPointsCount = 0;

   for( ; setPointsCount < inPointsCount; setPointsCount++ ) {

      meshBuilder.setPoint( setPointsCount, inMesh.getPoint( setPointsCount ) );

   }

   auto copyPoint = [&] ( GlobalIndexType i ) mutable {

      meshBuilder.setPoint( i, inMesh.getPoint( i ) );

   };

   ParallelFor< Devices::Host >::exec( 0, inPointsCount, copyPoint );



   // Lambda for creating new points

   auto createPointFunc = [&] ( const PointType & point ) {

      const auto pointIdx = setPointsCount++;

   // Decompose each cell

   auto decomposeCell = [&] ( GlobalIndexType i ) mutable {

      const auto cell = inMesh.template getEntity< TetrahedronMesh::getMeshDimension() >( i );

      const auto indexPair = indeces[ i ];



      // Lambda for adding new points

      GlobalIndexType setPointIndex = inPointsCount + indexPair.first;

      auto addPoint = [&] ( const PointType& point ) {

         const auto pointIdx = setPointIndex++;

         meshBuilder.setPoint( pointIdx, point );

         return pointIdx;

      };



   // Lambda for setting decomposed cells in meshBuilder

   GlobalIndexType setCellsCount = 0;

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

      auto entitySeed = meshBuilder.getCellSeed( setCellsCount++ );

      // Lambda for adding new cells

      GlobalIndexType setCellIndex = indexPair.second;

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

         auto entitySeed = meshBuilder.getCellSeed( setCellIndex++ );

         entitySeed.setCornerId( 0, v0 );

         entitySeed.setCornerId( 1, v1 );

         entitySeed.setCornerId( 2, v2 );

         entitySeed.setCornerId( 3, v3 );

      };



   // Decompose each cell

   for( GlobalIndexType i = 0; i < inCellsCount; i++ ) {

      const auto cell = inMesh.template getEntity< 3 >( i );

      EntityDecomposer::decompose( cell, createPointFunc, setDecomposedCellFunc );

      EntityDecomposer::decompose( cell, addPoint, addCell );

   };

   ParallelFor< Devices::Host >::exec( 0, inCellsCount, decomposeCell );



   return meshBuilder;

}



template< EntityDecomposerVersion DecomposerVersion,

          EntityDecomposerVersion SubDecomposerVersion,

          typename MeshConfig,

          std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Polyhedron >::value, bool > = true >

auto // returns Mesh

getDecomposedMesh( const Mesh< MeshConfig, Devices::Host > & inMesh )

{

   using TetrahedronMeshConfig = TetrahedronConfig< MeshConfig >;

   using TetrahedronMesh = Mesh< TetrahedronMeshConfig, Devices::Host >;

   

   TetrahedronMesh outMesh;

   auto meshBuilder = decomposeMesh< DecomposerVersion, SubDecomposerVersion >( inMesh );

   meshBuilder.build( outMesh );

   return outMesh;

}