From c82d63d59c5fc8c8e9bc767d13f8727f50b4a8a0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jakub=20Klinkovsk=C3=BD?= <klinkjak@fjfi.cvut.cz>
Date: Fri, 28 Oct 2016 23:23:18 +0200
Subject: [PATCH] Storing subentities globally in multimaps, just like
 superentities

---
 src/TNL/Meshes/Mesh.h                         |   6 +-
 src/TNL/Meshes/MeshDetails/MeshEntity_impl.h  |  10 +-
 src/TNL/Meshes/MeshDetails/Mesh_impl.h        |   8 +-
 .../initializer/MeshEntityInitializer.h       |  30 +-
 .../MeshDetails/initializer/MeshInitializer.h | 105 ++--
 .../MeshSuperentityStorageInitializer.h       |   5 +-
 ...Rebinder.h => MeshEntityStorageRebinder.h} |  37 +-
 .../MeshDetails/layers/MeshStorageLayer.h     | 115 ++--
 .../MeshDetails/layers/MeshSubentityAccess.h  | 515 ++++++++++++++++++
 .../layers/MeshSubentityStorageLayer.h        | 346 +++---------
 .../layers/MeshSuperentityAccess.h            |   5 +-
 .../layers/MeshSuperentityStorageLayer.h      |  91 ++--
 .../MeshDetails/traits/MeshEntityTraits.h     |   6 +-
 .../MeshDetails/traits/MeshSubentityTraits.h  |   7 +
 src/TNL/Meshes/MeshEntity.h                   |  21 +-
 src/UnitTests/Meshes/MeshEntityTest.h         | 119 +++-
 16 files changed, 931 insertions(+), 495 deletions(-)
 rename src/TNL/Meshes/MeshDetails/layers/{MeshSuperentityStorageRebinder.h => MeshEntityStorageRebinder.h} (53%)
 create mode 100644 src/TNL/Meshes/MeshDetails/layers/MeshSubentityAccess.h

diff --git a/src/TNL/Meshes/Mesh.h b/src/TNL/Meshes/Mesh.h
index 518dc6b917..e1f7b662fe 100644
--- a/src/TNL/Meshes/Mesh.h
+++ b/src/TNL/Meshes/Mesh.h
@@ -23,7 +23,7 @@
 #include <TNL/Meshes/MeshDetails/layers/MeshStorageLayer.h>
 #include <TNL/Meshes/MeshDetails/config/MeshConfigValidator.h>
 #include <TNL/Meshes/MeshDetails/initializer/MeshInitializer.h>
-#include <TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageRebinder.h>
+#include <TNL/Meshes/MeshDetails/layers/MeshEntityStorageRebinder.h>
 
 namespace TNL {
 namespace Meshes {
@@ -94,7 +94,7 @@ class Mesh
 
    protected:
       // Methods for the mesh initializer
-      using StorageBaseType::getEntitiesArray;
+      using StorageBaseType::getSubentityStorageNetwork;
       using StorageBaseType::getSuperentityStorageNetwork;
 
       MeshConfigValidator< MeshConfig > configValidator;
@@ -102,7 +102,7 @@ class Mesh
       friend MeshInitializer< MeshConfig >;
 
       template< typename Mesh, typename DimensionsTag, typename SuperdimensionsTag >
-      friend struct MeshSuperentityStorageRebinderWorker;
+      friend struct MeshEntityStorageRebinderWorker;
 };
 
 template< typename MeshConfig >
diff --git a/src/TNL/Meshes/MeshDetails/MeshEntity_impl.h b/src/TNL/Meshes/MeshDetails/MeshEntity_impl.h
index 90e510680e..f037a87101 100644
--- a/src/TNL/Meshes/MeshDetails/MeshEntity_impl.h
+++ b/src/TNL/Meshes/MeshDetails/MeshEntity_impl.h
@@ -45,7 +45,7 @@ bool
 MeshEntity< MeshConfig, EntityTopology >::
 save( File& file ) const
 {
-   if( ! MeshSubentityStorageLayers< MeshConfig, EntityTopology >::save( file ) /*||
+   if( ! MeshSubentityAccess< MeshConfig, EntityTopology >::save( file ) /*||
        ! MeshSuperentityStorageLayers< MeshConfig, EntityTopology >::save( file )*/ )
       return false;
    return true;
@@ -57,7 +57,7 @@ bool
 MeshEntity< MeshConfig, EntityTopology >::
 load( File& file )
 {
-   if( ! MeshSubentityStorageLayers< MeshConfig, EntityTopology >::load( file ) /*||
+   if( ! MeshSubentityAccess< MeshConfig, EntityTopology >::load( file ) /*||
        ! MeshSuperentityStorageLayers< MeshConfig, EntityTopology >::load( file ) */ )
       return false;
    return true;
@@ -69,8 +69,8 @@ void
 MeshEntity< MeshConfig, EntityTopology >::
 print( std::ostream& str ) const
 {
-   str << "\t Mesh entity dimension: " << EntityTopology::dimensions << std::endl;
-   MeshSubentityStorageLayers< MeshConfig, EntityTopology >::print( str );
+   str << "\t Mesh entity dimensions: " << EntityTopology::dimensions << std::endl;
+   MeshSubentityAccess< MeshConfig, EntityTopology >::print( str );
    MeshSuperentityAccess< MeshConfig, EntityTopology >::print( str );
 }
 
@@ -80,7 +80,7 @@ bool
 MeshEntity< MeshConfig, EntityTopology >::
 operator==( const MeshEntity& entity ) const
 {
-   return ( MeshSubentityStorageLayers< MeshConfig, EntityTopology >::operator==( entity ) &&
+   return ( MeshSubentityAccess< MeshConfig, EntityTopology >::operator==( entity ) &&
             MeshSuperentityAccess< MeshConfig, EntityTopology >::operator==( entity ) &&
             MeshEntityId< typename MeshConfig::IdType,
                           typename MeshConfig::GlobalIndexType >::operator==( entity ) );
diff --git a/src/TNL/Meshes/MeshDetails/Mesh_impl.h b/src/TNL/Meshes/MeshDetails/Mesh_impl.h
index 4ba0a0505b..1ba98b18b5 100644
--- a/src/TNL/Meshes/MeshDetails/Mesh_impl.h
+++ b/src/TNL/Meshes/MeshDetails/Mesh_impl.h
@@ -131,7 +131,7 @@ load( File& file )
       return false;
    }
    // update pointers from entities into the superentity storage network
-   MeshSuperentityStorageRebinder< Mesh< MeshConfig > >::exec( *this );
+   MeshEntityStorageRebinder< Mesh< MeshConfig > >::exec( *this );
    return true;
 }
 
@@ -154,14 +154,14 @@ operator==( const Mesh& mesh ) const
 template< typename MeshConfig >
 bool
 Mesh< MeshConfig >::
-init( const typename Mesh< MeshConfig >::MeshTraitsType::PointArrayType& points,
-      const typename Mesh< MeshConfig >::MeshTraitsType::CellSeedArrayType& cellSeeds )
+init( const typename MeshTraitsType::PointArrayType& points,
+      const typename MeshTraitsType::CellSeedArrayType& cellSeeds )
 {
    MeshInitializer< MeshConfig> meshInitializer;
    if( ! meshInitializer.createMesh( points, cellSeeds, *this ) )
       return false;
    // update pointers from entities into the superentity storage network
-   MeshSuperentityStorageRebinder< Mesh< MeshConfig > >::exec( *this );
+   MeshEntityStorageRebinder< Mesh< MeshConfig > >::exec( *this );
    return true;
 }
 
diff --git a/src/TNL/Meshes/MeshDetails/initializer/MeshEntityInitializer.h b/src/TNL/Meshes/MeshDetails/initializer/MeshEntityInitializer.h
index 1ee6bc64de..9f59fea6a2 100644
--- a/src/TNL/Meshes/MeshDetails/initializer/MeshEntityInitializer.h
+++ b/src/TNL/Meshes/MeshDetails/initializer/MeshEntityInitializer.h
@@ -56,11 +56,10 @@ class MeshEntityInitializer
                                                                        EntityTopology,
                                                                        typename MeshTraits< MeshConfig >::DimensionsTag >;
 
-   static constexpr int Dimensions = DimensionsTag::value;
    using MeshTraitsType   = MeshTraits< MeshConfig >;
    using GlobalIndexType  = typename MeshTraitsType::GlobalIndexType;
    using LocalIndexType   = typename MeshTraitsType::LocalIndexType;
-   using EntityTraitsType = typename MeshTraitsType::template EntityTraits< Dimensions >;
+   using EntityTraitsType = typename MeshTraitsType::template EntityTraits< DimensionsTag::value >;
 
    using EntityType       = typename EntityTraitsType::EntityType;
    using SubvertexTraits  = typename MeshTraitsType::template SubentityTraits< EntityTopology, 0 >;
@@ -76,7 +75,7 @@ public:
    static void initEntity( EntityType& entity, GlobalIndexType entityIndex, const SeedType& entitySeed, InitializerType& initializer)
    {
       for( LocalIndexType i = 0; i < entitySeed.getCornerIds().getSize(); i++ )
-         initializer.template setSubentityIndex< 0 >( entity, i, entitySeed.getCornerIds()[ i ] );
+         initializer.template setSubentityIndex< 0 >( entity, entityIndex, i, entitySeed.getCornerIds()[ i ] );
       BaseType::initSubentities( entity, entityIndex, entitySeed, initializer );
    }
 
@@ -134,9 +133,8 @@ class MeshEntityInitializerLayer< MeshConfig,
                                                 EntityTopology,
                                                 typename DimensionsTag::Decrement >;
 
-   static constexpr int Dimensions = DimensionsTag::value;
    using MeshTraitsType            = MeshTraits< MeshConfig >;
-   using SubentityTraitsType       = typename MeshTraitsType::template SubentityTraits< EntityTopology, Dimensions >;
+   using SubentityTraitsType       = typename MeshTraitsType::template SubentityTraits< EntityTopology, DimensionsTag::value >;
    using GlobalIndexType           = typename SubentityTraitsType::GlobalIndexType;
 
    using InitializerType           = MeshInitializer< MeshConfig >;
@@ -157,7 +155,7 @@ protected:
       for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )
       {
          const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( subentitySeeds[ i ] );
-         InitializerType::template setSubentityIndex< DimensionsTag::value >( entity, i, subentityIndex );
+         meshInitializer.template setSubentityIndex< DimensionsTag::value >( entity, entityIndex, i, subentityIndex );
          //std::cout << "    Adding subentity " << subentityIndex << std::endl;
          meshInitializer.
             template getSuperentityInitializer< DimensionsTag >().
@@ -214,7 +212,7 @@ protected:
       for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )
       {
          const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( subentitySeeds[ i ] );
-         InitializerType::template setSubentityIndex< DimensionsTag::value >( entity, i, subentityIndex );
+         meshInitializer.template setSubentityIndex< DimensionsTag::value >( entity, entityIndex, i, subentityIndex );
          //std::cout << "    Adding subentity " << subentityIndex << std::endl;
          subentityOrientationsArray[ i ] = meshInitializer.template getReferenceOrientation< DimensionsTag >( subentityIndex ).createOrientation( subentitySeeds[ i ] );
          //std::cout << "    Subentity orientation = " << subentityOrientationsArray[ i ].getSubvertexPermutation() << std::endl;
@@ -273,7 +271,7 @@ protected:
       for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )
       {
          const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( subentitySeeds[ i ] );
-         InitializerType::template setSubentityIndex< DimensionsTag::value >( entity, i, subentityIndex );
+         meshInitializer.template setSubentityIndex< DimensionsTag::value >( entity, entityIndex, i, subentityIndex );
          //std::cout << "    Adding subentity " << subentityIndex << std::endl;
          subentityOrientationsArray[ i ] = meshInitializer.template getReferenceOrientation< DimensionsTag >( subentitySeeds[ i ] ).createOrientation( subentitySeeds[ i ] );
       }
@@ -318,7 +316,7 @@ protected:
       for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++)
       {
          const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( subentitySeeds[ i ] );
-         InitializerType::template setSubentityIndex< DimensionsTag::value >( entity, i, subentityIndex );
+         meshInitializer.template setSubentityIndex< DimensionsTag::value >( entity, entityIndex, i, subentityIndex );
       }
 
       BaseType::initSubentities(entity, entityIndex, entitySeed, meshInitializer);
@@ -411,7 +409,7 @@ protected:
       //std::cout << "   Initiating subentities with " << DimensionsTag::value << " dimensions ... " << std::endl;
       for( LocalIndexType i = 0; i < entity.template getNumberOfSubentities< DimensionsTag::value >(); i++ )
       {
-         const GlobalIndexType subentityIndex = entity.template getSubentityIndex< DimensionsTag::value >( i );
+         const GlobalIndexType subentityIndex = meshInitializer.template getSubentityIndex< DimensionsTag::value >( entity, entityIndex, i );
          meshInitializer.template getSuperentityInitializer< DimensionsTag >().addSuperentity( EntityDimensionsTag(), subentityIndex, entityIndex);
       }
    }
