diff --git a/src/TNL/Meshes/MeshDetails/initializer/EntityInitializer.h b/src/TNL/Meshes/MeshDetails/initializer/EntityInitializer.h
index 5671ffd779a08459b4537ba750f3b5bba4171933..d03a6e5f422aa8a74369c1b2646b765be7d13625 100644
--- a/src/TNL/Meshes/MeshDetails/initializer/EntityInitializer.h
+++ b/src/TNL/Meshes/MeshDetails/initializer/EntityInitializer.h
@@ -124,6 +124,7 @@ class EntityInitializerLayer< MeshConfig,
    using SubentitySeedsCreatorType  = SubentitySeedsCreator< MeshConfig, SuperentityTopology, SubdimensionTag >;
    using SuperentityMatrixType      = typename MeshTraitsType::SuperentityMatrixType;
    using NeighborCountsArray        = typename MeshTraitsType::NeighborCountsArray;
+   using SeedType                   = EntitySeed< MeshConfig, SubentityTopology >;
 
 public:
    static void initSuperentities( InitializerType& meshInitializer, MeshType& mesh )
@@ -150,13 +151,12 @@ public:
 
       for( GlobalIndexType superentityIndex = 0; superentityIndex < superentitiesCount; superentityIndex++ )
       {
-         auto subentitySeeds = SubentitySeedsCreatorType::create( meshInitializer, mesh, superentityIndex );
-         for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )
-         {
-            const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( subentitySeeds[ i ] );
-            meshInitializer.template setSubentityIndex< SuperdimensionTag::value, SubdimensionTag::value >( superentityIndex, i, subentityIndex );
+         LocalIndexType i = 0;
+         SubentitySeedsCreatorType::iterate( meshInitializer, mesh, superentityIndex, [&] ( SeedType& seed ) {
+            const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( seed );
+            meshInitializer.template setSubentityIndex< SuperdimensionTag::value, SubdimensionTag::value >( superentityIndex, i++, subentityIndex );
             superentitiesCounts[ subentityIndex ]++;
-         }
+         });
       }
 
       // allocate superentities storage
@@ -210,10 +210,13 @@ class EntityInitializerLayer< MeshConfig,
 
    using GlobalIndexType           = typename MeshTraitsType::GlobalIndexType;
    using LocalIndexType            = typename MeshTraitsType::LocalIndexType;
+   using SubentityTraitsType        = typename MeshTraitsType::template EntityTraits< SubdimensionTag::value >;
+   using SubentityTopology          = typename SubentityTraitsType::EntityTopology;
    using SuperentityTraitsType     = typename MeshTraitsType::template EntityTraits< SuperdimensionTag::value >;
    using SuperentityTopology       = typename SuperentityTraitsType::EntityTopology;
    using SubentitySeedsCreatorType = SubentitySeedsCreator< MeshConfig, SuperentityTopology, SubdimensionTag >;
    using NeighborCountsArray       = typename MeshTraitsType::NeighborCountsArray;
+   using SeedType                  = EntitySeed< MeshConfig, SubentityTopology >;
 
 public:
    static void initSuperentities( InitializerType& meshInitializer, MeshType& mesh )
@@ -236,12 +239,11 @@ public:
            superentityIndex < mesh.template getEntitiesCount< SuperdimensionTag::value >();
            superentityIndex++ )
       {
-         auto subentitySeeds = SubentitySeedsCreatorType::create( meshInitializer, mesh, superentityIndex );
-         for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )
-         {
-            const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( subentitySeeds[ i ] );
-            meshInitializer.template setSubentityIndex< SuperdimensionTag::value, SubdimensionTag::value >( superentityIndex, i, subentityIndex );
-         }
+         LocalIndexType i = 0;
+         SubentitySeedsCreatorType::iterate( meshInitializer, mesh, superentityIndex, [&] ( SeedType& seed ) {
+            const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( seed );
+            meshInitializer.template setSubentityIndex< SuperdimensionTag::value, SubdimensionTag::value >( superentityIndex, i++, subentityIndex );
+         });
       }
 
       BaseType::initSuperentities( meshInitializer, mesh );
@@ -282,6 +284,7 @@ class EntityInitializerLayer< MeshConfig,
    using SuperentityTopology        = typename SuperentityTraitsType::EntityTopology;
    using SubentitySeedsCreatorType  = SubentitySeedsCreator< MeshConfig, SuperentityTopology, SubdimensionTag >;
    using SuperentityMatrixType      = typename MeshTraitsType::SuperentityMatrixType;
