Meshes/Geometry: added function getPlanarMesh + refactoring

- added function getPlanarMesh for correcting non-planar polygons in 3D polygon and polyhedral meshes
- added unit tests for getPlanarMesh function
- refactored getDecomposedMesh specialization for polygons by creating PolygonDecomposer struct to reuse decomposing logic in getPlanarMesh function
- other minor refactoring in Meshes/MeshGeometryTest.h

src/TNL/Meshes/Geometry/PolygonDecomposer.h

+#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

 #pragma once
 
-#include <TNL/Cuda/CudaCallable.h>
 #include <TNL/Meshes/Mesh.h>
 #include <TNL/Meshes/MeshEntity.h>
-#include <TNL/Meshes/Topologies/Vertex.h>
-#include <TNL/Meshes/Topologies/Edge.h>
+#include <TNL/Meshes/MeshBuilder.h>
 #include <TNL/Meshes/Topologies/Triangle.h>
 #include <TNL/Meshes/Topologies/Tetrahedron.h>
 #include <TNL/Meshes/Topologies/Polygon.h>
 #include <TNL/Meshes/Topologies/Polyhedron.h>
-#include <TNL/Meshes/MeshBuilder.h>
 #include <TNL/Meshes/Geometry/getEntityCenter.h>
 #include <TNL/Meshes/Geometry/getEntityMeasure.h>