@@ -426,12 +424,12 @@ class MeshEntityInitializerLayer< MeshConfig,
                                   false,
                                   false >
 {
-   using InitializerType     = MeshInitializer< MeshConfig >;
-   using DimensionsTag       = MeshDimensionsTag< 0 >;
-   using SubentitiesTraits   = MeshSubentityTraits< MeshConfig, EntityTopology, DimensionsTag::value >;
-   using GlobalIndexType     = typename SubentitiesTraits::GlobalIndexType;
-   using EntityType          = MeshEntity< MeshConfig, EntityTopology >;
-   using SeedType            = MeshEntitySeed< MeshConfig, EntityTopology >;
+   using InitializerType   = MeshInitializer< MeshConfig >;
+   using DimensionsTag     = MeshDimensionsTag< 0 >;
+   using SubentitiesTraits = MeshSubentityTraits< MeshConfig, EntityTopology, DimensionsTag::value >;
+   using GlobalIndexType   = typename SubentitiesTraits::GlobalIndexType;
+   using EntityType        = MeshEntity< MeshConfig, EntityTopology >;
+   using SeedType          = MeshEntitySeed< MeshConfig, EntityTopology >;
 
 protected:
    static void initSubentities( EntityType& entity, GlobalIndexType entityIndex, const SeedType& entitySeed,
diff --git a/src/TNL/Meshes/MeshDetails/initializer/MeshInitializer.h b/src/TNL/Meshes/MeshDetails/initializer/MeshInitializer.h
index 99757900a9..9b92e94ac7 100644
--- a/src/TNL/Meshes/MeshDetails/initializer/MeshInitializer.h
+++ b/src/TNL/Meshes/MeshDetails/initializer/MeshInitializer.h
@@ -54,18 +54,12 @@ class MeshInitializer
    public:
       using MeshType          = Mesh< MeshConfig >;
       using MeshTraitsType    = MeshTraits< MeshConfig >;
-      static constexpr int Dimensions = MeshTraitsType::meshDimensions;
-      using DimensionsTag     = MeshDimensionsTag< Dimensions >;
+      using DimensionsTag     = MeshDimensionsTag< MeshTraitsType::meshDimensions >;
       using BaseType          = MeshInitializerLayer< MeshConfig, DimensionsTag >;
       using PointArrayType    = typename MeshTraitsType::PointArrayType;
       using CellSeedArrayType = typename MeshTraitsType::CellSeedArrayType;
       using GlobalIndexType   = typename MeshTraitsType::GlobalIndexType;
 
-      template< typename DimensionsTag, typename SuperdimensionsTag >
-      using SuperentityStorageNetwork = typename MeshTraitsType::template SuperentityTraits<
-            typename MeshTraitsType::template EntityTraits< DimensionsTag::value >::EntityTopology,
-            SuperdimensionsTag::value >::StorageNetworkType;
-
 
       MeshInitializer()
       : verbose( false ), mesh( 0 )
@@ -90,16 +84,33 @@ class MeshInitializer
          BaseType::createEntityReferenceOrientations();
 
          if( verbose ) std::cout << "====== Initiating entities ==============" << std::endl;
-         BaseType::initEntities( *this, points, cellSeeds );
+         BaseType::initEntities( *this, points, cellSeeds, mesh );
 
          return true;
       }
 
-      template<int Subdimensions, typename EntityType, typename LocalIndex, typename GlobalIndex >
-      static void
-      setSubentityIndex( EntityType& entity, const LocalIndex& localIndex, const GlobalIndex& globalIndex )
+      template< int Dimensions >
+      bool setNumberOfEntities( const GlobalIndexType& entitiesCount )
+      {
+         return mesh->template setNumberOfEntities< Dimensions >( entitiesCount );
+      }
+
+      template< int Subdimensions, typename EntityType, typename LocalIndex, typename GlobalIndex >
+      void
+      setSubentityIndex( const EntityType& entity, const GlobalIndex& entityIndex, const LocalIndex& localIndex, const GlobalIndex& globalIndex )
+      {
+         // The mesh entities are not yet bound to the storage network at this point,
+         // so we operate directly on the storage.
+         mesh->template getSubentityStorageNetwork< EntityType::EntityTopology::dimensions, Subdimensions >().getValues( entityIndex )[ localIndex ] = globalIndex;
+      }
+
+      template< int Subdimensions, typename EntityType, typename LocalIndex, typename GlobalIndex >
+      GlobalIndex
+      getSubentityIndex( const EntityType& entity, const GlobalIndex& entityIndex, const LocalIndex& localIndex )
       {
-         entity.template setSubentityIndex< Subdimensions >( localIndex, globalIndex );
+         // The mesh entities are not yet bound to the storage network at this point,
+         // so we operate directly on the storage.
+         return mesh->template getSubentityStorageNetwork< EntityType::EntityTopology::dimensions, Subdimensions >().getValues( entityIndex )[ localIndex ];
       }
 
       template<typename SubDimensionsTag, typename MeshEntity >
@@ -109,18 +120,11 @@ class MeshInitializer
          return entity.template subentityOrientationsArray< SubDimensionTag::value >();
       }
 
-      template< int Dimensions >
-      typename MeshTraitsType::template EntityTraits< Dimensions >::StorageArrayType&
-      meshEntitiesArray()
-      {
-         return mesh->template getEntitiesArray< Dimensions >();
-      }
-
-      template< typename EntityTopology, typename SuperdimensionsTag >
-      typename MeshTraitsType::template SuperentityTraits< EntityTopology, SuperdimensionsTag::value >::StorageNetworkType&
+      template< typename EntityTopology, int Superdimensions >
+      typename MeshTraitsType::template SuperentityTraits< EntityTopology, Superdimensions >::StorageNetworkType&
       meshSuperentityStorageNetwork()
       {
-         return mesh->template getSuperentityStorageNetwork< EntityTopology, SuperdimensionsTag >();
+         return mesh->template getSuperentityStorageNetwork< EntityTopology::dimensions, Superdimensions >();
       }
 
       static void
@@ -139,9 +143,9 @@ class MeshInitializer
 
       template< typename DimensionsTag >
       const MeshEntityReferenceOrientation< MeshConfig, typename MeshTraitsType::template EntityTraits< DimensionsTag::value >::EntityTopology >&
-      getReferenceOrientation( GlobalIndexType index) const
+      getReferenceOrientation( GlobalIndexType index ) const
       {
-         return BaseType::getReferenceOrientation( DimensionTag(), index);
+         return BaseType::getReferenceOrientation( DimensionsTag(), index );
       }
 
    protected:
@@ -191,21 +195,20 @@ class MeshInitializerLayer< MeshConfig,
          BaseType::createEntitySeedsFromCellSeeds( cellSeeds );
       }
 
-      void initEntities( InitializerType& initializer, const PointArrayType& points, const CellSeedArrayType& cellSeeds)
+      void initEntities( InitializerType& initializer, const PointArrayType& points, const CellSeedArrayType& cellSeeds, MeshType& mesh )
       {
-         StorageArrayType& entityArray = initializer.template meshEntitiesArray< Dimensions >();
-         //cout << " Initiating entities with " << DimensionsTag::value << " dimensions ... " << std::endl;
-         entityArray.setSize( cellSeeds.getSize() );
-         for( GlobalIndexType i = 0; i < entityArray.getSize(); i++ )
+         //std::cout << " Initiating entities with " << DimensionsTag::value << " dimensions ... " << std::endl;
+         initializer.template setNumberOfEntities< Dimensions >( cellSeeds.getSize() );
+         for( GlobalIndexType i = 0; i < cellSeeds.getSize(); i++ )
          {
-            //cout << "  Initiating entity " << i << std::endl;
-            EntityInitializerType::initEntity( entityArray[i], i, cellSeeds[i], initializer );
+            //std::cout << "  Initiating entity " << i << std::endl;
+            EntityInitializerType::initEntity( mesh.template getEntity< Dimensions >( i ), i, cellSeeds[i], initializer );
          }
          /***
           * There are no superentities in this layer storing mesh cells.
           */
 
-         BaseType::initEntities( initializer, points );
+         BaseType::initEntities( initializer, points, mesh );
       }
 
       using BaseType::findEntitySeedIndex;
@@ -270,11 +273,10 @@ class MeshInitializerLayer< MeshConfig,
                                   typename DimensionsTag::Decrement >
 {
    using MeshTraitsType               = MeshTraits< MeshConfig >;
-   static constexpr int Dimensions = DimensionsTag::value;
    using BaseType                     = MeshInitializerLayer< MeshConfig, typename DimensionsTag::Decrement >;
 
    using MeshType                     = Mesh< MeshConfig >;
-   using EntityTraitsType             = typename MeshTraitsType::template EntityTraits< Dimensions >;
+   using EntityTraitsType             = typename MeshTraitsType::template EntityTraits< DimensionsTag::value >;
    using EntityTopology               = typename EntityTraitsType::EntityTopology;
    using GlobalIndexType              = typename MeshTraitsType::GlobalIndexType;
    using CellTopology                 = typename MeshTraitsType::CellTopology;
@@ -325,24 +327,23 @@ class MeshInitializerLayer< MeshConfig,
          return this->superentityInitializer;
       }
 
-      void initEntities( InitializerType& initializer, const PointArrayType& points )
+      void initEntities( InitializerType& initializer, const PointArrayType& points, MeshType& mesh )
       {
-         StorageArrayType& entityArray = initializer.template meshEntitiesArray< Dimensions >();
-         //cout << " Initiating entities with " << DimensionsTag::value << " dimensions ... " << std::endl;
-         entityArray.setSize( this->seedsIndexedSet.getSize() );
+         //std::cout << " Initiating entities with " << DimensionsTag::value << " dimensions ... " << std::endl;
+         initializer.template setNumberOfEntities< DimensionsTag::value >( this->seedsIndexedSet.getSize() );
          EntitySeedArrayType seedsArray;
          seedsArray.setSize( this->seedsIndexedSet.getSize() );
          this->seedsIndexedSet.toArray( seedsArray );
          for( GlobalIndexType i = 0; i < this->seedsIndexedSet.getSize(); i++ )
          {
-            //cout << "  Initiating entity " << i << std::endl;
-            EntityInitializerType::initEntity( entityArray[ i ], i, seedsArray[ i ], initializer );
+            //std::cout << "  Initiating entity " << i << std::endl;
+            EntityInitializerType::initEntity( mesh.template getEntity< DimensionsTag::value >( i ), i, seedsArray[ i ], initializer );
          }
          this->seedsIndexedSet.reset();
 
          this->superentityInitializer.initSuperentities( initializer );
 
-         BaseType::initEntities(initializer, points);
+         BaseType::initEntities( initializer, points, mesh );
       }
 
       void createEntityReferenceOrientations()
@@ -427,28 +428,27 @@ class MeshInitializerLayer< MeshConfig,
          return this->superentityInitializer;
       }
 
-      void initEntities( InitializerType& initializer, const PointArrayType& points )
+      void initEntities( InitializerType& initializer, const PointArrayType& points, MeshType& mesh )
       {
-         EntityArrayType& entityArray = initializer.template meshEntitiesArray< DimensionsTag::value >();
          //cout << " Initiating entities with " << DimensionsTag::value << " dimensions ... " << std::endl;
-         entityArray.setSize( this->seedsIndexedSet.getSize() );
+         initializer.template setNumberOfEntities< DimensionsTag::value >( this->seedsIndexedSet.getSize() );
          SeedArrayType seedsArray;
          seedsArray.setSize( this->seedsIndexedSet.getSize() );
          this->seedsIndexedSet.toArray( seedsArray );
          for( GlobalIndexType i = 0; i < this->seedsIndexedSet.getSize(); i++ )
          {
-            //cout << "  Initiating entity " << i << std::endl;
-            EntityInitializerType::initEntity( entityArray[ i ], i, seedsArray[ i ], initializer );
+            //std::cout << "  Initiating entity " << i << std::endl;
+            EntityInitializerType::initEntity( mesh.template getEntity< DimensionsTag::value >( i ), i, seedsArray[ i ], initializer );
          }
          this->seedsIndexedSet.reset();
 
          this->superentityInitializer.initSuperentities( initializer );
 
-         BaseType::initEntities(initializer, points);
+         BaseType::initEntities( initializer, points, mesh );
       }
 
       using BaseType::getReferenceOrientation;
-      const ReferenceOrientationType& getReferenceOrientation( DimensionTag, GlobalIndexType index) const
+      const ReferenceOrientationType& getReferenceOrientation( DimensionsTag, GlobalIndexType index ) const
       {
          return this->referenceOrientations[ index ];
       }
@@ -534,12 +534,11 @@ class MeshInitializerLayer< MeshConfig,
 
       void createEntitySeedsFromCellSeeds( const CellSeedArrayType& cellSeeds ) {}
 
-      void initEntities( InitializerType& initializer, const PointArrayType& points )
+      void initEntities( InitializerType& initializer, const PointArrayType& points, MeshType& mesh )
       {
-         EntityArrayType& vertexArray = initializer.template meshEntitiesArray< DimensionsTag::value >();
-         vertexArray.setSize( points.getSize() );
-         for( GlobalIndexType i = 0; i < vertexArray.getSize(); i++ )
-            EntityInitializerType::setVertexPoint( vertexArray[i], points[i], initializer );
+         initializer.template setNumberOfEntities< 0 >( points.getSize() );
+         for( GlobalIndexType i = 0; i < points.getSize(); i++ )
+            EntityInitializerType::setVertexPoint( mesh.template getEntity< 0 >( i ), points[i], initializer );
 
          superentityInitializer.initSuperentities( initializer );
       }
diff --git a/src/TNL/Meshes/MeshDetails/initializer/MeshSuperentityStorageInitializer.h b/src/TNL/Meshes/MeshDetails/initializer/MeshSuperentityStorageInitializer.h
index c6c20d91fe..c54f3ec746 100644
--- a/src/TNL/Meshes/MeshDetails/initializer/MeshSuperentityStorageInitializer.h
+++ b/src/TNL/Meshes/MeshDetails/initializer/MeshSuperentityStorageInitializer.h
@@ -55,7 +55,6 @@ class MeshSuperentityStorageInitializerLayer< MeshConfig,
                                                             EntityTopology,
                                                             typename DimensionsTag::Decrement >;
 
-   static const int Dimensions = DimensionsTag::value;
    using EntityDimensions          = MeshDimensionsTag< EntityTopology::dimensions >;
    using EntityType                = MeshEntity< MeshConfig, EntityTopology >;
 
@@ -63,7 +62,7 @@ class MeshSuperentityStorageInitializerLayer< MeshConfig,
    using GlobalIndexType           = typename MeshTraitsType::GlobalIndexType;
    using LocalIndexType            = typename MeshTraitsType::LocalIndexType;
    using MeshInitializerType       = MeshInitializer< MeshConfig >;
-   using SuperentityTraitsType     = typename MeshTraitsType::template SuperentityTraits< EntityTopology, Dimensions >;
+   using SuperentityTraitsType     = typename MeshTraitsType::template SuperentityTraits< EntityTopology, DimensionsTag::value >;
    using SuperentityStorageNetwork = typename SuperentityTraitsType::StorageNetworkType;
 
    public:
@@ -96,7 +95,7 @@ class MeshSuperentityStorageInitializerLayer< MeshConfig,
             /****
              * Network initializer
              */
-            SuperentityStorageNetwork& superentityStorageNetwork = meshInitializer.template meshSuperentityStorageNetwork< EntityTopology, DimensionsTag >();
+            SuperentityStorageNetwork& superentityStorageNetwork = meshInitializer.template meshSuperentityStorageNetwork< EntityTopology, DimensionsTag::value >();
             superentityStorageNetwork.setKeysRange( maxEntityIndex + 1 );
             typename SuperentityStorageNetwork::ValuesAllocationVectorType storageNetworkAllocationVector;
             storageNetworkAllocationVector.setSize( maxEntityIndex + 1 );