+   using SeedType                   = EntitySeed< MeshConfig, SubentityTopology >;
 
 public:
    static void initSuperentities( InitializerType& meshInitializer, MeshType& mesh )
@@ -298,12 +301,10 @@ public:
 
       for( GlobalIndexType superentityIndex = 0; superentityIndex < superentitiesCount; superentityIndex++ )
       {
-         auto subentitySeeds = SubentitySeedsCreatorType::create( meshInitializer, mesh, superentityIndex );
-         for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )
-         {
-            const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( subentitySeeds[ i ] );
+         SubentitySeedsCreatorType::iterate( meshInitializer, mesh, superentityIndex, [&] ( SeedType& seed ) {
+            const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( seed );
             superentitiesCounts[ subentityIndex ]++;
-         }
+         });
       }
 
       // allocate superentities storage
@@ -315,13 +316,11 @@ public:
       // initialize superentities storage
       for( GlobalIndexType superentityIndex = 0; superentityIndex < superentitiesCount; superentityIndex++ )
       {
-         auto subentitySeeds = SubentitySeedsCreatorType::create( meshInitializer, mesh, superentityIndex );
-         for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )
-         {
-            const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( subentitySeeds[ i ] );
+         SubentitySeedsCreatorType::iterate( meshInitializer, mesh, superentityIndex, [&] ( SeedType& seed ) {
+            const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( seed );
             auto row = matrix.getRow( subentityIndex );
             row.setElement( superentitiesCounts[ subentityIndex ]++, superentityIndex, true );
-         }
+         });
       }
 
       BaseType::initSuperentities( meshInitializer, mesh );
diff --git a/src/TNL/Meshes/MeshDetails/initializer/Initializer.h b/src/TNL/Meshes/MeshDetails/initializer/Initializer.h
index 5f6134c83e86c2172d160e30b2be29fc390cd5e3..eb8e54564c160644aebab7a4206d27fe579fdc17 100644
--- a/src/TNL/Meshes/MeshDetails/initializer/Initializer.h
+++ b/src/TNL/Meshes/MeshDetails/initializer/Initializer.h
@@ -251,47 +251,20 @@ protected:
    using EntityInitializerType = EntityInitializer< MeshConfig, EntityTopology >;
    using SeedType              = EntitySeed< MeshConfig, EntityTopology >;
    using SeedIndexedSet        = typename MeshTraits< MeshConfig >::template EntityTraits< DimensionTag::value >::SeedIndexedSetType;
-   using SeedSet               = typename MeshTraits< MeshConfig >::template EntityTraits< DimensionTag::value >::SeedSetType;
    using CellSeedArrayType     = typename MeshTraitsType::CellSeedArrayType;
    using EntitySeedArrayType   = Containers::Array< SeedType, typename MeshTraitsType::DeviceType, GlobalIndexType >;
    using NeighborCountsArray   = typename MeshTraitsType::NeighborCountsArray;
-
+   
    public:
 