+#include <TNL/Meshes/Geometry/PolygonDecomposer.h>
 
 namespace TNL {
 namespace Meshes {
 
-
 // Polygon Mesh
 template< typename ParentConfig >
 struct TriangleConfig: public ParentConfig
@@ -24,155 +21,59 @@ struct TriangleConfig: public ParentConfig
    using CellTopology = Topologies::Triangle;
 };
 
-enum class GetTriangleMeshVersion
-{ 
-   V1, V2 
-};
-
-/**
- * Triangles are made by connecting each edge to the centroid of cell.
- */
-template< typename MeshConfig,
+template< PolygonDecomposerVersion version,
+          typename MeshConfig,
           std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Polygon >::value, bool > = true >
 auto
-getTriangleMesh_v1( const Mesh< MeshConfig, Devices::Host > & inMesh )
+getDecomposedMesh( const Mesh< MeshConfig, Devices::Host > & inMesh )
 {
    using TriangleMeshConfig = TriangleConfig< MeshConfig >;
    using TriangleMesh = Mesh< TriangleMeshConfig, Devices::Host >;
+   using MeshBuilder = MeshBuilder< TriangleMesh >;
    using GlobalIndexType = typename TriangleMesh::GlobalIndexType;
    using LocalIndexType = typename TriangleMesh::LocalIndexType;
-
+   using PolygonDecomposer = PolygonDecomposer< MeshConfig, version >;
+   
    TriangleMesh outMesh;
-   MeshBuilder< TriangleMesh > meshBuilder;
+   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
+   // 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 );
-      const auto verticesCount = cell.template getSubentitiesCount< 0 >();
-      if( verticesCount == 3 ) { // cell is not decomposed as it's already a triangle
-         outCellsCount++; // cell is just copied
-      }
-      else { // cell is decomposed as it has got more than 3 vertices
-         outPointsCount++; // each decomposed cell has 1 new center point
-         outCellsCount += verticesCount; // cell is decomposed into verticesCount number of triangles
-      }
+      outPointsCount += PolygonDecomposer::getExtraPointsCount( cell );
+      outCellsCount += PolygonDecomposer::getEntitiesCount( cell );
    }
 
    meshBuilder.setPointsCount( outPointsCount );
    meshBuilder.setCellsCount( outCellsCount );
 
-   // copy the points from inMesh to outMesh
-   GlobalIndexType pointSetIdx = 0;
-   for( ; pointSetIdx < inPointsCount; pointSetIdx++ ) {
-      meshBuilder.setPoint( pointSetIdx, inMesh.getPoint( pointSetIdx ) );
+   // Copy the points from inMesh to outMesh
+   GlobalIndexType setPointsCount = 0;
+   for( ; setPointsCount < inPointsCount; setPointsCount++ ) {
+      meshBuilder.setPoint( setPointsCount, inMesh.getPoint( setPointsCount ) );
    }
 
-   for( GlobalIndexType i = 0, cellSeedIdx = 0; i < inCellsCount; i++ ) {
+   // Decompose each cell into triangles
+   auto cellSeedGetter = [&] ( const GlobalIndexType i ) -> auto& { return meshBuilder.getCellSeed( i ); };
+   for( GlobalIndexType i = 0, setCellsCount = 0; i < inCellsCount; i++ ) {
       const auto cell = inMesh.template getEntity< 2 >( i );
-      const auto verticesCount = cell.template getSubentitiesCount< 0 >();
-      if( verticesCount == 3 ) { // cell is not decomposed as it's already a triangle
-         // copy cell
-         auto & cellSeed = meshBuilder.getCellSeed( cellSeedIdx++ );
-         cellSeed.setCornerId( 0, cell.template getSubentityIndex< 0 >( 0 ) );
-         cellSeed.setCornerId( 1, cell.template getSubentityIndex< 0 >( 1 ) );
-         cellSeed.setCornerId( 2, cell.template getSubentityIndex< 0 >( 2 ) );
-      }
-      else { // cell is decomposed as it has got more than 3 vertices
-         // add centroid of cell to outMesh
-         const auto cellCenter = getEntityCenter( inMesh, cell );
-         const auto cellCenterIdx = pointSetIdx++;
-         meshBuilder.setPoint( cellCenterIdx, cellCenter );
-         // decompose cell into triangles by connecting each edge to the centroid
-         for( LocalIndexType j = 0, k = 1; k < verticesCount; j++, k++ ) {
-            auto & cellSeed = meshBuilder.getCellSeed( cellSeedIdx++ );
-            cellSeed.setCornerId( 0, cell.template getSubentityIndex< 0 >( j ) );
-            cellSeed.setCornerId( 1, cell.template getSubentityIndex< 0 >( k ) );
-            cellSeed.setCornerId( 2, cellCenterIdx );
-         }
-         { // wrap around term
-            auto & cellSeed = meshBuilder.getCellSeed( cellSeedIdx++ );
-            cellSeed.setCornerId( 0, cell.template getSubentityIndex< 0 >( verticesCount - 1 ) );
-            cellSeed.setCornerId( 1, cell.template getSubentityIndex< 0 >( 0 ) );
-            cellSeed.setCornerId( 2, cellCenterIdx );
-         }
-      }
+      PolygonDecomposer::decompose( meshBuilder,
+                                    cellSeedGetter, 
+                                    setPointsCount,
+                                    setCellsCount,
+                                    cell );
    }
    
    meshBuilder.build( outMesh );
    return outMesh;
 }
 
-/**
- * Triangles are made by choosing the 0th point of cell and connecting each non-adjacent edge to it.
- */
-template< typename MeshConfig,
-          std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Polygon >::value, bool > = true >
-auto
-getTriangleMesh_v2( const Mesh< MeshConfig, Devices::Host > & inMesh )
-{
-   using TriangleMeshConfig = TriangleConfig< MeshConfig >;
-   using TriangleMesh = Mesh< TriangleMeshConfig, Devices::Host >;
-   using GlobalIndexType = typename TriangleMesh::GlobalIndexType;
-   using LocalIndexType = typename TriangleMesh::LocalIndexType;
-
-   TriangleMesh outMesh;
-   MeshBuilder< TriangleMesh > meshBuilder;
-   
-   const GlobalIndexType inPointsCount = inMesh.template getEntitiesCount< 0 >();
-   const GlobalIndexType inCellsCount = inMesh.template getEntitiesCount< 2 >();
-
-   // outMesh keeps all the points of inMesh
-   meshBuilder.setPointsCount( inPointsCount );
-   for( GlobalIndexType i = 0; i < inPointsCount; i++ ) {
-      meshBuilder.setPoint( i, inMesh.getPoint( i ) );
-   }
-   
-   // find the number of cells in the outMesh
-   GlobalIndexType outCellsCount = 0;
-   for( GlobalIndexType i = 0; i < inCellsCount; i++ ) {
-      const auto cell = inMesh.template getEntity< 2 >( i );
-      const auto verticesCount = cell.template getSubentitiesCount< 0 >();
-      outCellsCount += verticesCount - 2;
-   }
-   meshBuilder.setCellsCount( outCellsCount );
-
-   for( GlobalIndexType i = 0, cellSeedIdx = 0; i < inCellsCount; i++ ) {
-      const auto cell = inMesh.template getEntity< 2 >( i );
-      const auto verticesCount = cell.template getSubentitiesCount< 0 >();
-      const auto v0 = cell.template getSubentityIndex< 0 >( 0 );
-      for( LocalIndexType j = 1, k = 2; k < verticesCount; j++, k++ ) {
-         auto & cellSeed = meshBuilder.getCellSeed( cellSeedIdx++ );
-         const auto v1 = cell.template getSubentityIndex< 0 >( j );
-         const auto v2 = cell.template getSubentityIndex< 0 >( k );
-         cellSeed.setCornerId( 0, v0 );
-         cellSeed.setCornerId( 1, v1 );
-         cellSeed.setCornerId( 2, v2 );
-      }
-   }
-   
-   meshBuilder.build( outMesh );
-   return outMesh;
-}
-
-template< GetTriangleMeshVersion version,
-          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 )
-{
-   switch( version )
-   {
-      case GetTriangleMeshVersion::V1: return getTriangleMesh_v1( inMesh );
-      case GetTriangleMeshVersion::V2: return getTriangleMesh_v2( inMesh );
-   }
-}
-
 // Polyhedral Mesh
 template< typename ParentConfig >
 struct TetrahedronConfig: public ParentConfig

src/TNL/Meshes/Geometry/getPlanarMesh.h

+#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/Tetrahedron.h>
+#include <TNL/Meshes/Topologies/Polygon.h>
+#include <TNL/Meshes/Topologies/Polyhedron.h>
+#include <TNL/Meshes/Geometry/isPlanar.h>
+#include <TNL/Meshes/Geometry/PolygonDecomposer.h>
+
+namespace TNL {
+namespace Meshes {
+
+// 3D Polygon Mesh
+template< PolygonDecomposerVersion version,
+          typename MeshConfig,
+          std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Polygon >::value, bool > = true,
+          std::enable_if_t< MeshConfig::spaceDimension == 3, bool > = true >
+Mesh< MeshConfig, Devices::Host >
+getPlanarMesh( const Mesh< MeshConfig, Devices::Host > & inMesh )
+{
+   using PolygonMesh = Mesh< MeshConfig, Devices::Host >;
+   using MeshBuilder = MeshBuilder< PolygonMesh >;
+   using GlobalIndexType = typename PolygonMesh::GlobalIndexType;
+   using LocalIndexType = typename PolygonMesh::LocalIndexType;
+   using RealType = typename PolygonMesh::RealType;
+   using PolygonDecomposer = PolygonDecomposer< MeshConfig, version >;
+
+   constexpr RealType precision{ 1e-6 };
+
+   PolygonMesh 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 );
+      if( isPlanar( inMesh, cell, precision ) ) { // Cell is not decomposed
+         outCellsCount++;
+      }
+      else { // Cell is decomposed
+         outPointsCount += PolygonDecomposer::getExtraPointsCount( cell );
+         outCellsCount += PolygonDecomposer::getEntitiesCount( cell );
+      }
+   }
+
+   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 cellSeedGetter = [&] ( const GlobalIndexType i ) -> auto& { return meshBuilder.getCellSeed( i ); };
+   for( GlobalIndexType i = 0, setCellsCount = 0; i < inCellsCount; i++ ) {
+      const auto cell = inMesh.template getEntity< 2 >( i );
+      if( isPlanar( inMesh, cell, precision ) ) { // Copy planar cells
+         auto & cellSeed = meshBuilder.getCellSeed( setCellsCount++ );
+         const auto verticesCount = cell.template getSubentitiesCount< 0 >();
+         cellSeed.setCornersCount( verticesCount );
+         for( LocalIndexType j = 0; j < verticesCount; j++ ) {
+            cellSeed.setCornerId( j, cell.template getSubentityIndex< 0 >( j ) );
+         }
+      }
+      else { // decompose non-planar cells
+         PolygonDecomposer::decompose( meshBuilder,
+                                       cellSeedGetter,
+                                       setPointsCount,
+                                       setCellsCount,
+                                       cell );
+      }
+   }
+
+   meshBuilder.build( outMesh );
+   return outMesh;
+}
+
+// Polyhedral Mesh
+template< PolygonDecomposerVersion version,
+          typename MeshConfig,
+          std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Polyhedron >::value, bool > = true >
+Mesh< MeshConfig, Devices::Host >
+getPlanarMesh( const Mesh< MeshConfig, Devices::Host > & inMesh )
+{
+   using PolyhedronMesh = Mesh< MeshConfig, Devices::Host >;
+   using GlobalIndexType = typename PolyhedronMesh::GlobalIndexType;
+   using LocalIndexType = typename PolyhedronMesh::LocalIndexType;
+   using RealType = typename PolyhedronMesh::RealType;
+   using FaceMapArray = Containers::Array< std::pair< GlobalIndexType, GlobalIndexType >, Devices::Host, GlobalIndexType >;
+   using PolygonDecomposer = PolygonDecomposer< MeshConfig, version >;
+
+   constexpr RealType precision{ 1e-6 };
+
+   PolyhedronMesh outMesh;
+   MeshBuilder< PolyhedronMesh > meshBuilder;
+
+   const GlobalIndexType inPointsCount = inMesh.template getEntitiesCount< 0 >();
+   const GlobalIndexType inFacesCount = inMesh.template getEntitiesCount< 2 >();
+   const GlobalIndexType inCellsCount = inMesh.template getEntitiesCount< 3 >();
+
+   FaceMapArray faceMap( inFacesCount );
+
+   // find the number of points and faces in the outMesh
+   GlobalIndexType outPointsCount = inPointsCount;
+   GlobalIndexType outFacesCount = 0;
+   for( GlobalIndexType i = 0; i < inFacesCount; i++ ) {
+      const auto face = inMesh.template getEntity< 2 >( i );
+      if( isPlanar( inMesh, face, precision ) ) { //
+         const auto startFaceIdx = outFacesCount;
+         const auto endFaceIdx = ++outFacesCount;
+         faceMap[ i ] = { startFaceIdx, endFaceIdx };
+      }
+      else {
+         outPointsCount += PolygonDecomposer::getExtraPointsCount( face );
+         const auto startFaceIdx = outFacesCount;
+         outFacesCount += PolygonDecomposer::getEntitiesCount( face );
+         const auto endFaceIdx = outFacesCount;
+         faceMap[ i ] = { startFaceIdx, endFaceIdx };
+      }
+   }
+
+   meshBuilder.setPointsCount( outPointsCount );
+   meshBuilder.setFacesCount( outFacesCount );
+   meshBuilder.setCellsCount( inCellsCount ); // The number of cells stays the same
+
+   // copy the points from inMesh to outMesh
+   GlobalIndexType setPointsCount = 0;
+   for( ; setPointsCount < inPointsCount; setPointsCount++ ) {
+      meshBuilder.setPoint( setPointsCount, inMesh.getPoint( setPointsCount ) );
+   }
+
+   for( GlobalIndexType i = 0; i < inCellsCount; i++ ) {
+      const auto cell = inMesh.template getEntity< 3 >( i );
+      const LocalIndexType cellFacesCount = cell.template getSubentitiesCount< 2 >();
+
+      // Find the number of faces in the output cell
+      LocalIndexType cellSeedFacesCount = 0;
+      for( LocalIndexType j = 0; j < cellFacesCount; j++ ) {
+         const GlobalIndexType faceIdx = cell.template getSubentityIndex< 2 >( j );
+         const auto & faceMapping = faceMap[ faceIdx ];
+         cellSeedFacesCount += faceMapping.second - faceMapping.first;
+      }
+
+      // Set up the cell seed
+      auto & cellSeed = meshBuilder.getCellSeed( i );
+      cellSeed.setCornersCount( cellSeedFacesCount );
+      for( LocalIndexType j = 0, faceSetIdx = 0; j < cellFacesCount; j++ ) {
+         const GlobalIndexType faceIdx = cell.template getSubentityIndex< 2 >( j );
+         const auto & faceMapping = faceMap[ faceIdx ];
+         for( GlobalIndexType k = faceMapping.first; k < faceMapping.second; k++ ) {
+            cellSeed.setCornerId( faceSetIdx++, k );
+         }
+      }
+   }
+
+   auto faceSeedGetter = [&] ( const GlobalIndexType i ) -> auto& { return meshBuilder.getFaceSeed( i ); };
+   for( GlobalIndexType i = 0, setFacesCount = 0; i < inFacesCount; i++ ) {
+      const auto face = inMesh.template getEntity< 2 >( i );
+      const auto verticesCount = face.template getSubentitiesCount< 0 >();
+      const auto & faceMapping = faceMap[ i ];
+      const bool isPlanarRes = ( faceMapping.second - faceMapping.first ) == 1; // face was planar if face maps only onto 1 face
+      if( isPlanarRes ) { // Copy planar faces
+         auto & faceSeed = meshBuilder.getFaceSeed( setFacesCount++ );
+         faceSeed.setCornersCount( verticesCount );
+         for( LocalIndexType j = 0; j < verticesCount; j++ ) {
+            faceSeed.setCornerId( j, face.template getSubentityIndex< 0 >( j ) );
+         }
+      }
+      else { // decompose non-planar cells
+         PolygonDecomposer::decompose( meshBuilder,
+                                       faceSeedGetter,
+                                       setPointsCount,
+                                       setFacesCount,
+                                       face );
+      }
+   }
+
+   faceMap.reset();
+   meshBuilder.build( outMesh );
+   return outMesh;
+}
+
+} // namespace Meshes
+} // namespace TNL

src/UnitTests/Meshes/CMakeLists.txt

@@ -23,7 +23,7 @@ if( ${BUILD_CUDA} AND ${CUDA_VERSION_MAJOR} GREATER_EQUAL 9 )
    target_link_options( MeshOrderingTest PRIVATE ${TESTS_LINKER_FLAGS} )
 
    CUDA_ADD_EXECUTABLE( MeshGeometryTest MeshGeometryTest.cu
-                        OPTIONS ${CXX_TESTS_FLAGS} )
+                        OPTIONS ${CUDA_TESTS_FLAGS} )
    TARGET_LINK_LIBRARIES( MeshGeometryTest ${TESTS_LIBRARIES} )
    target_link_options( MeshGeometryTest PRIVATE ${TESTS_LINKER_FLAGS} )
 else()