diff --git a/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageRebinder.h b/src/TNL/Meshes/MeshDetails/layers/MeshEntityStorageRebinder.h
similarity index 53%
rename from src/TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageRebinder.h
rename to src/TNL/Meshes/MeshDetails/layers/MeshEntityStorageRebinder.h
index 1339e295e1..265f7bcbcd 100644
--- a/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageRebinder.h
+++ b/src/TNL/Meshes/MeshDetails/layers/MeshEntityStorageRebinder.h
@@ -1,5 +1,5 @@
 /***************************************************************************
-                          MeshSuperentityStorageRebinder.h  -  description
+                          MeshEntityStorageRebinder.h  -  description
                              -------------------
     begin                : Oct 22, 2016
     copyright            : (C) 2014 by Tomas Oberhuber et al.
@@ -31,56 +31,63 @@ namespace TNL {
 namespace Meshes {
 
 template< typename Mesh, typename DimensionsTag, typename SuperdimensionsTag >
-struct MeshSuperentityStorageRebinderWorker
+struct MeshEntityStorageRebinderWorker
 {
    static void exec( Mesh& mesh )
    {
+      for( typename Mesh::GlobalIndexType i = 0; i < mesh.template getNumberOfEntities< SuperdimensionsTag::value >(); i++ )
+      {
+         auto& superentity = mesh.template getEntity< SuperdimensionsTag::value >( i );
+         auto& subentitiesStorage = mesh.template getSubentityStorageNetwork< SuperdimensionsTag::value, DimensionsTag::value >();
+         superentity.template bindSubentitiesStorageNetwork< DimensionsTag::value >( subentitiesStorage.getValues( i ) );
+      }
+
       for( typename Mesh::GlobalIndexType i = 0; i < mesh.template getNumberOfEntities< DimensionsTag::value >(); i++ )
       {
-         auto& entity = mesh.template getEntity< DimensionsTag::value >( i );
-         auto& storage = mesh.template getSuperentityStorageNetwork< typename Mesh::template EntityTraits< DimensionsTag::value >::EntityTopology, SuperdimensionsTag >();
-         entity.template bindSuperentitiesStorageNetwork< SuperdimensionsTag::value >( storage.getValues( i ) );
+         auto& subentity = mesh.template getEntity< DimensionsTag::value >( i );
+         auto& superentitiesStorage = mesh.template getSuperentityStorageNetwork< DimensionsTag::value, SuperdimensionsTag::value >();
+         subentity.template bindSuperentitiesStorageNetwork< SuperdimensionsTag::value >( superentitiesStorage.getValues( i ) );
       }
    }
 };
 
 
 template< typename Mesh, typename DimensionsTag, typename SuperdimensionsTag >
-struct MeshSuperentityStorageRebinderInner
+struct MeshEntityStorageRebinderInner
 {
    static void exec( Mesh& mesh )
    {
-      MeshSuperentityStorageRebinderWorker< Mesh, DimensionsTag, SuperdimensionsTag >::exec( mesh );
-      MeshSuperentityStorageRebinderInner< Mesh, DimensionsTag, typename SuperdimensionsTag::Decrement >::exec( mesh );
+      MeshEntityStorageRebinderWorker< Mesh, DimensionsTag, SuperdimensionsTag >::exec( mesh );
+      MeshEntityStorageRebinderInner< Mesh, DimensionsTag, typename SuperdimensionsTag::Decrement >::exec( mesh );
    }
 };
 
 template< typename Mesh, typename SuperdimensionsTag >
-struct MeshSuperentityStorageRebinderInner< Mesh, typename SuperdimensionsTag::Decrement, SuperdimensionsTag >
+struct MeshEntityStorageRebinderInner< Mesh, typename SuperdimensionsTag::Decrement, SuperdimensionsTag >
 {
    static void exec( Mesh& mesh )
    {
-      MeshSuperentityStorageRebinderWorker< Mesh, typename SuperdimensionsTag::Decrement, SuperdimensionsTag >::exec( mesh );
+      MeshEntityStorageRebinderWorker< Mesh, typename SuperdimensionsTag::Decrement, SuperdimensionsTag >::exec( mesh );
    }
 };
 
 
 template< typename Mesh, typename DimensionsTag = typename MeshDimensionsTag< Mesh::MeshTraitsType::meshDimensions >::Decrement >
-struct MeshSuperentityStorageRebinder
+struct MeshEntityStorageRebinder
 {
    static void exec( Mesh& mesh )
    {
-      MeshSuperentityStorageRebinderInner< Mesh, DimensionsTag, MeshDimensionsTag< Mesh::MeshTraitsType::meshDimensions > >::exec( mesh );
-      MeshSuperentityStorageRebinder< Mesh, typename DimensionsTag::Decrement >::exec( mesh );
+      MeshEntityStorageRebinderInner< Mesh, DimensionsTag, MeshDimensionsTag< Mesh::MeshTraitsType::meshDimensions > >::exec( mesh );
+      MeshEntityStorageRebinder< Mesh, typename DimensionsTag::Decrement >::exec( mesh );
    }
 };
 
 template< typename Mesh >
-struct MeshSuperentityStorageRebinder< Mesh, MeshDimensionsTag< 0 > >
+struct MeshEntityStorageRebinder< Mesh, MeshDimensionsTag< 0 > >
 {
    static void exec( Mesh& mesh )
    {
-      MeshSuperentityStorageRebinderInner< Mesh, MeshDimensionsTag< 0 >, MeshDimensionsTag< Mesh::MeshTraitsType::meshDimensions > >::exec( mesh );
+      MeshEntityStorageRebinderInner< Mesh, MeshDimensionsTag< 0 >, MeshDimensionsTag< Mesh::MeshTraitsType::meshDimensions > >::exec( mesh );
    }
 };
 
diff --git a/src/TNL/Meshes/MeshDetails/layers/MeshStorageLayer.h b/src/TNL/Meshes/MeshDetails/layers/MeshStorageLayer.h
index 706d65e1e5..32c56f01e1 100644
--- a/src/TNL/Meshes/MeshDetails/layers/MeshStorageLayer.h
+++ b/src/TNL/Meshes/MeshDetails/layers/MeshStorageLayer.h
@@ -19,6 +19,7 @@
 #include <TNL/File.h>
 #include <TNL/Meshes/MeshDetails/traits/MeshTraits.h>
 #include <TNL/Meshes/MeshDetails/traits/MeshEntityTraits.h>
+#include <TNL/Meshes/MeshDetails/layers/MeshSubentityStorageLayer.h>
 #include <TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageLayer.h>
 
 namespace TNL {
@@ -41,16 +42,29 @@ class MeshStorageLayers
 
 protected:
    template< int Dimensions >
-   typename EntityTraits< Dimensions >::StorageArrayType& getEntitiesArray()
+   bool setNumberOfEntities( const typename EntityTraits< Dimensions >::GlobalIndexType& entitiesCount )
    {
-      return BaseType::getEntitiesArray( MeshDimensionsTag< Dimensions >() );
+      return BaseType::setNumberOfEntities( MeshDimensionsTag< Dimensions >(), entitiesCount );
    }
 
-   template< typename EntityTopology, typename SuperdimensionsTag >
-   typename MeshTraitsType::template SuperentityTraits< EntityTopology, SuperdimensionsTag::value >::StorageNetworkType&
+   template< int Dimensions, int Subdimensions >
+   typename MeshTraitsType::template SubentityTraits< typename EntityTraits< Dimensions >::EntityTopology, Subdimensions >::StorageNetworkType&
+   getSubentityStorageNetwork()
+   {
+      static_assert( Dimensions > Subdimensions, "Invalid combination of Dimensions and Subdimensions." );
+      using BaseType = MeshSubentityStorageLayers< MeshConfig,
+                                                   typename MeshTraits< MeshConfig >::template EntityTraits< Dimensions >::EntityTopology >;
+      return BaseType::template getSubentityStorageNetwork< Subdimensions >();
+   }
+
+   template< int Dimensions, int Superdimensions >
+   typename MeshTraitsType::template SuperentityTraits< typename EntityTraits< Dimensions >::EntityTopology, Superdimensions >::StorageNetworkType&
    getSuperentityStorageNetwork()
    {
-      return BaseType::template getSuperentityStorageNetwork< SuperdimensionsTag >( MeshDimensionsTag< EntityTopology::dimensions >() );
+      static_assert( Dimensions < Superdimensions, "Invalid combination of Dimensions and Superdimensions." );
+      using BaseType = MeshSuperentityStorageLayers< MeshConfig,
+                                                     typename MeshTraits< MeshConfig >::template EntityTraits< Dimensions >::EntityTopology >;
+      return BaseType::template getSuperentityStorageNetwork< Superdimensions >( MeshDimensionsTag< Dimensions >() );
    }
 };
 
@@ -61,6 +75,8 @@ class MeshStorageLayer< MeshConfig,
                         DimensionsTag,
                         true >
    : public MeshStorageLayer< MeshConfig, typename DimensionsTag::Decrement >,
+     public MeshSubentityStorageLayers< MeshConfig,
+                                        typename MeshTraits< MeshConfig >::template EntityTraits< DimensionsTag::value >::EntityTopology >,
      public MeshSuperentityStorageLayers< MeshConfig,
                                           typename MeshTraits< MeshConfig >::template EntityTraits< DimensionsTag::value >::EntityTopology >
 {
@@ -69,23 +85,34 @@ public:
    using MeshTraitsType   = MeshTraits< MeshConfig >;
    using EntityTraitsType = typename MeshTraitsType::template EntityTraits< DimensionsTag::value >;
    using StorageArrayType = typename EntityTraitsType::StorageArrayType;
-   using AccessArrayType  = typename EntityTraitsType::AccessArrayType;
    using GlobalIndexType  = typename EntityTraitsType::GlobalIndexType;
    using EntityType       = typename EntityTraitsType::EntityType;
    using EntityTopology   = typename EntityTraitsType::EntityTopology;
+   using SubentityStorageBaseType = MeshSubentityStorageLayers< MeshConfig, EntityTopology >;
    using SuperentityStorageBaseType = MeshSuperentityStorageLayers< MeshConfig, EntityTopology >;
 
    /****
      * Make visible getters of the lower layer
      */
+   using BaseType::setNumberOfEntities;
    using BaseType::getNumberOfEntities;
    using BaseType::getEntity;
-   using BaseType::getEntities;
 
    MeshStorageLayer()
    {
    }
 
+   bool setNumberOfEntities( DimensionsTag, const GlobalIndexType& entitiesCount )
+   {
+      if( ! this->entities.setSize( entitiesCount ) )
+         return false;
+      if( ! SubentityStorageBaseType::setNumberOfEntities( entitiesCount ) )
+         return false;
+      if( ! SuperentityStorageBaseType::setNumberOfEntities( entitiesCount ) )
+         return false;
+      return true;
+   }
+
    GlobalIndexType getNumberOfEntities( DimensionsTag ) const
    {
       return this->entities.getSize();
@@ -103,19 +130,10 @@ public:
       return this->entities[ entityIndex ];
    }
 
-   AccessArrayType& getEntities( DimensionsTag )
-   {
-      return this->entitiesAccess;
-   }
-
-   const AccessArrayType& getEntities( DimensionsTag ) const
-   {
-      return this->entitiesAccess;
-   }
-
    bool save( File& file ) const
    {
       if( ! BaseType::save( file ) ||
+          ! SubentityStorageBaseType::save( file ) ||
           ! SuperentityStorageBaseType::save( file ) ||
           ! this->entities.save( file ) )
       {
@@ -128,13 +146,13 @@ public:
    bool load( File& file )
    {
       if( ! BaseType::load( file ) ||
+          ! SubentityStorageBaseType::load( file ) ||
           ! SuperentityStorageBaseType::load( file ) ||
           ! this->entities.load( file ) )
       {
          std::cerr << "Loading of the mesh entities with " << DimensionsTag::value << " dimensions failed." << std::endl;
          return false;
       }
-      this->entitiesAccess.bind( this->entities );
       return true;
    }
 
@@ -144,34 +162,21 @@ public:
       str << "The entities with " << DimensionsTag::value << " dimensions are: " << std::endl;
       for( GlobalIndexType i = 0; i < entities.getSize();i ++ )
          str << i << " " << entities[ i ] << std::endl;
+      SubentityStorageBaseType::print( str );
       SuperentityStorageBaseType::print( str );
       str << std::endl;
    }
 
    bool operator==( const MeshStorageLayer& meshLayer ) const
    {
-      return ( BaseType::operator==( meshLayer ) && SuperentityStorageBaseType::operator==( meshLayer ) && entities == meshLayer.entities );
+      return ( BaseType::operator==( meshLayer ) &&
+               SubentityStorageBaseType::operator==( meshLayer ) &&
+               SuperentityStorageBaseType::operator==( meshLayer ) &&
+               entities == meshLayer.entities );
    }
 
 protected:
    StorageArrayType entities;
-
-   AccessArrayType entitiesAccess;
-
-   // Methods for the mesh initializer
-   using BaseType::getEntitiesArray;
-   typename EntityTraitsType::StorageArrayType& getEntitiesArray( DimensionsTag )
-   {
-      return entities;
-   }
-
-   using BaseType::getSuperentityStorageNetwork;
-   template< typename SuperdimensionsTag >
-   typename MeshTraitsType::template SuperentityTraits< EntityTopology, SuperdimensionsTag::value >::StorageNetworkType&
-   getSuperentityStorageNetwork( MeshDimensionsTag< EntityTopology::dimensions > )
-   {
-      return SuperentityStorageBaseType::getSuperentityStorageNetwork( SuperdimensionsTag() );
-   }
 };
 
 template< typename MeshConfig,
@@ -193,7 +198,6 @@ public:
    using MeshTraitsType             = MeshTraits< MeshConfig >;
    using EntityTraitsType           = typename MeshTraitsType::template EntityTraits< 0 >;
    using StorageArrayType           = typename EntityTraitsType::StorageArrayType;
-   using AccessArrayType            = typename EntityTraitsType::AccessArrayType;
    using GlobalIndexType            = typename EntityTraitsType::GlobalIndexType;
    using VertexType                 = typename EntityTraitsType::EntityType;
    using PointType                  = typename VertexType::PointType;
@@ -235,7 +239,16 @@ public:
     * with higher dimensions entities storage layers.
     */
 
-   GlobalIndexType getNumberOfEntities( DimensionTag ) const
+   bool setNumberOfEntities( DimensionsTag, const GlobalIndexType& entitiesCount )
+   {
+      if( ! this->vertices.setSize( entitiesCount ) )
+         return false;
+      if( ! SuperentityStorageBaseType::setNumberOfEntities( entitiesCount ) )
+         return false;
+      return true;
+   }
+
+   GlobalIndexType getNumberOfEntities( DimensionsTag ) const
    {
       return this->vertices.getSize();
    }
@@ -252,16 +265,6 @@ public:
       return this->vertices.getElement( entityIndex );
    }
 
-   AccessArrayType& getEntities( DimensionTag )
-   {
-      return this->verticesAccess;
-   }
-
-   const AccessArrayType& getEntities( DimensionTag ) const
-   {
-      return this->verticesAccess;
-   }
-
    bool save( File& file ) const
    {
       if( ! SuperentityStorageBaseType::save( file ) ||
@@ -281,7 +284,6 @@ public:
          std::cerr << "Loading of the mesh entities with " << DimensionTag::value << " dimensions failed." << std::endl;
          return false;
       }
-      this->verticesAccess.bind( this->vertices );
       return true;
    }
 
@@ -301,21 +303,6 @@ public:
 
 protected:
    StorageArrayType vertices;