-      GlobalIndexType getEntitiesCount( InitializerType& initializer, MeshType& mesh )
+      void createSeeds( InitializerType& initializer, MeshType& mesh )
       {
          using SubentitySeedsCreator = SubentitySeedsCreator< MeshConfig, typename MeshTraitsType::CellTopology, DimensionTag >;
-         SeedSet seedSet;
-
-         for( GlobalIndexType i = 0; i < mesh.template getEntitiesCount< MeshType::getMeshDimension() >(); i++ )
-         {
-            auto subentitySeeds = SubentitySeedsCreator::create( initializer, mesh, i );
-            for( LocalIndexType j = 0; j < subentitySeeds.getSize(); j++ )
-               seedSet.insert( subentitySeeds[ j ] );
+         for( GlobalIndexType i = 0; i < mesh.template getEntitiesCount< MeshType::getMeshDimension() >(); i++ ) {
+            SubentitySeedsCreator::iterate( initializer, mesh, i, [&] ( SeedType& seed ) {
+               this->seedsIndexedSet.insert( seed );
+            });
          }
-
-         return seedSet.size();
-      }
-
-      NeighborCountsArray getCapacities( InitializerType& initializer, MeshType& mesh, const GlobalIndexType numberOfEntities )
-      {
-         using SubentitySeedsCreator = SubentitySeedsCreator< MeshConfig, typename MeshTraitsType::CellTopology, DimensionTag >;
-         SeedSet seedSet;
-         NeighborCountsArray capacities( numberOfEntities );
-         GlobalIndexType capSize = 0;
-
-         for( GlobalIndexType i = 0; i < mesh.template getEntitiesCount< MeshType::getMeshDimension() >(); i++ )
-         {
-            auto subentitySeeds = SubentitySeedsCreator::create( initializer, mesh, i );
-            for( LocalIndexType j = 0; j < subentitySeeds.getSize(); j++ )
-            {
-               const auto& subentitySeed = subentitySeeds[ j ];
-               bool inserted = seedSet.insert( subentitySeed ).second;
-               if( inserted )
-                  capacities.setElement( capSize++, subentitySeed.getCornersCount() );
-            }
-         }
-         return capacities;
       }
 
       using BaseType::findEntitySeedIndex;
@@ -303,51 +276,56 @@ protected:
       void initEntities( InitializerType& initializer, MeshType& mesh )
       {
          //std::cout << " Initiating entities with dimension " << DimensionTag::value << " ... " << std::endl;
-         const GlobalIndexType numberOfEntities = getEntitiesCount( initializer, mesh );
+
+         // create seeds and set entities count
+         createSeeds( initializer, mesh );
+         const GlobalIndexType numberOfEntities = this->seedsIndexedSet.size();
          initializer.template setEntitiesCount< DimensionTag::value >( numberOfEntities );
-         NeighborCountsArray capacities = getCapacities( initializer, mesh, numberOfEntities );
+
+         // allocate the subvertex matrix
+         NeighborCountsArray capacities( numberOfEntities );
+         for( auto& pair : this->seedsIndexedSet ) {
+            const auto& seed = pair.first;
+            const auto& entityIndex = pair.second;
+            capacities.setElement( entityIndex, seed.getCornersCount() );
+         }
          EntityInitializerType::initSubvertexMatrix( capacities, initializer );
-         
-         using SubentitySeedsCreator = SubentitySeedsCreator< MeshConfig, typename MeshTraitsType::CellTopology, DimensionTag >;
-         for( GlobalIndexType i = 0; i < mesh.template getEntitiesCount< MeshType::getMeshDimension() >(); i++ )
-         {
-            auto subentitySeeds = SubentitySeedsCreator::create( initializer, mesh, i );
-            for( LocalIndexType j = 0; j < subentitySeeds.getSize(); j++ )
-            {
-               auto& seed = subentitySeeds[ j ];
-               const auto pair = this->seedsIndexedSet.try_insert( seed );
-               const GlobalIndexType& entityIndex = pair.first;
-               if( pair.second ) {
-                  // insertion took place, initialize the entity
-                  EntityInitializerType::initEntity( entityIndex, seed, initializer );
-               }
-            }
+
+         // initialize the entities
+         for( auto& pair : this->seedsIndexedSet ) {
+            const auto& seed = pair.first;
+            const auto& entityIndex = pair.second;
+            EntityInitializerType::initEntity( entityIndex, seed, initializer );
          }
 
+         // initialize links between the entities and all superentities
          EntityInitializerType::initSuperentities( initializer, mesh );
-         this->seedsIndexedSet.clear();
 
+         // deallocate the indexed set and continue with the next dimension
+         this->seedsIndexedSet.clear();
          BaseType::initEntities( initializer, mesh );
       }
 
       void initEntities( InitializerType& initializer, EntitySeedArrayType& faceSeeds, MeshType& mesh )
       {
          //std::cout << " Initiating entities with dimension " << DimensionTag::value << " ... " << std::endl;
+
          initializer.template setEntitiesCount< DimensionTag::value >( faceSeeds.getSize() );
 
+         // allocate the subvertex matrix
          NeighborCountsArray capacities( faceSeeds.getSize() );
-
-         for( LocalIndexType i = 0; i < capacities.getSize(); i++ )
-            capacities[ i ] = faceSeeds[ i ].getCornersCount();
-
+         for( LocalIndexType i = 0; i < capacities.getSize(); i++ ) {
+            capacities.setElement( i, faceSeeds[ i ].getCornersCount() );
+         }
          EntityInitializerType::initSubvertexMatrix( capacities, initializer );
 
-         for( GlobalIndexType i = 0; i < faceSeeds.getSize(); i++ )
-         {
+         // initialize the entities
+         for( GlobalIndexType i = 0; i < faceSeeds.getSize(); i++ ) {
             const auto& seed = faceSeeds[ i ];
             GlobalIndexType entityIndex = this->seedsIndexedSet.insert( seed );
             EntityInitializerType::initEntity( entityIndex, seed, initializer );
          }
+
          faceSeeds.reset();
          EntityInitializerType::initSuperentities( initializer, mesh );
          this->seedsIndexedSet.clear();
diff --git a/src/TNL/Meshes/MeshDetails/initializer/SubentitySeedsCreator.h b/src/TNL/Meshes/MeshDetails/initializer/SubentitySeedsCreator.h
index 57a17297f66b58cf7887bc633f87ec4bdb1e6c6b..f091294727eb477b6f26b0d99e409484f5aa09c6 100644
--- a/src/TNL/Meshes/MeshDetails/initializer/SubentitySeedsCreator.h
+++ b/src/TNL/Meshes/MeshDetails/initializer/SubentitySeedsCreator.h
@@ -37,9 +37,11 @@ class SubentitySeedsCreator
    using EntityTraitsType      = typename MeshTraitsType::template EntityTraits< EntityTopology::dimension >;
    using SubentityTraits       = typename MeshTraitsType::template SubentityTraits< EntityTopology, SubentityDimensionTag::value >;
    using SubentityTopology     = typename SubentityTraits::SubentityTopology;
-
+   
 public:
-   using SubentitySeedArray = Containers::StaticArray< SubentityTraits::count, EntitySeed< MeshConfig, SubentityTopology > >;
+   using SubentitySeed = EntitySeed< MeshConfig, SubentityTopology >;
+   using SubentitySeedArray = Containers::StaticArray< SubentityTraits::count, SubentitySeed >;
+   using FunctorType = std::function< void( SubentitySeed& ) >;
 
    static SubentitySeedArray create( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
@@ -64,6 +66,27 @@ public:
       return subentitySeeds;
    }
 
+   static void iterate( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex, FunctorType&& functor )
+   {
+      const auto& subvertices = mesh.template getSubentitiesMatrix< EntityTopology::dimension, 0 >().getRow( entityIndex );
+
+      Algorithms::staticFor< LocalIndexType, 0, SubentitySeedArray::getSize() >(
+         [&] ( auto subentityIndex ) {
+            constexpr LocalIndexType subentityVerticesCount = Topologies::SubentityVertexCount< EntityTopology, SubentityTopology, subentityIndex >::count;
+            SubentitySeed subentitySeed;
+            subentitySeed.setCornersCount( subentityVerticesCount );
+            Algorithms::staticFor< LocalIndexType, 0, subentityVerticesCount >(
+               [&] ( auto subentityVertexIndex ) {
+                  // subentityIndex cannot be captured as constexpr, so we need to create another instance of its type
+                  static constexpr LocalIndexType VERTEX_INDEX = SubentityTraits::template Vertex< decltype(subentityIndex){}, subentityVertexIndex >::index;
+                  subentitySeed.setCornerId( subentityVertexIndex, subvertices.getColumnIndex( VERTEX_INDEX ) );
+               }
+            );
+            std::forward< FunctorType >( functor )( subentitySeed );
+         }
+      );
+   }
+
    constexpr static LocalIndexType getSubentitiesCount( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
       return SubentityTraits::count;
@@ -84,7 +107,9 @@ class SubentitySeedsCreator< MeshConfig, EntityTopology, DimensionTag< 0 > >
    using SubentityTopology     = typename SubentityTraits::SubentityTopology;
 
 public:
-   using SubentitySeedArray = Containers::StaticArray< SubentityTraits::count, EntitySeed< MeshConfig, SubentityTopology > >;
+   using SubentitySeed = EntitySeed< MeshConfig, SubentityTopology >;
+   using SubentitySeedArray = Containers::StaticArray< SubentityTraits::count, SubentitySeed >;
+   using FunctorType = std::function< void( SubentitySeed& ) >;
 
    static SubentitySeedArray create( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
@@ -96,6 +121,17 @@ public:
       return seeds;
    }
 
+   static void iterate( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex, FunctorType&& functor )
+   {
+      const auto& subvertices = mesh.template getSubentitiesMatrix< EntityTopology::dimension, 0 >().getRow( entityIndex );
+
+      for( LocalIndexType i = 0; i < SubentitySeedArray::getSize(); i++ ) {
+         SubentitySeed seed;
+         seed.setCornerId( 0, subvertices.getColumnIndex( i ) );
+         std::forward< FunctorType >( functor )( seed );
+      }
+   }
+
    constexpr static LocalIndexType getSubentitiesCount( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
       return SubentityTraits::count;
@@ -119,6 +155,7 @@ class SubentitySeedsCreator< MeshConfig, Topologies::Polygon, DimensionTag< 1 >
 public:
    using SubentitySeed = EntitySeed< MeshConfig, SubentityTopology >;
    using SubentitySeedArray = Containers::Array< SubentitySeed, DeviceType, LocalIndexType >;
+   using FunctorType = std::function< void( SubentitySeed& ) >;
    
    static SubentitySeedArray create( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
@@ -138,6 +175,20 @@ public:
       return seeds;
    }
 
+   static void iterate( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex, FunctorType&& functor )
+   {
+      const auto& subvertices = mesh.template getSubentitiesMatrix< EntityTopology::dimension, 0 >().getRow( entityIndex );
+      const LocalIndexType subverticesCount = mesh.template getSubentitiesCount< EntityTopology::dimension, 0 >( entityIndex );
+
+      for( LocalIndexType i = 0; i < subverticesCount; i++ )
+      {
+         SubentitySeed seed;
+         seed.setCornerId( 0, subvertices.getColumnIndex( i ) );
+         seed.setCornerId( 1, subvertices.getColumnIndex( (i + 1) % subverticesCount ) );
+         std::forward< FunctorType >( functor )( seed );
+      }
+   }
+
    static LocalIndexType getSubentitiesCount( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
       return mesh.template getSubentitiesCount< EntityTopology::dimension, 0 >( entityIndex );
@@ -150,6 +201,7 @@ class SubentitySeedsCreator< MeshConfig, Topologies::Polygon, DimensionTag< 0 >
    using MeshType              = Mesh< MeshConfig >;
    using MeshTraitsType        = MeshTraits< MeshConfig >;
    using InitializerType       = Initializer< MeshConfig >;
+   using DeviceType            = typename MeshTraitsType::DeviceType;
    using GlobalIndexType       = typename MeshTraitsType::GlobalIndexType;
    using LocalIndexType        = typename MeshTraitsType::LocalIndexType;
    using EntityTopology        = Topologies::Polygon;
@@ -158,7 +210,9 @@ class SubentitySeedsCreator< MeshConfig, Topologies::Polygon, DimensionTag< 0 >
    using SubentityTopology     = typename SubentityTraits::SubentityTopology;
 
 public:
-   using SubentitySeedArray = Containers::Array< EntitySeed< MeshConfig, SubentityTopology >, Devices::Host, LocalIndexType >;
+   using SubentitySeed = EntitySeed< MeshConfig, SubentityTopology >;
+   using SubentitySeedArray = Containers::Array< SubentitySeed, Devices::Host, LocalIndexType >;
+   using FunctorType = std::function< void( SubentitySeed& ) >;
 
    static SubentitySeedArray create( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
@@ -174,6 +228,18 @@ public:
       return seeds;
    }
 
+   static void iterate( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex, FunctorType&& functor )
+   {
+      const auto& subvertices = mesh.template getSubentitiesMatrix< EntityTopology::dimension, 0 >().getRow( entityIndex );
+      const LocalIndexType subverticesCount = mesh.template getSubentitiesCount< EntityTopology::dimension, 0 >( entityIndex );
+
+      for( LocalIndexType i = 0; i < subverticesCount; i++ ) {
+         SubentitySeed seed;
+         seed.setCornerId( 0, subvertices.getColumnIndex( i ) );
+         std::forward< FunctorType >( functor )( seed );
+      }
+   }
+
    static LocalIndexType getSubentitiesCount( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
       return mesh.template getSubentitiesCount< EntityTopology::dimension, 0 >( entityIndex );
@@ -197,6 +263,7 @@ class SubentitySeedsCreator< MeshConfig, Topologies::Polyhedron, DimensionTag< 2
 public:
    using SubentitySeed = EntitySeed< MeshConfig, SubentityTopology >;
    using SubentitySeedArray = Containers::Array< SubentitySeed, DeviceType, LocalIndexType >;
+   using FunctorType = std::function< void( SubentitySeed& ) >;
    
    static SubentitySeedArray create( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
@@ -222,6 +289,25 @@ public:
       return seeds;
    }
 
+   static void iterate( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex, FunctorType&& functor )
+   {
+      const auto& cellSeeds = initializer.getCellSeeds();
+      const auto& faces = cellSeeds[ entityIndex ].getCornerIds();
+
+      for( LocalIndexType i = 0; i < faces.getSize(); i++ )
+      {
+         GlobalIndexType faceIdx = faces[ i ];
+         const auto& subvertices = mesh.template getSubentitiesMatrix< 2, 0 >().getRow( faceIdx );
+         const LocalIndexType subverticesCount = mesh.template getSubentitiesCount< 2, 0 >( faceIdx );
+         SubentitySeed seed;
+         seed.setCornersCount( subverticesCount );
+         for( LocalIndexType j = 0; j < subverticesCount; j++ ) {
+            seed.setCornerId( j, subvertices.getColumnIndex( j ) );
+         }
+         std::forward< FunctorType >( functor )( seed );
+      }
+   }
+
    static LocalIndexType getSubentitiesCount( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
       auto& cellSeeds = initializer.getCellSeeds();
@@ -248,6 +334,7 @@ class SubentitySeedsCreator< MeshConfig, Topologies::Polyhedron, DimensionTag< 1
 public:
    using SubentitySeed = EntitySeed< MeshConfig, SubentityTopology >;
    using SubentitySeedArray = Containers::Array< SubentitySeed, DeviceType, LocalIndexType >;
+   using FunctorType = std::function< void( SubentitySeed& ) >;
    
    static SubentitySeedArray create( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
@@ -276,6 +363,22 @@ public:
       return seeds;
    }
 
+   static void iterate( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex, FunctorType&& functor )
+   {
+      SeedSet seedSet;
+      const auto& faces = mesh.template getSubentitiesMatrix< EntityTopology::dimension, 2 >().getRow( entityIndex );
+      const LocalIndexType facesCount = mesh.template getSubentitiesCount< EntityTopology::dimension, 2 >( entityIndex );
+
+      for( LocalIndexType i = 0; i < facesCount; i++ ) {
+         GlobalIndexType faceIdx = faces.getColumnIndex( i );
+         FaceSubentitySeedsCreator::iterate( initializer, mesh, faceIdx, [&] ( SubentitySeed& seed ) {
+            const bool inserted = seedSet.insert( seed ).second;
+            if( inserted )
+               std::forward< FunctorType >( functor )( seed );
+         });
+      }
+   }
+
    static LocalIndexType getSubentitiesCount( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
       SeedSet seedSet;
@@ -312,6 +415,7 @@ class SubentitySeedsCreator< MeshConfig, Topologies::Polyhedron, DimensionTag< 0
 public:
    using SubentitySeed = EntitySeed< MeshConfig, SubentityTopology >;
    using SubentitySeedArray = Containers::Array< SubentitySeed, DeviceType, LocalIndexType >;
+   using FunctorType = std::function< void( SubentitySeed& ) >;
    
    static SubentitySeedArray create( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
@@ -340,6 +444,22 @@ public:
       return seeds;
    }
 
+   static void iterate( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex, FunctorType&& functor )
+   {
+      SeedSet seedSet;
+      const auto& faces = mesh.template getSubentitiesMatrix< EntityTopology::dimension, 2 >().getRow( entityIndex );
+      const LocalIndexType facesCount = mesh.template getSubentitiesCount< EntityTopology::dimension, 2 >( entityIndex );
+
+      for( LocalIndexType i = 0; i < facesCount; i++ ) {
+         GlobalIndexType faceIdx = faces.getColumnIndex( i );
+         FaceSubentitySeedsCreator::iterate( initializer, mesh, faceIdx, [&] ( SubentitySeed& seed ) {
+            const bool inserted = seedSet.insert( seed ).second;
+            if( inserted )
+               std::forward< FunctorType >( functor )( seed );
+         });
+      }
+   }
+
    static LocalIndexType getSubentitiesCount( InitializerType& initializer, MeshType& mesh, const GlobalIndexType entityIndex )
    {
       SeedSet seedSet;