-
-   AccessArrayType verticesAccess;
-
-   // Methods for the mesh initializer
-   typename EntityTraitsType::StorageArrayType& getEntitiesArray( DimensionsTag )
-   {
-      return vertices;
-   }
-
-   template< typename SuperdimensionsTag >
-   typename MeshTraitsType::template SuperentityTraits< EntityTopology, SuperdimensionsTag::value >::StorageNetworkType&
-   getSuperentityStorageNetwork( MeshDimensionsTag< EntityTopology::dimensions > )
-   {
-      return SuperentityStorageBaseType::getSuperentityStorageNetwork( SuperdimensionsTag() );
-   }
 };
 
 /****
diff --git a/src/TNL/Meshes/MeshDetails/layers/MeshSubentityAccess.h b/src/TNL/Meshes/MeshDetails/layers/MeshSubentityAccess.h
new file mode 100644
index 0000000000..7c243fbe95
--- /dev/null
+++ b/src/TNL/Meshes/MeshDetails/layers/MeshSubentityAccess.h
@@ -0,0 +1,515 @@
+/***************************************************************************
+                          MeshSubentityAccess.h  -  description
+                             -------------------
+    begin                : Oct 26, 2016
+    copyright            : (C) 2014 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/File.h>
+#include <TNL/Meshes/MeshDimensionsTag.h>
+#include <TNL/Meshes/MeshDetails/traits/MeshTraits.h>
+#include <TNL/Meshes/MeshDetails/traits/MeshSubentityTraits.h>
+#include <TNL/Meshes/MeshDetails/MeshEntityOrientation.h>
+
+namespace TNL {
+namespace Meshes {
+
+template< typename MeshConfig,
+          typename EntityTopology,
+          typename DimensionsTag,
+          bool SubentityStorage =
+               MeshConfig::subentityStorage( EntityTopology(), DimensionsTag::value ),
+          bool SubentityOrientationStorage =
+               MeshConfig::subentityOrientationStorage( EntityTopology(), DimensionsTag::value ) >
+class MeshSubentityAccessLayer;
+
+
+template< typename MeshConfig,
+          typename EntityTopology >
+class MeshSubentityAccess
+   : public MeshSubentityAccessLayer< MeshConfig,
+                                      EntityTopology,
+                                      MeshDimensionsTag< 0 > >
+{
+   using BaseType = MeshSubentityAccessLayer< MeshConfig,
+                                              EntityTopology,
+                                              MeshDimensionsTag< 0 > >;
+
+   using MeshTraitsType = MeshTraits< MeshConfig >;
+
+   template< int Subdimensions >
+   using SubentityTraits = typename MeshTraitsType::template SubentityTraits< EntityTopology, Subdimensions >;
+
+public:
+   template< int Subdimensions >
+   void bindSubentitiesStorageNetwork( const typename SubentityTraits< Subdimensions >::SubentityAccessorType& storage )
+   {
+      static_assert( SubentityTraits< Subdimensions >::storageEnabled, "You try to bind subentities which are not configured for storage." );
+      BaseType::bindSubentitiesStorageNetwork( MeshDimensionsTag< Subdimensions >(),
+                                               storage );
+   }
+
+   template< int Subdimensions >
+   constexpr typename SubentityTraits< Subdimensions >::LocalIndexType getNumberOfSubentities() const
+   {
+      return SubentityTraits< Subdimensions >::count;
+   }
+
+   template< int Subdimensions >
+   void setSubentityIndex( const typename SubentityTraits< Subdimensions >::LocalIndexType& localIndex,
+                           const typename SubentityTraits< Subdimensions >::GlobalIndexType& globalIndex )
+   {
+      static_assert( SubentityTraits< Subdimensions >::storageEnabled, "You try to set subentity which is not configured for storage." );
+      BaseType::setSubentityIndex( MeshDimensionsTag< Subdimensions >(),
+                                   localIndex,
+                                   globalIndex );
+   }
+
+   template< int Subdimensions >
+   typename SubentityTraits< Subdimensions >::GlobalIndexType
+   getSubentityIndex( const typename SubentityTraits< Subdimensions >::LocalIndexType localIndex ) const
+   {
+      static_assert( SubentityTraits< Subdimensions >::storageEnabled, "You try to get subentity which is not configured for storage." );
+      return BaseType::getSubentityIndex( MeshDimensionsTag< Subdimensions >(),
+                                          localIndex );
+   }
+
+   template< int Subdimensions >
+   typename SubentityTraits< Subdimensions >::OrientationArrayType& subentityOrientationsArray()
+   {
+      static_assert( SubentityTraits< Subdimensions >::orientationEnabled, "You try to get subentity orientation which is not configured for storage." );
+      return BaseType::subentityOrientationsArray( MeshDimensionsTag< Subdimensions >() );
+   }
+
+   template< int Subdimensions >
+   typename SubentityTraits< Subdimensions >::IdPermutationArrayType getSubentityOrientation( typename SubentityTraits< Subdimensions >::LocalIndexType index ) const
+   {
+      static_assert( SubentityTraits< Subdimensions >::orientationEnabled, "You try to get subentity orientation which is not configured for storage." );
+      return BaseType::getSubentityOrientation( MeshDimensionsTag< Subdimensions >(), index );
+   }
+
+   bool operator==( const MeshSubentityAccess& other ) const
+   {
+      return BaseType::operator==( other );
+   }
+
+   void print( std::ostream& str ) const
+   {
+      BaseType::print( str );
+   }
+};
+
+
+template< typename MeshConfig,
+          typename EntityTopology,
+          typename DimensionsTag >
+class MeshSubentityAccessLayer< MeshConfig,
+                                EntityTopology,
+                                DimensionsTag,
+                                true,
+                                true >
+   : public MeshSubentityAccessLayer< MeshConfig,
+                                      EntityTopology,
+                                      typename DimensionsTag::Increment >
+{
+   using BaseType = MeshSubentityAccessLayer< MeshConfig,
+                                              EntityTopology,
+                                              typename DimensionsTag::Increment >;
+
+   using MeshTraitsType         = MeshTraits< MeshConfig >;
+   using SubentityTraitsType    = typename MeshTraitsType::template SubentityTraits< EntityTopology, DimensionsTag::value >;
+
+protected:
+   using GlobalIndexType        = typename SubentityTraitsType::GlobalIndexType;
+   using LocalIndexType         = typename SubentityTraitsType::LocalIndexType;
+   using StorageNetworkType     = typename SubentityTraitsType::StorageNetworkType;
+   using SubentityAccessorType  = typename SubentityTraitsType::SubentityAccessorType;
+   using OrientationArrayType   = typename SubentityTraitsType::OrientationArrayType;
+   using IdPermutationArrayType = typename SubentityTraitsType::IdPermutationArrayType;
+
+   MeshSubentityAccessLayer() = default;
+
+   explicit MeshSubentityAccessLayer( const MeshSubentityAccessLayer& layer )
+      : BaseType( layer )
+   {
+      this->subentityIndices.bind( layer.subentityIndices );
+   }
+
+   MeshSubentityAccessLayer& operator=( const MeshSubentityAccessLayer& layer )
+   {
+      BaseType::operator=( layer );
+      this->subentityIndices.bind( layer.subentityIndices );
+      return *this;
+   }
+
+   bool save( File& file ) const
+   {
+      if( ! BaseType::save( file ) )
+      {
+         std::cerr << "Saving of the entity subentities layer with " << DimensionsTag::value << " failed." << std::endl;
+         return false;
+      }
+      return true;
+   }
+
+   bool load( File& file )
+   {
+      if( ! BaseType::load( file ) )
+      {
+         std::cerr << "Loading of the entity subentities layer with " << DimensionsTag::value << " failed." << std::endl;
+         return false;
+      }
+      return true;
+   }
+
+   void print( std::ostream& str ) const
+   {
+      BaseType::print( str );
+      str << "\t Subentities with " << DimensionsTag::value << " dimensions are: " << subentityIndices << "." << std::endl;
+   }
+
+   bool operator==( const MeshSubentityAccessLayer& layer ) const
+   {
+      return ( BaseType::operator==( layer ) &&
+               subentityIndices == layer.subentityIndices );
+   }
+
+   /****
+    * Make visible setters and getters of the lower subentities
+    */
+   using BaseType::bindSubentitiesStorageNetwork;
+   using BaseType::getSubentityIndex;
+   using BaseType::setSubentityIndex;
+   using BaseType::getSubentityIndices;
+
+   /****
+    * Define setter/getter for the current level of the subentities
+    */
+   void bindSubentitiesStorageNetwork( DimensionsTag,
+                                       const SubentityAccessorType& storage )
+   {
+      this->subentityIndices.bind( storage );
+   }
+
+   void setSubentityIndex( DimensionsTag,
+                           const LocalIndexType localIndex,
+                           const GlobalIndexType globalIndex )
+   {
+      this->subentityIndices[ localIndex ] = globalIndex;
+   }
+
+   GlobalIndexType getSubentityIndex( DimensionsTag,
+                                      const LocalIndexType localIndex ) const
+   {
+      return this->subentityIndices[ localIndex ];
+   }
+
+   const SubentityAccessorType& getSubentityIndices( DimensionsTag ) const
+   {
+      return this->subentityIndices;
+   }
+
+   SubentityAccessorType& getSubentityIndices( DimensionsTag )
+   {
+      return this->subentityIndices;
+   }
+
+   using BaseType::getSubentityOrientation;
+   const IdPermutationArrayType& getSubentityOrientation( DimensionsTag, LocalIndexType index) const
+   {
+      Assert( 0 <= index && index < SubentityTraitsType::count, );
+      return this->subentityOrientations[ index ].getSubvertexPermutation();
+   }
+
+   using BaseType::subentityOrientationsArray;
+	OrientationArrayType& subentityOrientationsArray( DimensionsTag ) { return this->subentityOrientations; }
+
+private:
+   SubentityAccessorType subentityIndices;
+
+   OrientationArrayType subentityOrientations;
+};
+
+template< typename MeshConfig,
+          typename EntityTopology,
+          typename DimensionsTag >
+class MeshSubentityAccessLayer< MeshConfig,
+                                EntityTopology,
+                                DimensionsTag,
+                                true,
+                                false >
+   : public MeshSubentityAccessLayer< MeshConfig,
+                                      EntityTopology,
+                                      typename DimensionsTag::Increment >
+{
+   static_assert( DimensionsTag::value < EntityTopology::dimensions, "" );
+   using BaseType = MeshSubentityAccessLayer< MeshConfig,
+                                              EntityTopology,
+                                              typename DimensionsTag::Increment >;
+
+   using MeshTraitsType        = MeshTraits< MeshConfig >;
+   using SubentityTraitsType   = typename MeshTraitsType::template SubentityTraits< EntityTopology, DimensionsTag::value >;
+
+protected:
+   using GlobalIndexType       = typename SubentityTraitsType::GlobalIndexType;
+   using LocalIndexType        = typename SubentityTraitsType::LocalIndexType;
+   using StorageNetworkType    = typename SubentityTraitsType::StorageNetworkType;
+   using SubentityAccessorType = typename SubentityTraitsType::SubentityAccessorType;
+
+   MeshSubentityAccessLayer() = default;
+
+   explicit MeshSubentityAccessLayer( const MeshSubentityAccessLayer& layer )
+      : BaseType( layer )
+   {
+      this->subentityIndices.bind( layer.subentityIndices );
+   }
+
+   MeshSubentityAccessLayer& operator=( const MeshSubentityAccessLayer& layer )
+   {
+      BaseType::operator=( layer );
+      this->subentityIndices.bind( layer.subentityIndices );
+      return *this;
+   }
+
+   bool save( File& file ) const
+   {
+      if( ! BaseType::save( file ) )
+      {
+         std::cerr << "Saving of the entity subentities layer with " << DimensionsTag::value << " failed." << std::endl;
+         return false;
+      }
+      return true;
+   }
+
+   bool load( File& file )
+   {
+      if( ! BaseType::load( file ) )
+      {
+         std::cerr << "Loading of the entity subentities layer with " << DimensionsTag::value << " failed." << std::endl;
+         return false;
+      }
+      return true;
+   }
+
+   void print( std::ostream& str ) const
+   {
+      BaseType::print( str );
+      str << "\t Subentities with " << DimensionsTag::value << " dimensions are: " << subentityIndices << "." << std::endl;
+   }
+
+   bool operator==( const MeshSubentityAccessLayer& layer ) const
+   {
+      return ( BaseType::operator==( layer ) &&
+               subentityIndices == layer.subentityIndices );
+   }
+
+   /****
+    * Make visible setters and getters of the lower subentities
+    */
+   using BaseType::bindSubentitiesStorageNetwork;
+   using BaseType::getSubentityIndex;
+   using BaseType::setSubentityIndex;
+   using BaseType::getSubentityIndices;
+
+   /****
+    * Define setter/getter for the current level of the subentities
+    */
+   void bindSubentitiesStorageNetwork( DimensionsTag,
+                                       const SubentityAccessorType& storage )
+   {
+      this->subentityIndices.bind( storage );
+   }
+
+   void setSubentityIndex( DimensionsTag,
+                           const LocalIndexType localIndex,
+                           const GlobalIndexType globalIndex )
+   {
+      this->subentityIndices[ localIndex ] = globalIndex;
+   }
+
+   GlobalIndexType getSubentityIndex( DimensionsTag,
+                                      const LocalIndexType localIndex ) const
+   {
+      return this->subentityIndices[ localIndex ];
+   }
+
+   const SubentityAccessorType& getSubentityIndices( DimensionsTag ) const
+   {
+      return this->subentityIndices;
+   }
+
+   SubentityAccessorType& getSubentityIndices( DimensionsTag )
+   {
+      return this->subentityIndices;
+   }
+
+private:
+   SubentityAccessorType subentityIndices;
+};
+
+
+template< typename MeshConfig,
+          typename EntityTopology >
+class MeshSubentityAccessLayer< MeshConfig,
+                                EntityTopology,
+                                MeshDimensionsTag< EntityTopology::dimensions >,
+                                true,
+                                true >
+{
+   using DimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
+
+protected:
+   /***
+    * Necessary because of 'using BaseType::...;' in the derived classes
+    */
+   template< typename SubentityAccessorType >
+   void bindSubentitiesStorageNetwork( DimensionsTag,
+                                       const SubentityAccessorType& storage ) {}
+   void getNumberOfSubentities( DimensionsTag ) const {}
+   template< typename LocalIndexType >
+   void getSubentityIndex( DimensionsTag,
+                           const LocalIndexType localIndex ) const {}
+   template< typename LocalIndexType, typename GlobalIndexType >
+   void setSubentityIndex( DimensionsTag,
+                           const LocalIndexType& localIndex,
+                           const GlobalIndexType& globalIndex ) {}
+   void getSubentityIndices() {}
+
+   template< typename LocalIndexType >
+   void getSubentityOrientation( DimensionsTag, LocalIndexType index) const {}
+	void subentityOrientationsArray( DimensionsTag ) {}
+
+   bool save( File& file ) const
+   {
+      return true;
+   }
+
+   bool load( File& file )
+   {
+      return true;
+   }
+
+   bool operator==( const MeshSubentityAccessLayer& other ) const
+   {
+      return true;
+   }
+
+   void print( std::ostream& str ) const {}
+};
+
+template< typename MeshConfig,
+          typename EntityTopology >
+class MeshSubentityAccessLayer< MeshConfig,
+                                EntityTopology,
+                                MeshDimensionsTag< EntityTopology::dimensions >,
+                                true,
+                                false >
+{
+   using DimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
+
+protected:
+   /***
+    * Necessary because of 'using BaseType::...;' in the derived classes
+    */
+   template< typename SubentityAccessorType >
+   void bindSubentitiesStorageNetwork( DimensionsTag,
+                                       const SubentityAccessorType& storage ) {}
+   void getNumberOfSubentities( DimensionsTag ) const {}
+   template< typename LocalIndexType >
+   void getSubentityIndex( DimensionsTag,
+                           const LocalIndexType localIndex ) const {}
+   template< typename LocalIndexType, typename GlobalIndexType >
+   void setSubentityIndex( DimensionsTag,
+                           const LocalIndexType& localIndex,
+                           const GlobalIndexType& globalIndex ) {}
+   void getSubentityIndices() {}
+
+   template< typename LocalIndexType >
+   void getSubentityOrientation( DimensionsTag, LocalIndexType index) const {}
+	void subentityOrientationsArray( DimensionsTag ) {}
+
+   bool save( File& file ) const
+   {
+      return true;
+   }
+
+   bool load( File& file )
+   {
+      return true;
+   }
+
+   bool operator==( const MeshSubentityAccessLayer& other ) const
+   {
+      return true;
+   }
+
+   void print( std::ostream& str ) const {}
+};
+
+template< typename MeshConfig,
+          typename EntityTopology >
+class MeshSubentityAccessLayer< MeshConfig,
+                                EntityTopology,
+                                MeshDimensionsTag< EntityTopology::dimensions >,
+                                false,
+                                true >
+{
+   using DimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
+
+protected:
+   /***
+    * Necessary because of 'using BaseType::...;' in the derived classes
+    */
+   template< typename SubentityAccessorType >
+   void bindSubentitiesStorageNetwork( DimensionsTag,
+                                       const SubentityAccessorType& storage ) {}
+   void getNumberOfSubentities( DimensionsTag ) const {}
+   template< typename LocalIndexType >
+   void getSubentityIndex( DimensionsTag,
+                           const LocalIndexType localIndex ) const {}
+   template< typename LocalIndexType, typename GlobalIndexType >
+   void setSubentityIndex( DimensionsTag,
+                           const LocalIndexType& localIndex,
+                           const GlobalIndexType& globalIndex ) {}
+   void getSubentityIndices() {}
+
+   template< typename LocalIndexType >
+   void getSubentityOrientation( DimensionsTag, LocalIndexType index) const {}
+	void subentityOrientationsArray( DimensionsTag ) {}
+
+   bool save( File& file ) const
+   {
+      return true;
+   }
+
+   bool load( File& file )
+   {
+      return true;
+   }
+
+   bool operator==( const MeshSubentityAccessLayer& other ) const
+   {
+      return true;
+   }
+
+   void print( std::ostream& str ) const {}
+};
+
+template< typename MeshConfig,
+          typename EntityTopology,
+          typename DimensionsTag >
+class MeshSubentityAccessLayer< MeshConfig,
+                                EntityTopology,
+                                DimensionsTag,
+                                false,
+                                false >
+{
+};
+
+} // namespace Meshes
+} // namespace TNL
diff --git a/src/TNL/Meshes/MeshDetails/layers/MeshSubentityStorageLayer.h b/src/TNL/Meshes/MeshDetails/layers/MeshSubentityStorageLayer.h
index b9551a50d8..e0b71e011a 100644
--- a/src/TNL/Meshes/MeshDetails/layers/MeshSubentityStorageLayer.h
+++ b/src/TNL/Meshes/MeshDetails/layers/MeshSubentityStorageLayer.h
@@ -2,7 +2,7 @@
                           MeshSubentityStorageLayer.h  -  description
                              -------------------
     begin                : Feb 11, 2014
-    copyright            : (C) 2014 by Tomas Oberhuber
+    copyright            : (C) 2014 by Tomas Oberhuber et al.
     email                : tomas.oberhuber@fjfi.cvut.cz
  ***************************************************************************/
 
@@ -17,126 +17,79 @@
 #pragma once
 
 #include <TNL/File.h>
-#include <TNL/Meshes/MeshDimensionTag.h>
+#include <TNL/Meshes/MeshDimensionsTag.h>
+#include <TNL/Meshes/MeshDetails/traits/MeshTraits.h>
 #include <TNL/Meshes/MeshDetails/traits/MeshSubentityTraits.h>
-#include <TNL/Meshes/MeshDetails/MeshEntityOrientation.h>
 
 namespace TNL {
 namespace Meshes {
 
 template< typename MeshConfig,
           typename EntityTopology,
-          typename DimensionTag,
+          typename SubdimensionsTag,
           bool SubentityStorage =
-               MeshTraits< MeshConfig >::template SubentityTraits< EntityTopology, DimensionsTag::value >::storageEnabled,
-          bool SubentityOrientationStorage =
-               MeshTraits< MeshConfig >::template SubentityTraits< EntityTopology, DimensionsTag::value >::orientationEnabled >
+               MeshConfig::subentityStorage( EntityTopology(), SubdimensionsTag::value ) >
 class MeshSubentityStorageLayer;
 
-
 template< typename MeshConfig,
           typename EntityTopology >
 class MeshSubentityStorageLayers
    : public MeshSubentityStorageLayer< MeshConfig,
                                        EntityTopology,
-                                       MeshDimensionsTag< EntityTopology::dimensions - 1 > >
+                                       MeshDimensionsTag< 0 > >
 {
    using BaseType = MeshSubentityStorageLayer< MeshConfig,
                                                EntityTopology,
-                                               MeshDimensionsTag< EntityTopology::dimensions - 1 > >;
-
-   static constexpr int Dimensions = EntityTopology::dimensions;
+                                               MeshDimensionsTag< 0 > >;
    using MeshTraitsType = MeshTraits< MeshConfig >;
 
-   template< int Subdimensions >
-   using SubentityTraits = typename MeshTraitsType::template SubentityTraits< EntityTopology, Subdimensions >;
-
 public:
    template< int Subdimensions >
-   constexpr typename SubentityTraits< Subdimensions >::LocalIndexType getNumberOfSubentities() const
-   {
-      return SubentityTraits< Subdimensions >::count;
-   }
-
-   template< int Subdimensions >
-   void setSubentityIndex( const typename SubentityTraits< Subdimensions >::LocalIndexType& localIndex,
-                           const typename SubentityTraits< Subdimensions >::GlobalIndexType& globalIndex )
+   typename MeshTraitsType::template SubentityTraits< EntityTopology, Subdimensions >::StorageNetworkType&
+   getSubentityStorageNetwork()
    {
-      static_assert( SubentityTraits< Subdimensions >::storageEnabled, "You try to set subentity which is not configured for storage." );
-      BaseType::setSubentityIndex( MeshDimensionsTag< Subdimensions >(),
-                                   localIndex,
-                                   globalIndex );
-   }
-
-   template< int Subdimensions >
-   typename SubentityTraits< Subdimensions >::GlobalIndexType
-   getSubentityIndex( const typename SubentityTraits< Subdimensions >::LocalIndexType localIndex ) const
-   {
-      static_assert( SubentityTraits< Subdimensions >::storageEnabled, "You try to get subentity which is not configured for storage." );
-      return BaseType::getSubentityIndex( MeshDimensionsTag< Subdimensions >(),
-                                          localIndex );
-   }
-
-   template< int Subdimensions >
-   typename SubentityTraits< Subdimensions >::OrientationArrayType& subentityOrientationsArray()
-   {
-      static_assert( SubentityTraits< Subdimensions >::orientationEnabled, "You try to get subentity orientation which is not configured for storage." );
-      return BaseType::subentityOrientationsArray( MeshDimensionsTag< Subdimensions >() );
-   }
-
-   template< int Subdimensions >
-   typename SubentityTraits< Subdimensions >::IdPermutationArrayType getSubentityOrientation( typename SubentityTraits< Subdimensions >::LocalIndexType index ) const
-   {
-      static_assert( SubentityTraits< Subdimensions >::orientationEnabled, "You try to get subentity orientation which is not configured for storage." );
-      return BaseType::getSubentityOrientation( MeshDimensionsTag< Subdimensions >(), index );
+      static_assert( EntityTopology::dimensions > Subdimensions, "Invalid combination of Dimensions and Subdimensions." );
+      return BaseType::getSubentityStorageNetwork( MeshDimensionsTag< EntityTopology::dimensions >(), MeshDimensionsTag< Subdimensions >() );
    }
 };
 
-
 template< typename MeshConfig,
           typename EntityTopology,
-          typename DimensionTag >
+          typename SubdimensionsTag >
 class MeshSubentityStorageLayer< MeshConfig,
                                  EntityTopology,
-                                 DimensionsTag,
-                                 true,
+                                 SubdimensionsTag,
                                  true >
    : public MeshSubentityStorageLayer< MeshConfig,
                                        EntityTopology,
-                                       typename DimensionsTag::Decrement >
+                                       typename SubdimensionsTag::Increment >
 {
    using BaseType = MeshSubentityStorageLayer< MeshConfig,
                                                EntityTopology,
-                                               typename DimensionsTag::Decrement >;
+                                               typename SubdimensionsTag::Increment >;
+
+   using MeshTraitsType      = MeshTraits< MeshConfig >;
+   using SubentityTraitsType = typename MeshTraitsType::template SubentityTraits< EntityTopology, SubdimensionsTag::value >;
 
 protected:
-   static constexpr int Dimensions = DimensionsTag::value;
-   using MeshTraitsType         = MeshTraits< MeshConfig >;
-   using SubentityTraitsType    = typename MeshTraitsType::template SubentityTraits< EntityTopology, Dimensions >;
-   using GlobalIndexType        = typename MeshTraitsType::GlobalIndexType;
-   using LocalIndexType         = typename MeshTraitsType::LocalIndexType;
-   using IdArrayType            = typename SubentityTraitsType::IdArrayType;
-   using OrientationArrayType   = typename SubentityTraitsType::OrientationArrayType;
-   using IdPermutationArrayType = typename SubentityTraitsType::IdPermutationArrayType;
-
-   MeshSubentityStorageLayer()
-   {
-      this->subentitiesIndices.setValue( -1 );
-   }
+   using GlobalIndexType    = typename SubentityTraitsType::GlobalIndexType;
+   using LocalIndexType     = typename SubentityTraitsType::LocalIndexType;
+   using StorageNetworkType = typename SubentityTraitsType::StorageNetworkType;
 
-   MeshSubentityStorageLayer& operator=( const MeshSubentityStorageLayer& layer )
+   bool setNumberOfEntities( const GlobalIndexType& entitiesCount )
    {
-      BaseType::operator=( layer );
-      this->subentitiesIndices = layer.subentitiesIndices;
-      return *this;
+      if( ! BaseType::setNumberOfEntities( entitiesCount ) )
+         return false;
+      this->storageNetwork.setKeysRange( entitiesCount );
+      return this->storageNetwork.allocate();
    }
 
    bool save( File& file ) const
    {
       if( ! BaseType::save( file ) ||
-          ! this->subentitiesIndices.save( file ) )
+          ! this->storageNetwork.save( file ) )
       {
-         std::cerr << "Saving of the entity subentities layer with " << DimensionTag::value << " failed." << std::endl;
+         std::cerr << "Saving of the entity subentities layer with " << SubdimensionsTag::value << " dimensions failed." << std::endl;
          return false;
       }
       return true;
@@ -145,9 +98,9 @@ protected:
    bool load( File& file )
    {
       if( ! BaseType::load( file ) ||
-          ! this->subentitiesIndices.load( file ) )
+          ! this->storageNetwork.load( file ) )
       {
-         std::cerr << "Loading of the entity subentities layer with " << DimensionTag::value << " failed." << std::endl;
+         std::cerr << "Loading of the entity subentities layer with " << SubdimensionsTag::value << " dimensions failed." << std::endl;
          return false;
       }
       return true;
@@ -156,258 +109,105 @@ protected:
    void print( std::ostream& str ) const
    {
       BaseType::print( str );
-      str << "\t Subentities with " << DimensionsTag::value << " dimensions are: " << subentitiesIndices << "." << std::endl;
+      str << "Storage network for subentities with " << SubdimensionsTag::value << " dimensions of entities with " << EntityTopology::dimensions << " dimensions is: " << std::endl;
+      str << this->storageNetwork << std::endl;
    }
 
    bool operator==( const MeshSubentityStorageLayer& layer ) const
    {
       return ( BaseType::operator==( layer ) &&
-               subentitiesIndices == layer.subentitiesIndices );
-   }
-
-   /****
-    * Make visible setters and getters of the lower subentities
-    */
-   using BaseType::getSubentityIndex;
-   using BaseType::setSubentityIndex;
-
-   /****
-    * Define setter/getter for the current level of the subentities
-    */
-   void setSubentityIndex( DimensionTag,
-                           const LocalIndexType localIndex,
-                           const GlobalIndexType globalIndex )
-   {
-      this->subentitiesIndices[ localIndex ] = globalIndex;
-   }
-
-   GlobalIndexType getSubentityIndex( DimensionTag,
-                                      const LocalIndexType localIndex ) const
-   {
-      return this->subentitiesIndices[ localIndex ];
+               storageNetwork == layer.storageNetwork );
    }
 
-   using BaseType::getSubentityOrientation;
-   const IdPermutationArrayType& getSubentityOrientation( DimensionsTag, LocalIndexType index) const
+   using BaseType::getSubentityStorageNetwork;
+   StorageNetworkType& getSubentityStorageNetwork( MeshDimensionsTag< EntityTopology::dimensions >, SubdimensionsTag )
    {
-      Assert( 0 <= index && index < SubentityTraitsType::count, );
-      return this->subentityOrientations[ index ].getSubvertexPermutation();
+      return this->storageNetwork;
    }
 
-   using BaseType::subentityOrientationsArray;
-	OrientationArrayType& subentityOrientationsArray( DimensionsTag ) { return this->subentityOrientations; }
-
 private:
-   IdArrayType subentitiesIndices;
-
-   OrientationArrayType subentityOrientations;
+   StorageNetworkType storageNetwork;
 };
 
-
 template< typename MeshConfig,
           typename EntityTopology,
-          typename DimensionTag >
+          typename SubdimensionsTag >
 class MeshSubentityStorageLayer< MeshConfig,
                                  EntityTopology,
-                                 DimensionsTag,
-                                 true,
+                                 SubdimensionsTag,
                                  false >
    : public MeshSubentityStorageLayer< MeshConfig,
                                        EntityTopology,
-                                       typename DimensionsTag::Decrement >
+                                       typename SubdimensionsTag::Increment >
 {
-   using BaseType = MeshSubentityStorageLayer< MeshConfig,
-                                               EntityTopology,
-                                               typename DimensionsTag::Decrement >;
-
-protected:
-   static constexpr int Dimensions = DimensionsTag::value;
-   using MeshTraitsType      = MeshTraits< MeshConfig >;
-   using SubentityTraitsType = typename MeshTraitsType::template SubentityTraits< EntityTopology, Dimensions >;
-   using GlobalIndexType     = typename MeshTraitsType::GlobalIndexType;
-   using LocalIndexType      = typename MeshTraitsType::LocalIndexType;
-   using IdArrayType         = typename SubentityTraitsType::IdArrayType;
+};
 
-   MeshSubentityStorageLayer()
-   {
-      this->subentitiesIndices.setValue( -1 );
-   }
 
-   MeshSubentityStorageLayer& operator=( const MeshSubentityStorageLayer& layer )
-   {
-      BaseType::operator=( layer );
-      this->subentitiesIndices = layer.subentitiesIndices;
-      return *this;
-   }
+template< typename MeshConfig,
+          typename EntityTopology >
+class MeshSubentityStorageLayer< MeshConfig,
+                                 EntityTopology,
+                                 MeshDimensionsTag< EntityTopology::dimensions >,
+                                 true >
+{
+   using SubdimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
 
-   bool save( File& file ) const
+protected:
+   /****
+    * These methods are due to 'using BaseType::...;' in the derived classes.
+    */
+   template< typename GlobalIndexType >
+   bool setNumberOfEntities( const GlobalIndexType& entitiesCount )
    {
-      if( ! BaseType::save( file ) ||
-          ! this->subentitiesIndices.save( file ) )
-      {
-         std::cerr << "Saving of the entity subentities layer with " << DimensionTag::value << " failed." << std::endl;
-         return false;
-      }
       return true;
    }
 
-   bool load( File& file )
-   {
-      if( ! BaseType::load( file ) ||
-          ! this->subentitiesIndices.load( file ) )
-      {
-         std::cerr << "Loading of the entity subentities layer with " << DimensionTag::value << " failed." << std::endl;
-         return false;
-      }
-      return true;
-   }
-
-   void print( std::ostream& str ) const
-   {
-      BaseType::print( str );
-      str << "\t Subentities with " << DimensionsTag::value << " dimensions are: " << subentitiesIndices << "." << std::endl;
-   }
+   void print( std::ostream& str ) const {}
 
    bool operator==( const MeshSubentityStorageLayer& layer ) const
    {
-      return ( BaseType::operator==( layer ) &&
-               subentitiesIndices == layer.subentitiesIndices );
+      return true;
    }
 
-   /****
-    * Make visible setters and getters of the lower subentities
-    */
-   using BaseType::getSubentityIndex;
-   using BaseType::setSubentityIndex;
-
-   /****
-    * Define setter/getter for the current level of the subentities
-    */
-   void setSubentityIndex( DimensionTag,
-                           const LocalIndexType localIndex,
-                           const GlobalIndexType globalIndex )
+   bool save( File& file ) const
    {
-      this->subentitiesIndices[ localIndex ] = globalIndex;
+      return true;
    }
 
-   GlobalIndexType getSubentityIndex( DimensionTag,
-                                      const LocalIndexType localIndex ) const
+   bool load( File& file )
    {
-      return this->subentitiesIndices[ localIndex ];
+      return true;
    }
-
-   using BaseType::subentityOrientationsArray;
-   void subentityOrientationsArray() {}
-
-private:
-   IdArrayType subentitiesIndices;
+ 
+   void getSubentityStorageNetwork( MeshDimensionsTag< EntityTopology::dimensions >, SubdimensionsTag ) {}
 };
 
-template< typename MeshConfig,
-          typename EntityTopology,
-          typename DimensionTag >
-class MeshSubentityStorageLayer< MeshConfig,
-                                 EntityTopology,
-                                 DimensionsTag,
-                                 false,
-                                 false >
-   : public MeshSubentityStorageLayer< MeshConfig,
-                                       EntityTopology,
-                                       typename DimensionsTag::Decrement >
-{
-};
-
-
 template< typename MeshConfig,
           typename EntityTopology >
 class MeshSubentityStorageLayer< MeshConfig,
                                  EntityTopology,
-                                 MeshDimensionsTag< 0 >,
-                                 true,
+                                 MeshDimensionsTag< EntityTopology::dimensions >,
                                  false >
 {
-   using DimensionsTag = MeshDimensionsTag< 0 >;
+   using SubdimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
 
 protected:
-   static constexpr int Dimensions = 0;
-   using MeshTraitsType      = MeshTraits< MeshConfig >;
-   using SubentityTraitsType = typename MeshTraitsType::template SubentityTraits< EntityTopology, Dimensions >;
-   using GlobalIndexType     = typename MeshTraitsType::GlobalIndexType;
-   using LocalIndexType      = typename MeshTraitsType::LocalIndexType;
-   using IdArrayType         = typename SubentityTraitsType::IdArrayType;
-
-   MeshSubentityStorageLayer()
-   {
-      this->verticesIndices.setValue( -1 );
-   }
-
-   MeshSubentityStorageLayer& operator=( const MeshSubentityStorageLayer& layer )
-   {
-      this->verticesIndices = layer.verticesIndices;
-      return *this;
-   }
-
-   bool save( File& file ) const
-   {
-      if( ! this->verticesIndices.save( file ) )
-      {
-         std::cerr << "Saving of the entity subentities layer with " << DimensionTag::value << " failed." << std::endl;
-         return false;
-      }
-      return true;
-   }
-
-   bool load( File& file )
+   /****
+    * These methods are due to 'using BaseType::...;' in the derived classes.
+    */
+   template< typename GlobalIndexType >
+   bool setNumberOfEntities( const GlobalIndexType& entitiesCount )
    {
-      if( ! this->verticesIndices.load( file ) )
-      {
-         std::cerr << "Loading of the entity subentities layer with " << DimensionTag::value << " failed." << std::endl;
-         return false;
-      }
       return true;
    }
 
-   void print( std::ostream& str ) const
-   {
-      str << "\t Subentities with " << DimensionsTag::value << " dimensions are: " << this->verticesIndices << "." << std::endl;
-   }
+   void print( std::ostream& str ) const {}
 
    bool operator==( const MeshSubentityStorageLayer& layer ) const
    {
-      return ( verticesIndices == layer.verticesIndices );
-   }
-
-   GlobalIndexType getSubentityIndex( DimensionTag,
-                                      const LocalIndexType localIndex ) const
-   {
-      return this->verticesIndices[ localIndex ];
-   }
-   void setSubentityIndex( DimensionTag,
-                           const LocalIndexType localIndex,
-                           const GlobalIndexType globalIndex )
-   {
-      this->verticesIndices[ localIndex ] = globalIndex;
+      return true;
    }
 
-protected:
-   /***
-    *  Necessary because of 'using BaseType::...;' in the derived classes
-    */
-   void getSubentityOrientation() {}
-   void subentityOrientationsArray() {}
-
-   IdArrayType verticesIndices;
-};
-
-template< typename MeshConfig,
-          typename EntityTopology >
-class MeshSubentityStorageLayer< MeshConfig,
-                                 EntityTopology,
-                                 MeshDimensionsTag< 0 >,
-                                 false,
-                                 false >
-{
-public:
    bool save( File& file ) const
    {
       return true;
@@ -417,6 +217,8 @@ public:
    {
       return true;
    }
+ 
+   void getSubentityStorageNetwork( MeshDimensionsTag< EntityTopology::dimensions >, SubdimensionsTag ) {}
 };
 
 } // namespace Meshes
diff --git a/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityAccess.h b/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityAccess.h
index ca1463fbed..e352da7f07 100644
--- a/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityAccess.h
+++ b/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityAccess.h
@@ -145,6 +145,7 @@ public:
    {
       BaseType::operator=( layer );
       this->superentityIndices.bind( layer.superentityIndices );
+      return *this;
    }
 
    /****
@@ -223,7 +224,6 @@ class MeshSuperentityAccessLayer< MeshConfig,
                                   MeshDimensionsTag< EntityTopology::dimensions >,
                                   false >
 {
-   static constexpr int Dimensions = EntityTopology::dimensions;
    using DimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
    using MeshTraitsType = MeshTraits< MeshConfig >;
    using SuperentityTraitsType = typename MeshTraitsType::template SuperentityTraits< EntityTopology, DimensionsTag::value >;
@@ -246,7 +246,6 @@ protected:
    void setSuperentityIndex( DimensionsTag,
                              const LocalIndexType& localIndex,
                              const GlobalIndexType& globalIndex ) {}
-
    void getSuperentityIndices() {}
 
    bool operator==( const MeshSuperentityAccess< MeshConfig, EntityTopology >& other ) const
@@ -264,7 +263,6 @@ class MeshSuperentityAccessLayer< MeshConfig,
                                   MeshDimensionsTag< EntityTopology::dimensions >,
                                   true >
 {
-   static constexpr int Dimensions = EntityTopology::dimensions;
    using DimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
    using MeshTraitsType = MeshTraits< MeshConfig >;
    using SuperentityTraitsType = typename MeshTraitsType::template SuperentityTraits< EntityTopology, DimensionsTag::value >;
@@ -287,7 +285,6 @@ protected:
    void setSuperentityIndex( DimensionsTag,
                              const LocalIndexType& localIndex,
                              const GlobalIndexType& globalIndex ) {}
-
    void getSuperentityIndices() {}
 
    bool operator==( const MeshSuperentityAccessLayer& other ) const
diff --git a/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageLayer.h b/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageLayer.h
index 58b4c264e1..8f7cfe1e36 100644
--- a/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageLayer.h
+++ b/src/TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageLayer.h
@@ -26,9 +26,9 @@ namespace Meshes {
 
 template< typename MeshConfig,
           typename EntityTopology,
-          typename DimensionTag,
+          typename SuperdimensionsTag,
           bool SuperentityStorage =
-               MeshSuperentityTraits< MeshConfig, EntityTopology, DimensionsTag::value >::storageEnabled >
+               MeshTraits< MeshConfig >::template SuperentityTraits< EntityTopology, SuperdimensionsTag::value >::storageEnabled >
 class MeshSuperentityStorageLayer;
 
 template< typename MeshConfig,
@@ -38,19 +38,30 @@ class MeshSuperentityStorageLayers
                                          EntityTopology,
                                          MeshDimensionsTag< MeshTraits< MeshConfig >::meshDimensions > >
 {
+   using BaseType = MeshSuperentityStorageLayer< MeshConfig,
+                                                 EntityTopology,
+                                                 MeshDimensionsTag< MeshTraits< MeshConfig >::meshDimensions > >;
+   using MeshTraitsType = MeshTraits< MeshConfig >;
+
+public:
+   template< int Superdimensions >
+   typename MeshTraitsType::template SuperentityTraits< EntityTopology, Superdimensions >::StorageNetworkType&
+   getSuperentityStorageNetwork( MeshDimensionsTag< EntityTopology::dimensions > )
+   {
+      return BaseType::getSuperentityStorageNetwork( MeshDimensionsTag< Superdimensions >() );
+   }
 };
 
 template< typename MeshConfig,
           typename EntityTopology,
-          typename DimensionTag >
-class MeshSuperentityStorageLayer< MeshConfig, EntityTopology, DimensionTag, true >
-   : public MeshSuperentityStorageLayer< MeshConfig, EntityTopology, typename DimensionTag::Decrement >
+          typename SuperdimensionsTag >
+class MeshSuperentityStorageLayer< MeshConfig, EntityTopology, SuperdimensionsTag, true >
+   : public MeshSuperentityStorageLayer< MeshConfig, EntityTopology, typename SuperdimensionsTag::Decrement >
 {
-   using BaseType = MeshSuperentityStorageLayer< MeshConfig, EntityTopology, typename DimensionsTag::Decrement >;
+   using BaseType = MeshSuperentityStorageLayer< MeshConfig, EntityTopology, typename SuperdimensionsTag::Decrement >;
 
-   static constexpr int Dimensions = DimensionsTag::value;
    using MeshTraitsType        = MeshTraits< MeshConfig >;
-   using SuperentityTraitsType = typename MeshTraitsType::template SuperentityTraits< EntityTopology, Dimensions >;
+   using SuperentityTraitsType = typename MeshTraitsType::template SuperentityTraits< EntityTopology, SuperdimensionsTag::value >;
 
 protected:
    using GlobalIndexType    = typename SuperentityTraitsType::GlobalIndexType;
@@ -59,12 +70,20 @@ protected:
  
    MeshSuperentityStorageLayer& operator=( const MeshSuperentityStorageLayer& layer ) = delete;
 
+   bool setNumberOfEntities( const GlobalIndexType& entitiesCount )
+   {
+      if( ! BaseType::setNumberOfEntities( entitiesCount ) )
+         return false;
+      this->storageNetwork.setKeysRange( entitiesCount );
+      return true;
+   }
+
    bool save( File& file ) const
    {
       if( ! BaseType::save( file ) ||
           ! this->storageNetwork.save( file ) )
       {
-         std::cerr << "Saving of the entity superentities layer with " << DimensionsTag::value << " dimensions failed." << std::endl;
+         std::cerr << "Saving of the entity superentities layer with " << SuperdimensionsTag::value << " dimensions failed." << std::endl;
          return false;
       }
       return true;
@@ -75,7 +94,7 @@ protected:
       if( ! BaseType::load( file ) ||
           ! this->storageNetwork.load( file ) )
       {
-         std::cerr << "Loading of the entity superentities layer with " << DimensionsTag::value << " dimensions failed." << std::endl;
+         std::cerr << "Loading of the entity superentities layer with " << SuperdimensionsTag::value << " dimensions failed." << std::endl;
          return false;
       }
       return true;
@@ -84,7 +103,7 @@ protected:
    void print( std::ostream& str ) const
    {
       BaseType::print( str );
-      str << "Storage network for superentities with " << DimensionsTag::value << " dimensions of entities with " << EntityTopology::dimensions << " dimensions is: " << std::endl;
+      str << "Storage network for superentities with " << SuperdimensionsTag::value << " dimensions of entities with " << EntityTopology::dimensions << " dimensions is: " << std::endl;
       str << this->storageNetwork << std::endl;
    }
 
@@ -94,23 +113,21 @@ protected:
                storageNetwork == layer.storageNetwork );
    }
 
+   using BaseType::getSuperentityStorageNetwork;
+   StorageNetworkType& getSuperentityStorageNetwork( SuperdimensionsTag )
+   {
+      return this->storageNetwork;
+   }
+
 private:
-    StorageNetworkType storageNetwork;
- 
-   // TODO: this is only for the mesh initializer - fix it
-   public:
-      using BaseType::getSuperentityStorageNetwork;
-      StorageNetworkType& getSuperentityStorageNetwork( DimensionsTag )
-      {
-         return this->storageNetwork;
-      }
+   StorageNetworkType storageNetwork;
 };
 
 template< typename MeshConfig,
           typename EntityTopology,
-          typename DimensionTag >
-class MeshSuperentityStorageLayer< MeshConfig, EntityTopology, DimensionTag, false >
-   : public MeshSuperentityStorageLayer< MeshConfig, EntityTopology, typename DimensionTag::Decrement >
+          typename SuperdimensionsTag >
+class MeshSuperentityStorageLayer< MeshConfig, EntityTopology, SuperdimensionsTag, false >
+   : public MeshSuperentityStorageLayer< MeshConfig, EntityTopology, typename SuperdimensionsTag::Decrement >
 {
 };
 
@@ -118,11 +135,10 @@ template< typename MeshConfig,
           typename EntityTopology >
 class MeshSuperentityStorageLayer< MeshConfig, EntityTopology, MeshDimensionTag< EntityTopology::dimensions >, false >
 {
-   using DimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
-   static constexpr int Dimensions = DimensionsTag::value;
+   using SuperdimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
    using MeshTraitsType = MeshTraits< MeshConfig >;
-   using SuperentityTraitsType = typename MeshTraitsType::template SuperentityTraits< EntityTopology, Dimensions >;
-   using ThisType = MeshSuperentityStorageLayer< MeshConfig, EntityTopology, DimensionsTag, false >;
+   using SuperentityTraitsType = typename MeshTraitsType::template SuperentityTraits< EntityTopology, SuperdimensionsTag::value >;
+   using ThisType = MeshSuperentityStorageLayer< MeshConfig, EntityTopology, SuperdimensionsTag, false >;
 
 protected:
    using GlobalIndexType    = typename SuperentityTraitsType::GlobalIndexType;
@@ -132,6 +148,11 @@ protected:
    /****
     * These methods are due to 'using BaseType::...;' in the derived classes.
     */
+   bool setNumberOfEntities( const GlobalIndexType& entitiesCount )
+   {
+      return true;
+   }
+
    void print( std::ostream& str ) const {}
 
    bool operator==( const ThisType& layer ) const
@@ -149,7 +170,7 @@ protected:
       return true;
    }
  
-   void getSuperentityStorageNetwork( DimensionsTag ) {}
+   void getSuperentityStorageNetwork( SuperdimensionsTag ) {}
 };
 
 template< typename MeshConfig,
@@ -159,11 +180,10 @@ class MeshSuperentityStorageLayer< MeshConfig,
                                    MeshDimensionsTag< EntityTopology::dimensions >,
                                    true >
 {
-   using DimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
-   static constexpr int Dimensions = DimensionsTag::value;
+   using SuperdimensionsTag = MeshDimensionsTag< EntityTopology::dimensions >;
    using MeshTraitsType = MeshTraits< MeshConfig >;
-   using SuperentityTraitsType = typename MeshTraitsType::template SuperentityTraits< EntityTopology, Dimensions >;
-   using ThisType = MeshSuperentityStorageLayer< MeshConfig, EntityTopology, DimensionsTag, true >;
+   using SuperentityTraitsType = typename MeshTraitsType::template SuperentityTraits< EntityTopology, SuperdimensionsTag::value >;
+   using ThisType = MeshSuperentityStorageLayer< MeshConfig, EntityTopology, SuperdimensionsTag, true >;
 
 protected:
    using GlobalIndexType    = typename SuperentityTraitsType::GlobalIndexType;
@@ -173,6 +193,11 @@ protected:
    /****
     * These methods are due to 'using BaseType::...;' in the derived classes.
     */
+   bool setNumberOfEntities( const GlobalIndexType& entitiesCount )
+   {
+      return true;
+   }
+
    void print( std::ostream& str ) const {}
 
    bool operator==( const ThisType& layer ) const
@@ -190,7 +215,7 @@ protected:
       return true;
    }
  
-   void getSuperentityStorageNetwork( DimensionsTag ) {}
+   void getSuperentityStorageNetwork( SuperdimensionsTag ) {}
 };
 
 } // namespace Meshes
diff --git a/src/TNL/Meshes/MeshDetails/traits/MeshEntityTraits.h b/src/TNL/Meshes/MeshDetails/traits/MeshEntityTraits.h
index ed4352b732..06662323c3 100644
--- a/src/TNL/Meshes/MeshDetails/traits/MeshEntityTraits.h
+++ b/src/TNL/Meshes/MeshDetails/traits/MeshEntityTraits.h
@@ -66,7 +66,10 @@ class MeshEntityTraits
 {
 public:
    static constexpr bool storageEnabled = MeshConfig::entityStorage( Dimensions );
-   static constexpr bool orientationNeeded = MeshEntityOrientationNeeded< MeshConfig, MeshDimensionsTag< Dimensions > >::value;
+   // FIXME
+//   static constexpr bool orientationNeeded = MeshEntityOrientationNeeded< MeshConfig, MeshDimensionsTag< Dimensions > >::value;
+//   static constexpr bool orientationNeeded = false;
+   static constexpr bool orientationNeeded = ( Dimensions > 0 && Dimensions < MeshConfig::meshDimensions );
 
    using GlobalIndexType               = typename MeshConfig::GlobalIndexType;
    using LocalIndexType                = typename MeshConfig::LocalIndexType;
@@ -78,7 +81,6 @@ public:
    using Key                           = MeshEntitySeedKey< MeshConfig, EntityTopology >;
 
    using StorageArrayType              = Containers::Array< EntityType, Devices::Host, GlobalIndexType >;
-   using AccessArrayType               = Containers::SharedArray< EntityType, Devices::Host, GlobalIndexType >;
    using UniqueContainerType           = Containers::IndexedSet< EntityType, GlobalIndexType, Key >;
    using SeedIndexedSetType            = Containers::IndexedSet< SeedType, GlobalIndexType, Key >;
    using SeedArrayType                 = Containers::Array< SeedType, Devices::Host, GlobalIndexType >;
diff --git a/src/TNL/Meshes/MeshDetails/traits/MeshSubentityTraits.h b/src/TNL/Meshes/MeshDetails/traits/MeshSubentityTraits.h
index bad18de3e8..20935145d5 100644
--- a/src/TNL/Meshes/MeshDetails/traits/MeshSubentityTraits.h
+++ b/src/TNL/Meshes/MeshDetails/traits/MeshSubentityTraits.h
@@ -21,6 +21,7 @@
 #include <TNL/Meshes/MeshEntity.h>
 #include <TNL/Meshes/MeshConfigBase.h>
 #include <TNL/Meshes/Topologies/MeshEntityTopology.h>
+#include <TNL/Experimental/Multimaps/StaticEllpackIndexMultimap.h>
 
 namespace TNL {
 namespace Meshes {
@@ -47,6 +48,12 @@ public:
 
    static constexpr int count = MeshSubtopology< EntityTopology, Dimensions >::count;
 
+   /****
+    * Type of container for storing of the superentities indices.
+    */
+   using StorageNetworkType     = StaticEllpackIndexMultimap< count, GlobalIndexType, Devices::Host, LocalIndexType >;
+   using SubentityAccessorType  = typename StorageNetworkType::ValuesAccessorType;
+
    using StorageArrayType       = Containers::StaticArray< count, GlobalIndexType >;
    using IdArrayType            = Containers::StaticArray< count, GlobalIndexType >;
    using SeedArrayType          = Containers::StaticArray< count, Seed >;
diff --git a/src/TNL/Meshes/MeshEntity.h b/src/TNL/Meshes/MeshEntity.h
index b530b3bc51..d6a6d1160b 100644
--- a/src/TNL/Meshes/MeshEntity.h
+++ b/src/TNL/Meshes/MeshEntity.h
@@ -22,9 +22,9 @@
 #include <TNL/Meshes/MeshDetails/traits/MeshTraits.h>
 #include <TNL/Meshes/MeshDimensionTag.h>
 #include <TNL/Meshes/Topologies/MeshVertexTopology.h>
-#include <TNL/Meshes/MeshDetails/layers/MeshSubentityStorageLayer.h>
+#include <TNL/Meshes/MeshDetails/layers/MeshSubentityAccess.h>
 #include <TNL/Meshes/MeshDetails/layers/MeshSuperentityAccess.h>
-#include <TNL/Meshes/MeshDetails/layers/MeshSuperentityStorageRebinder.h>
+#include <TNL/Meshes/MeshDetails/layers/MeshEntityStorageRebinder.h>
 #include <TNL/Meshes/MeshDetails/initializer/MeshEntitySeed.h>
 
 namespace TNL {
@@ -36,7 +36,7 @@ class MeshInitializer;
 template< typename MeshConfig,
           typename EntityTopology_ >
 class MeshEntity
-   : protected MeshSubentityStorageLayers< MeshConfig, EntityTopology_ >,
+   : protected MeshSubentityAccess< MeshConfig, EntityTopology_ >,
      protected MeshSuperentityAccess< MeshConfig, EntityTopology_ >,
      public MeshEntityId< typename MeshConfig::IdType,
                           typename MeshConfig::GlobalIndexType >
@@ -72,9 +72,9 @@ class MeshEntity
       /****
        * Subentities
        */
-      using MeshSubentityStorageLayers< MeshConfig, EntityTopology_ >::getNumberOfSubentities;
-      using MeshSubentityStorageLayers< MeshConfig, EntityTopology_ >::getSubentityIndex;
-      using MeshSubentityStorageLayers< MeshConfig, EntityTopology_ >::getSubentityOrientation;
+      using MeshSubentityAccess< MeshConfig, EntityTopology_ >::getNumberOfSubentities;
+      using MeshSubentityAccess< MeshConfig, EntityTopology_ >::getSubentityIndex;
+      using MeshSubentityAccess< MeshConfig, EntityTopology_ >::getSubentityOrientation;
 
       /****
        * Superentities
@@ -93,8 +93,9 @@ class MeshEntity
       /****
        * Methods for the mesh initialization
        */
-      using MeshSubentityStorageLayers< MeshConfig, EntityTopology_ >::setSubentityIndex;
-      using MeshSubentityStorageLayers< MeshConfig, EntityTopology_ >::subentityOrientationsArray;
+      using MeshSubentityAccess< MeshConfig, EntityTopology_ >::bindSubentitiesStorageNetwork;
+      using MeshSubentityAccess< MeshConfig, EntityTopology_ >::setSubentityIndex;
+      using MeshSubentityAccess< MeshConfig, EntityTopology_ >::subentityOrientationsArray;
 
       using MeshSuperentityAccess< MeshConfig, EntityTopology_ >::bindSuperentitiesStorageNetwork;
       using MeshSuperentityAccess< MeshConfig, EntityTopology_ >::setNumberOfSuperentities;
@@ -103,7 +104,7 @@ class MeshEntity
    friend MeshInitializer< MeshConfig >;
 
    template< typename Mesh, typename DimensionsTag, typename SuperdimensionsTag >
-   friend struct MeshSuperentityStorageRebinderWorker;
+   friend struct MeshEntityStorageRebinderWorker;
 };
 
 /****
@@ -164,7 +165,7 @@ class MeshEntity< MeshConfig, MeshVertexTopology >
    friend MeshInitializer< MeshConfig >;
 
    template< typename Mesh, typename DimensionsTag, typename SuperdimensionsTag >
-   friend struct MeshSuperentityStorageRebinderWorker;
+   friend struct MeshEntityStorageRebinderWorker;
 };
 
 template< typename MeshConfig,
diff --git a/src/UnitTests/Meshes/MeshEntityTest.h b/src/UnitTests/Meshes/MeshEntityTest.h
index 5ef3bcb89f..3594cc1fd9 100644
--- a/src/UnitTests/Meshes/MeshEntityTest.h
+++ b/src/UnitTests/Meshes/MeshEntityTest.h
@@ -50,7 +50,10 @@ public:
 };
 
 template< typename MeshConfig, typename EntityTopology, int Dimensions >
-using StorageNetwork = typename MeshSuperentityTraits< MeshConfig, EntityTopology, Dimensions >::StorageNetworkType;
+using SubentityStorage = typename MeshSubentityTraits< MeshConfig, EntityTopology, Dimensions >::StorageNetworkType;
+
+template< typename MeshConfig, typename EntityTopology, int Dimensions >
+using SuperentityStorage = typename MeshSuperentityTraits< MeshConfig, EntityTopology, Dimensions >::StorageNetworkType;
 
 // stupid wrapper around MeshEntity to expose protected members needed for tests
 template< typename MeshConfig, typename EntityTopology >
@@ -60,6 +63,12 @@ class TestMeshEntity
    using BaseType = MeshEntity< MeshConfig, EntityTopology >;
 
 public:
+   template< int Subdimensions, typename Storage >
+   void bindSubentitiesStorageNetwork( const Storage& storage )
+   {
+      BaseType::template bindSubentitiesStorageNetwork< Subdimensions >( storage );
+   }
+
    template< int Subdimensions >
    void setSubentityIndex( const typename BaseType::LocalIndexType& localIndex,
                            const typename BaseType::GlobalIndexType& globalIndex )
@@ -134,10 +143,17 @@ TEST( MeshEntityTest, EdgeMeshEntityTest )
    ASSERT_TRUE( vertexEntities[ 2 ].getPoint() == point2 );
 
    Containers::StaticArray< 3, EdgeMeshEntityType > edgeEntities;
+   SubentityStorage< TestTriangleMeshConfig, MeshEdgeTopology, 0 > edgeVertexSubentities;
+   edgeVertexSubentities.setKeysRange( 3 );
+   ASSERT_TRUE( edgeVertexSubentities.allocate() );
+
+   edgeEntities[ 0 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 0 ) );
    edgeEntities[ 0 ].template setSubentityIndex< 0 >( 0, 0 );
    edgeEntities[ 0 ].template setSubentityIndex< 0 >( 1, 1 );
+   edgeEntities[ 1 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 1 ) );
    edgeEntities[ 1 ].template setSubentityIndex< 0 >( 0, 1 );
    edgeEntities[ 1 ].template setSubentityIndex< 0 >( 1, 2 );
+   edgeEntities[ 2 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 2 ) );
    edgeEntities[ 2 ].template setSubentityIndex< 0 >( 0, 2 );
    edgeEntities[ 2 ].template setSubentityIndex< 0 >( 1, 0 );
    edgeEntities[ 0 ].setId( 0 );
@@ -182,10 +198,17 @@ TEST( MeshEntityTest, TriangleMeshEntityTest )
    ASSERT_TRUE( vertexEntities[ 2 ].getPoint() == point2 );
 
    Containers::StaticArray< 3, EdgeMeshEntityType > edgeEntities;
+   SubentityStorage< TestTriangleMeshConfig, MeshEdgeTopology, 0 > edgeVertexSubentities;
+   edgeVertexSubentities.setKeysRange( 3 );
+   ASSERT_TRUE( edgeVertexSubentities.allocate() );
+
+   edgeEntities[ 0 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 0 ) );
    edgeEntities[ 0 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTriangleTopology, MeshEdgeTopology, 0, 0 >::index );
    edgeEntities[ 0 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTriangleTopology, MeshEdgeTopology, 0, 1 >::index );
+   edgeEntities[ 1 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 1 ) );
    edgeEntities[ 1 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTriangleTopology, MeshEdgeTopology, 1, 0 >::index );
    edgeEntities[ 1 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTriangleTopology, MeshEdgeTopology, 1, 1 >::index );
+   edgeEntities[ 2 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 2 ) );
    edgeEntities[ 2 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTriangleTopology, MeshEdgeTopology, 2, 0 >::index );
    edgeEntities[ 2 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTriangleTopology, MeshEdgeTopology, 2, 1 >::index );
 
@@ -197,7 +220,14 @@ TEST( MeshEntityTest, TriangleMeshEntityTest )
    ASSERT_TRUE( edgeEntities[ 2 ].getVertexIndex( 1 ) == ( tnlSubentityVertex< MeshTriangleTopology, MeshEdgeTopology, 2, 1 >::index ) );
 
    TriangleMeshEntityType triangleEntity;
-
+   SubentityStorage< TestTriangleMeshConfig, MeshTriangleTopology, 0 > triangleVertexSubentities;
+   SubentityStorage< TestTriangleMeshConfig, MeshTriangleTopology, 1 > triangleEdgeSubentities;
+   triangleVertexSubentities.setKeysRange( 1 );
+   triangleEdgeSubentities.setKeysRange( 1 );
+   ASSERT_TRUE( triangleVertexSubentities.allocate() );
+   ASSERT_TRUE( triangleEdgeSubentities.allocate() );
+
+   triangleEntity.template bindSubentitiesStorageNetwork< 0 >( triangleVertexSubentities.getValues( 0 ) );
    triangleEntity.template setSubentityIndex< 0 >( 0 , 0 );
    triangleEntity.template setSubentityIndex< 0 >( 1 , 1 );
    triangleEntity.template setSubentityIndex< 0 >( 2 , 2 );
@@ -206,6 +236,7 @@ TEST( MeshEntityTest, TriangleMeshEntityTest )
    ASSERT_TRUE( triangleEntity.template getSubentityIndex< 0 >( 1 ) == 1 );
    ASSERT_TRUE( triangleEntity.template getSubentityIndex< 0 >( 2 ) == 2 );
 
+   triangleEntity.template bindSubentitiesStorageNetwork< 1 >( triangleEdgeSubentities.getValues( 0 ) );
    triangleEntity.template setSubentityIndex< 1 >( 0 , 0 );
    triangleEntity.template setSubentityIndex< 1 >( 1 , 1 );
    triangleEntity.template setSubentityIndex< 1 >( 2 , 2 );
@@ -255,17 +286,27 @@ TEST( MeshEntityTest, TetragedronMeshEntityTest )
    ASSERT_TRUE( vertexEntities[ 3 ].getPoint() == point3 );
 
    Containers::StaticArray< MeshSubtopology< MeshTetrahedronTopology, 1 >::count,
-                   EdgeMeshEntityType > edgeEntities;
+                            EdgeMeshEntityType > edgeEntities;
+   SubentityStorage< TestTriangleMeshConfig, MeshEdgeTopology, 0 > edgeVertexSubentities;
+   edgeVertexSubentities.setKeysRange( 6 );
+   ASSERT_TRUE( edgeVertexSubentities.allocate() );
+
+   edgeEntities[ 0 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 0 ) );
    edgeEntities[ 0 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 0, 0 >::index );
    edgeEntities[ 0 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 0, 1 >::index );
+   edgeEntities[ 1 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 1 ) );
    edgeEntities[ 1 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 1, 0 >::index );
    edgeEntities[ 1 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 1, 1 >::index );
+   edgeEntities[ 2 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 2 ) );
    edgeEntities[ 2 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 2, 0 >::index );
    edgeEntities[ 2 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 2, 1 >::index );
+   edgeEntities[ 3 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 3 ) );
    edgeEntities[ 3 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 3, 0 >::index );
    edgeEntities[ 3 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 3, 1 >::index );
+   edgeEntities[ 4 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 4 ) );
    edgeEntities[ 4 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 4, 0 >::index );
    edgeEntities[ 4 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 4, 1 >::index );
+   edgeEntities[ 5 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 5 ) );
    edgeEntities[ 5 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 5, 0 >::index );
    edgeEntities[ 5 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 5, 1 >::index );
 
@@ -283,16 +324,24 @@ TEST( MeshEntityTest, TetragedronMeshEntityTest )
    ASSERT_TRUE( edgeEntities[ 5 ].getVertexIndex( 1 ) == ( tnlSubentityVertex< MeshTetrahedronTopology, MeshEdgeTopology, 5, 1 >::index ) );
 
    Containers::StaticArray< MeshSubtopology< MeshTetrahedronTopology, 2 >::count,
-                   TriangleMeshEntityType > triangleEntities;
+                            TriangleMeshEntityType > triangleEntities;
+   SubentityStorage< TestTriangleMeshConfig, MeshTriangleTopology, 0 > triangleVertexSubentities;
+   triangleVertexSubentities.setKeysRange( 4 );
+   ASSERT_TRUE( triangleVertexSubentities.allocate() );
+
+   triangleEntities[ 0 ].template bindSubentitiesStorageNetwork< 0 >( triangleVertexSubentities.getValues( 0 ) );
    triangleEntities[ 0 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 0, 0 >::index );
    triangleEntities[ 0 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 0, 1 >::index );
    triangleEntities[ 0 ].template setSubentityIndex< 0 >( 2, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 0, 2 >::index );
+   triangleEntities[ 1 ].template bindSubentitiesStorageNetwork< 0 >( triangleVertexSubentities.getValues( 1 ) );
    triangleEntities[ 1 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 1, 0 >::index );
    triangleEntities[ 1 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 1, 1 >::index );
    triangleEntities[ 1 ].template setSubentityIndex< 0 >( 2, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 1, 2 >::index );
+   triangleEntities[ 2 ].template bindSubentitiesStorageNetwork< 0 >( triangleVertexSubentities.getValues( 2 ) );
    triangleEntities[ 2 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 2, 0 >::index );
    triangleEntities[ 2 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 2, 1 >::index );
    triangleEntities[ 2 ].template setSubentityIndex< 0 >( 2, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 2, 2 >::index );
+   triangleEntities[ 3 ].template bindSubentitiesStorageNetwork< 0 >( triangleVertexSubentities.getValues( 3 ) );
    triangleEntities[ 3 ].template setSubentityIndex< 0 >( 0, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 3, 0 >::index );
    triangleEntities[ 3 ].template setSubentityIndex< 0 >( 1, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 3, 1 >::index );
    triangleEntities[ 3 ].template setSubentityIndex< 0 >( 2, tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 3, 2 >::index );
@@ -308,6 +357,17 @@ TEST( MeshEntityTest, TetragedronMeshEntityTest )
    ASSERT_TRUE( triangleEntities[ 2 ].getVertexIndex( 2 ) == ( tnlSubentityVertex< MeshTetrahedronTopology, MeshTriangleTopology, 2, 2 >::index ) );
 
    TetrahedronMeshEntityType tetrahedronEntity;
+   SubentityStorage< TestTriangleMeshConfig, MeshTetrahedronTopology, 0 > tetrahedronVertexSubentities;
+   SubentityStorage< TestTriangleMeshConfig, MeshTetrahedronTopology, 1 > tetrahedronEdgeSubentities;
+   SubentityStorage< TestTriangleMeshConfig, MeshTetrahedronTopology, 2 > tetrahedronTriangleSubentities;
+   tetrahedronVertexSubentities.setKeysRange( 1 );
+   tetrahedronEdgeSubentities.setKeysRange( 1 );
+   tetrahedronTriangleSubentities.setKeysRange( 1 );
+   ASSERT_TRUE( tetrahedronVertexSubentities.allocate() );
+   ASSERT_TRUE( tetrahedronEdgeSubentities.allocate() );
+   ASSERT_TRUE( tetrahedronTriangleSubentities.allocate() );
+
+   tetrahedronEntity.template bindSubentitiesStorageNetwork< 0 >( tetrahedronVertexSubentities.getValues( 0 ) );
    tetrahedronEntity.template setSubentityIndex< 0 >( 0, 0 );
    tetrahedronEntity.template setSubentityIndex< 0 >( 1, 1 );
    tetrahedronEntity.template setSubentityIndex< 0 >( 2, 2 );
@@ -318,6 +378,7 @@ TEST( MeshEntityTest, TetragedronMeshEntityTest )
    ASSERT_TRUE( tetrahedronEntity.getVertexIndex( 2 ) == 2 );
    ASSERT_TRUE( tetrahedronEntity.getVertexIndex( 3 ) == 3 );
 
+   tetrahedronEntity.template bindSubentitiesStorageNetwork< 2 >( tetrahedronTriangleSubentities.getValues( 0 ) );
    tetrahedronEntity.template setSubentityIndex< 2 >( 0, 0 );
    tetrahedronEntity.template setSubentityIndex< 2 >( 1, 1 );
    tetrahedronEntity.template setSubentityIndex< 2 >( 2, 2 );
@@ -328,6 +389,7 @@ TEST( MeshEntityTest, TetragedronMeshEntityTest )
    ASSERT_TRUE( tetrahedronEntity.template getSubentityIndex< 2 >( 2 ) == 2 );
    ASSERT_TRUE( tetrahedronEntity.template getSubentityIndex< 2 >( 3 ) == 3 );
 
+   tetrahedronEntity.template bindSubentitiesStorageNetwork< 1 >( tetrahedronEdgeSubentities.getValues( 0 ) );
    tetrahedronEntity.template setSubentityIndex< 1 >( 0, 0 );
    tetrahedronEntity.template setSubentityIndex< 1 >( 1, 1 );
    tetrahedronEntity.template setSubentityIndex< 1 >( 2, 2 );
@@ -395,14 +457,23 @@ TEST( MeshEntityTest, TwoTrianglesMeshEntityTest )
    ASSERT_TRUE( vertexEntities[ 3 ].getPoint() == point3 );
 
    Containers::StaticArray< 5, EdgeMeshEntityType > edgeEntities;
+   SubentityStorage< TestTriangleMeshConfig, MeshEdgeTopology, 0 > edgeVertexSubentities;
+   edgeVertexSubentities.setKeysRange( 5 );
+   ASSERT_TRUE( edgeVertexSubentities.allocate() );
+
+   edgeEntities[ 0 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 0 ) );
    edgeEntities[ 0 ].template setSubentityIndex< 0 >( 0, 1 );
    edgeEntities[ 0 ].template setSubentityIndex< 0 >( 1, 2 );
+   edgeEntities[ 1 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 1 ) );
    edgeEntities[ 1 ].template setSubentityIndex< 0 >( 0, 2 );
    edgeEntities[ 1 ].template setSubentityIndex< 0 >( 1, 0 );
+   edgeEntities[ 2 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 2 ) );
    edgeEntities[ 2 ].template setSubentityIndex< 0 >( 0, 0 );
    edgeEntities[ 2 ].template setSubentityIndex< 0 >( 1, 1 );
+   edgeEntities[ 3 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 3 ) );
    edgeEntities[ 3 ].template setSubentityIndex< 0 >( 0, 2 );
    edgeEntities[ 3 ].template setSubentityIndex< 0 >( 1, 3 );
+   edgeEntities[ 4 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 4 ) );
    edgeEntities[ 4 ].template setSubentityIndex< 0 >( 0, 3 );
    edgeEntities[ 4 ].template setSubentityIndex< 0 >( 1, 1 );
 
@@ -418,16 +489,26 @@ TEST( MeshEntityTest, TwoTrianglesMeshEntityTest )
    ASSERT_TRUE( edgeEntities[ 4 ].getVertexIndex( 1 ) == 1 );
 
    Containers::StaticArray< 2, TriangleMeshEntityType > triangleEntities;
-
+   SubentityStorage< TestTriangleMeshConfig, MeshTriangleTopology, 0 > triangleVertexSubentities;
+   SubentityStorage< TestTriangleMeshConfig, MeshTriangleTopology, 1 > triangleEdgeSubentities;
+   triangleVertexSubentities.setKeysRange( 2 );
+   ASSERT_TRUE( triangleVertexSubentities.allocate() );
+   triangleEdgeSubentities.setKeysRange( 2 );
+   ASSERT_TRUE( triangleEdgeSubentities.allocate() );
+
+   triangleEntities[ 0 ].template bindSubentitiesStorageNetwork< 0 >( triangleVertexSubentities.getValues( 0 ) );
    triangleEntities[ 0 ].template setSubentityIndex< 0 >( 0 , 0 );
    triangleEntities[ 0 ].template setSubentityIndex< 0 >( 1 , 1 );
    triangleEntities[ 0 ].template setSubentityIndex< 0 >( 2 , 2 );
+   triangleEntities[ 0 ].template bindSubentitiesStorageNetwork< 1 >( triangleEdgeSubentities.getValues( 0 ) );
    triangleEntities[ 0 ].template setSubentityIndex< 1 >( 0 , 0 );
    triangleEntities[ 0 ].template setSubentityIndex< 1 >( 1 , 1 );
    triangleEntities[ 0 ].template setSubentityIndex< 1 >( 2 , 2 );
+   triangleEntities[ 1 ].template bindSubentitiesStorageNetwork< 0 >( triangleVertexSubentities.getValues( 1 ) );
    triangleEntities[ 1 ].template setSubentityIndex< 0 >( 0 , 1 );
    triangleEntities[ 1 ].template setSubentityIndex< 0 >( 1 , 2 );
    triangleEntities[ 1 ].template setSubentityIndex< 0 >( 2 , 3 );
+   triangleEntities[ 1 ].template bindSubentitiesStorageNetwork< 1 >( triangleEdgeSubentities.getValues( 1 ) );
    triangleEntities[ 1 ].template setSubentityIndex< 1 >( 0 , 3 );
    triangleEntities[ 1 ].template setSubentityIndex< 1 >( 1 , 4 );
    triangleEntities[ 1 ].template setSubentityIndex< 1 >( 2 , 0 );
@@ -449,7 +530,7 @@ TEST( MeshEntityTest, TwoTrianglesMeshEntityTest )
    /*
     * Tests for the superentities layer.
     */
-   StorageNetwork< TestTriangleMeshConfig, MeshVertexTopology, 1 > vertexEdgeSuperentities;
+   SuperentityStorage< TestTriangleMeshConfig, MeshVertexTopology, 1 > vertexEdgeSuperentities;
    vertexEdgeSuperentities.setKeysRange( 4 );
    vertexEdgeSuperentities.allocate( 3 );
 
@@ -474,7 +555,7 @@ TEST( MeshEntityTest, TwoTrianglesMeshEntityTest )
    ASSERT_EQ( vertexEntities[ 1 ].template getSuperentityIndex< 1 >( 2 ),    4 );
 
 
-   StorageNetwork< TestTriangleMeshConfig, MeshVertexTopology, 2 > vertexCellSuperentities;
+   SuperentityStorage< TestTriangleMeshConfig, MeshVertexTopology, 2 > vertexCellSuperentities;
    vertexCellSuperentities.setKeysRange( 4 );
    vertexCellSuperentities.allocate( 2 );
 
@@ -488,7 +569,7 @@ TEST( MeshEntityTest, TwoTrianglesMeshEntityTest )
    ASSERT_EQ( vertexEntities[ 1 ].template getSuperentityIndex< 2 >( 1 ),    1 );
 
 
-   StorageNetwork< TestTriangleMeshConfig, MeshEdgeTopology, 2 > edgeCellSuperentities;
+   SuperentityStorage< TestTriangleMeshConfig, MeshEdgeTopology, 2 > edgeCellSuperentities;
    edgeCellSuperentities.setKeysRange( 5 );
    edgeCellSuperentities.allocate( 2 );
 
@@ -528,23 +609,39 @@ TEST( MeshEntityTest, OneTriangleComparisonTest )
    vertices[ 2 ].setPoint( point2 );
 
    Containers::StaticArray< 3, EdgeMeshEntityType > edges;
+   SubentityStorage< TestTriangleMeshConfig, MeshEdgeTopology, 0 > edgeVertexSubentities;
+   edgeVertexSubentities.setKeysRange( 3 );
+   ASSERT_TRUE( edgeVertexSubentities.allocate() );
+
+   edges[ 0 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 0 ) );
    edges[ 0 ].template setSubentityIndex< 0 >( 0, 1 );
    edges[ 0 ].template setSubentityIndex< 0 >( 1, 2 );
+   edges[ 1 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 1 ) );
    edges[ 1 ].template setSubentityIndex< 0 >( 0, 2 );
    edges[ 1 ].template setSubentityIndex< 0 >( 1, 0 );
+   edges[ 2 ].template bindSubentitiesStorageNetwork< 0 >( edgeVertexSubentities.getValues( 2 ) );
    edges[ 2 ].template setSubentityIndex< 0 >( 0, 0 );
    edges[ 2 ].template setSubentityIndex< 0 >( 1, 1 );
 
    TriangleMeshEntityType triangle;
+   SubentityStorage< TestTriangleMeshConfig, MeshTriangleTopology, 0 > triangleVertexSubentities;
+   SubentityStorage< TestTriangleMeshConfig, MeshTriangleTopology, 1 > triangleEdgeSubentities;
+   triangleVertexSubentities.setKeysRange( 1 );
+   triangleEdgeSubentities.setKeysRange( 1 );
+   ASSERT_TRUE( triangleVertexSubentities.allocate() );
+   ASSERT_TRUE( triangleEdgeSubentities.allocate() );
+
+   triangle.template bindSubentitiesStorageNetwork< 0 >( triangleVertexSubentities.getValues( 0 ) );
    triangle.template setSubentityIndex< 0 >( 0 , 0 );
    triangle.template setSubentityIndex< 0 >( 1 , 1 );
    triangle.template setSubentityIndex< 0 >( 2 , 2 );
+   triangle.template bindSubentitiesStorageNetwork< 1 >( triangleVertexSubentities.getValues( 0 ) );
    triangle.template setSubentityIndex< 1 >( 0 , 0 );
    triangle.template setSubentityIndex< 1 >( 1 , 1 );
    triangle.template setSubentityIndex< 1 >( 2 , 2 );
 
 
-   StorageNetwork< TestTriangleMeshConfig, MeshVertexTopology, 1 > vertexEdgeSuperentities;
+   SuperentityStorage< TestTriangleMeshConfig, MeshVertexTopology, 1 > vertexEdgeSuperentities;
    vertexEdgeSuperentities.setKeysRange( 3 );
    vertexEdgeSuperentities.allocate( 2 );
 
@@ -564,7 +661,7 @@ TEST( MeshEntityTest, OneTriangleComparisonTest )
    vertices[ 2 ].template setSuperentityIndex< 1 >( 1, 1 );
 
 
-   StorageNetwork< TestTriangleMeshConfig, MeshVertexTopology, 2 > vertexCellSuperentities;
+   SuperentityStorage< TestTriangleMeshConfig, MeshVertexTopology, 2 > vertexCellSuperentities;
    vertexCellSuperentities.setKeysRange( 3 );
    vertexCellSuperentities.allocate( 1 );
 
@@ -581,7 +678,7 @@ TEST( MeshEntityTest, OneTriangleComparisonTest )
    vertices[ 2 ].template setSuperentityIndex< 2 >( 0, 0 );
 
 
-   StorageNetwork< TestTriangleMeshConfig, MeshEdgeTopology, 2 > edgeCellSuperentities;
+   SuperentityStorage< TestTriangleMeshConfig, MeshEdgeTopology, 2 > edgeCellSuperentities;
    edgeCellSuperentities.setKeysRange( 3 );
    edgeCellSuperentities.allocate( 1 );
 
-- 
GitLab