Commit 99cf6c69 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Replaced multimaps in Mesh with segments-based binary matrices

parent fae134f6

src/TNL/Matrices/MatrixPermutationApplier.h

0 → 100644
+97 −0
Original line number Diff line number Diff line 
/***************************************************************************

                          MatrixPermutationApplier.h  -  description

                             -------------------

    begin                : May 7, 2020

    copyright            : (C) 2020 by Tomas Oberhuber et al.

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/* See Copyright Notice in tnl/Copyright */



// Implemented by: Jakub Klinkovský



#pragma once



#include <TNL/Algorithms/ParallelFor.h>



namespace TNL {

namespace Matrices {



template< typename Matrix,

          typename PermutationArray >

void permuteMatrixRows( Matrix& matrix, const PermutationArray& perm )

{

   static_assert( std::is_same< typename Matrix::DeviceType, typename PermutationArray::DeviceType >::value,

                  "The matrix and permutation vector must be stored on the same device." );

   using IndexType = typename Matrix::IndexType;

   using DeviceType = typename Matrix::DeviceType;

   TNL_ASSERT_EQ( matrix.getRows(), perm.getSize(), "permutation size does not match the matrix size" );



   // FIXME: getConstView does not work

//   const auto matrix_view = matrix.getConstView();

   const auto matrix_view = matrix.getView();

   const auto perm_view = perm.getConstView();



   // create temporary matrix for the permuted data

   Matrix matrixCopy;

   matrixCopy.setLike( matrix );



   // permute the row capacities

   typename Matrix::RowsCapacitiesType capacities( matrix.getRows() );

   auto capacities_view = capacities.getView();



   auto kernel_capacities = [=] __cuda_callable__ ( IndexType i ) mutable

   {

      capacities_view[ i ] = matrix_view.getRowCapacity( perm_view[ i ] );

   };

   Algorithms::ParallelFor< DeviceType >::exec( (IndexType) 0, matrix.getRows(), kernel_capacities );



   matrixCopy.setRowCapacities( capacities );

   auto copy_view = matrixCopy.getView();



   auto kernel = [=] __cuda_callable__ ( IndexType i ) mutable

   {

      const auto srcRow = matrix_view.getRow( perm_view[ i ] );

      auto destRow = copy_view.getRow( i );

      for( IndexType c = 0; c < srcRow.getSize(); c++ )

         if( srcRow.isBinary() )

            destRow.setElement( c, srcRow.getColumnIndex( c ), true );  // the value does not matter

         else

            destRow.setElement( c, srcRow.getColumnIndex( c ), srcRow.getValue( c ) );

   };

   Algorithms::ParallelFor< DeviceType >::exec( (IndexType) 0, matrix.getRows(), kernel );



   // copy the permuted data back into the matrix

   matrix = matrixCopy;

}



template< typename Matrix,

          typename PermutationArray >

void permuteMatrixColumns( Matrix& matrix, const PermutationArray& iperm )

{

   static_assert( std::is_same< typename Matrix::DeviceType, typename PermutationArray::DeviceType >::value,

                  "The matrix and permutation vector must be stored on the same device." );

   using IndexType = typename Matrix::IndexType;

   using DeviceType = typename Matrix::DeviceType;



   auto matrix_view = matrix.getView();

   const auto iperm_view = iperm.getConstView();



   auto kernel = [=] __cuda_callable__ ( IndexType i ) mutable

   {

      auto row = matrix_view.getRow( i );

      for( IndexType c = 0; c < row.getSize(); c++ ) {

         const IndexType col = row.getColumnIndex( c );

         if( col == matrix_view.getPaddingIndex() )

            break;

         if( row.isBinary() )

            row.setElement( c, iperm_view[ col ], true );  // the value does not matter

         else

            row.setElement( c, iperm_view[ col ], row.getValue( c ) );

      }

   };

   Algorithms::ParallelFor< DeviceType >::exec( (IndexType) 0, matrix.getRows(), kernel );

}



} // namespace Matrices

} // namespace TNL

src/TNL/Meshes/Mesh.h

+2 −2
Original line number Diff line number Diff line
@@ -210,8 +210,8 @@ class Mesh 
      // Methods for the mesh initializer

      using StorageBaseType::getPoints;

      using StorageBaseType::setEntitiesCount;

      using StorageBaseType::getSubentityStorageNetwork;

      using StorageBaseType::getSuperentityStorageNetwork;

      using StorageBaseType::getSubentitiesMatrix;

      using StorageBaseType::getSuperentitiesMatrix;



      friend Initializer< MeshConfig >;

src/TNL/Meshes/Mesh.hpp

+23 −9
Original line number Diff line number Diff line
@@ -175,7 +175,7 @@ constexpr typename Mesh< MeshConfig, Device >::LocalIndexType 
Mesh< MeshConfig, Device >::

getSubentitiesCount( const GlobalIndexType entityIndex ) const

{

   return this->template getSubentityStorageNetwork< EntityDimension, SubentityDimension >().getValuesCount( entityIndex );

   return StorageBaseType::template getSubentitiesCount< EntityDimension, SubentityDimension >();

}



template< typename MeshConfig, typename Device >
@@ -185,7 +185,9 @@ typename Mesh< MeshConfig, Device >::GlobalIndexType 
Mesh< MeshConfig, Device >::

getSubentityIndex( const GlobalIndexType entityIndex, const LocalIndexType subentityIndex ) const

{

   return this->template getSubentityStorageNetwork< EntityDimension, SubentityDimension >().getValue( entityIndex, subentityIndex );

   const auto& row = this->template getSubentitiesMatrix< EntityDimension, SubentityDimension >().getRow( entityIndex );

   TNL_ASSERT_GE( row.getColumnIndex( subentityIndex ), 0, "padding index returned for given subentity index" );

   return row.getColumnIndex( subentityIndex );

}



template< typename MeshConfig, typename Device >
@@ -195,7 +197,7 @@ typename Mesh< MeshConfig, Device >::LocalIndexType 
Mesh< MeshConfig, Device >::

getSuperentitiesCount( const GlobalIndexType entityIndex ) const

{

   return this->template getSuperentityStorageNetwork< EntityDimension, SuperentityDimension >().getValuesCount( entityIndex );

   return this->template getSuperentitiesCountsArray< EntityDimension, SuperentityDimension >()[ entityIndex ];

}



template< typename MeshConfig, typename Device >
@@ -205,7 +207,9 @@ typename Mesh< MeshConfig, Device >::GlobalIndexType 
Mesh< MeshConfig, Device >::

getSuperentityIndex( const GlobalIndexType entityIndex, const LocalIndexType superentityIndex ) const

{

   return this->template getSuperentityStorageNetwork< EntityDimension, SuperentityDimension >().getValue( entityIndex, superentityIndex );

   const auto row = this->template getSuperentitiesMatrix< EntityDimension, SuperentityDimension >().getRow( entityIndex );

   TNL_ASSERT_GE( row.getColumnIndex( superentityIndex ), 0, "padding index returned for given superentity index" );

   return row.getColumnIndex( superentityIndex );

}



template< typename MeshConfig, typename Device >
@@ -229,7 +233,9 @@ getCellNeighborIndex( const GlobalIndexType cellIndex, const LocalIndexType neig 
                  "You try to access the dual graph which is disabled in the mesh configuration." );

   TNL_ASSERT_GE( neighborIndex, 0, "Invalid cell neighbor index." );

   TNL_ASSERT_LT( neighborIndex, getCellNeighborsCount( cellIndex ), "Invalid cell neighbor index." );

   return this->getDualGraph().getValues( cellIndex )[ neighborIndex ];

   const auto row = this->getDualGraph().getRow( cellIndex );

   TNL_ASSERT_GE( row.getColumnIndex( neighborIndex ), 0, "padding index returned for given neighbor index" );

   return row.getColumnIndex( neighborIndex );

}
@@ -275,17 +281,25 @@ void 
Mesh< MeshConfig, Device >::

save( File& file ) const

{

   // saving via host is necessary due to segment-based sparse matrices

   if( std::is_same< Device, Devices::Cuda >::value ) {

      Mesh< MeshConfig, Devices::Host > hostMesh;

      hostMesh = *this;

      hostMesh.save( file );

   }

   else {

      Object::save( file );

      StorageBaseType::save( file );

      EntityTagsLayerFamily::save( file );

   }

}



template< typename MeshConfig, typename Device >

void

Mesh< MeshConfig, Device >::

load( File& file )

{

   // loading via host is necessary for the initialization of the dual graph

   // loading via host is necessary for the initialization of the dual graph (and due to segment-based sparse matrices)

   if( std::is_same< Device, Devices::Cuda >::value ) {

      Mesh< MeshConfig, Devices::Host > hostMesh;

      hostMesh.load( file );

src/TNL/Meshes/MeshDetails/IndexPermutationApplier.h

+12 −11
Original line number Diff line number Diff line
@@ -12,7 +12,7 @@ 


#include <TNL/Meshes/DimensionTag.h>

#include <TNL/Meshes/Mesh.h>

#include <TNL/Containers/Multimaps/MultimapPermutationApplier.h>

#include <TNL/Matrices/MatrixPermutationApplier.h>



namespace TNL {

namespace Meshes {
@@ -32,8 +32,8 @@ private: 
   {

      static void exec( Mesh& mesh, const GlobalIndexArray& perm )

      {

         auto& subentitiesStorage = mesh.template getSubentityStorageNetwork< Dimension, Subdimension >();

         Containers::Multimaps::permuteMultimapKeys( subentitiesStorage, perm );

         auto& subentitiesStorage = mesh.template getSubentitiesMatrix< Dimension, Subdimension >();

         Matrices::permuteMatrixRows( subentitiesStorage, perm );

      }

   };
@@ -53,8 +53,9 @@ private: 
   {

      static void exec( Mesh& mesh, const GlobalIndexArray& perm )

      {

         auto& superentitiesStorage = mesh.template getSuperentityStorageNetwork< Dimension, Superdimension >();

         Containers::Multimaps::permuteMultimapKeys( superentitiesStorage, perm );

         permuteArray( mesh.template getSuperentitiesCountsArray< Dimension, Superdimension >(), perm );

         auto& superentitiesStorage = mesh.template getSuperentitiesMatrix< Dimension, Superdimension >();

         Matrices::permuteMatrixRows( superentitiesStorage, perm );

      }

   };
@@ -74,8 +75,8 @@ private: 
   {

      static void exec( Mesh& mesh, const GlobalIndexArray& iperm )

      {

         auto& superentitiesStorage = mesh.template getSuperentityStorageNetwork< Subdimension, Dimension >();

         Containers::Multimaps::permuteMultimapValues( superentitiesStorage, iperm );

         auto& superentitiesStorage = mesh.template getSuperentitiesMatrix< Subdimension, Dimension >();

         Matrices::permuteMatrixColumns( superentitiesStorage, iperm );

      }

   };
@@ -95,8 +96,8 @@ private: 
   {

      static void exec( Mesh& mesh, const GlobalIndexArray& iperm )

      {

         auto& subentitiesStorage = mesh.template getSubentityStorageNetwork< Superdimension, Dimension >();

         Containers::Multimaps::permuteMultimapValues( subentitiesStorage, iperm );

         auto& subentitiesStorage = mesh.template getSubentitiesMatrix< Superdimension, Dimension >();

         Matrices::permuteMatrixColumns( subentitiesStorage, iperm );

      }

   };
@@ -125,8 +126,8 @@ private: 
   {

      permuteArray( mesh.getNeighborCounts(), perm );

      auto& graph = mesh.getDualGraph();

      Containers::Multimaps::permuteMultimapKeys( graph, perm );

      Containers::Multimaps::permuteMultimapValues( graph, iperm );

      Matrices::permuteMatrixRows( graph, perm );

      Matrices::permuteMatrixColumns( graph, iperm );

   }



   template< typename Mesh_, std::enable_if_t< ! Mesh_::Config::dualGraphStorage(), bool > = true >

src/TNL/Meshes/MeshDetails/initializer/EntityInitializer.h

+63 −41
Original line number Diff line number Diff line
@@ -65,7 +65,12 @@ class EntityInitializer 
   using InitializerType  = Initializer< MeshConfig >;



public:

   static void initEntity( const GlobalIndexType& entityIndex, const SeedType& entitySeed, InitializerType& initializer)

   static void initSubvertexMatrix( const GlobalIndexType entitiesCount, InitializerType& initializer )

   {

      initializer.template initSubentityMatrix< EntityTopology::dimension, 0 >( entitiesCount );

   }



   static void initEntity( const GlobalIndexType entityIndex, const SeedType& entitySeed, InitializerType& initializer )

   {

      // this is necessary if we want to use existing entities instead of intermediate seeds to create subentity seeds

      for( LocalIndexType i = 0; i < entitySeed.getCornerIds().getSize(); i++ )
@@ -86,7 +91,8 @@ class EntityInitializer< MeshConfig, EntityTopology, false > 
   using SeedType         = EntitySeed< MeshConfig, EntityTopology >;

   using InitializerType  = Initializer< MeshConfig >;

public:

   static void initEntity( const GlobalIndexType& entityIndex, const SeedType& entitySeed, InitializerType& initializer) {}

   static void initSubvertexMatrix( const GlobalIndexType entitiesCount, InitializerType& initializer ) {}

   static void initEntity( const GlobalIndexType entityIndex, const SeedType& entitySeed, InitializerType& initializer ) {}

};
@@ -123,22 +129,24 @@ class EntityInitializerLayer< MeshConfig, 
   using SuperentityTraitsType      = typename MeshTraits< MeshConfig >::template EntityTraits< SuperdimensionTag::value >;

   using SuperentityTopology        = typename SuperentityTraitsType::EntityTopology;

   using SubentitySeedsCreatorType  = SubentitySeedsCreator< MeshConfig, SuperdimensionTag, SubdimensionTag >;

   using SuperentityStorageNetwork  = typename MeshTraits< MeshConfig >::template SuperentityTraits< SubentityTopology, SuperdimensionTag::value >::StorageNetworkType;

   using ValuesAllocationVectorType = typename SuperentityStorageNetwork::ValuesAllocationVectorType;

   using SuperentityMatrixType      = typename MeshTraits< MeshConfig >::SuperentityMatrixType;



public:

   static void initSuperentities( InitializerType& meshInitializer, MeshType& mesh )

   {

      //std::cout << "   Initiating superentities with dimension " << SuperdimensionTag::value << " for subentities with dimension " << SubdimensionTag::value << " ... " << std::endl;



      // counter for superentities of each subentity

      const GlobalIndexType subentitiesCount = mesh.template getEntitiesCount< SubdimensionTag::value >();

      ValuesAllocationVectorType superentitiesCounts( subentitiesCount );

      const GlobalIndexType superentitiesCount = mesh.template getEntitiesCount< SuperdimensionTag::value >();

      if( SubdimensionTag::value > 0 )

         meshInitializer.template initSubentityMatrix< SuperdimensionTag::value, SubdimensionTag::value >( superentitiesCount, subentitiesCount );



      // counter for superentities of each subentity

      auto& superentitiesCounts = meshInitializer.template getSuperentitiesCountsArray< SubdimensionTag::value, SuperdimensionTag::value >();

      superentitiesCounts.setSize( subentitiesCount );

      superentitiesCounts.setValue( 0 );



      for( GlobalIndexType superentityIndex = 0;

           superentityIndex < mesh.template getEntitiesCount< SuperdimensionTag::value >();

           superentityIndex++ )

      for( GlobalIndexType superentityIndex = 0; superentityIndex < superentitiesCount; superentityIndex++ )

      {

         auto subentitySeeds = SubentitySeedsCreatorType::create( meshInitializer.template getSubvertices< SuperdimensionTag::value >( superentityIndex ) );

         for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )
@@ -150,21 +158,21 @@ public: 
      }



      // allocate superentities storage

      SuperentityStorageNetwork& superentityStorageNetwork = meshInitializer.template meshSuperentityStorageNetwork< SubdimensionTag::value, SuperdimensionTag::value >();

      superentityStorageNetwork.allocate( superentitiesCounts );

      SuperentityMatrixType& matrix = meshInitializer.template getSuperentitiesMatrix< SubdimensionTag::value, SuperdimensionTag::value >();

      matrix.setDimensions( subentitiesCount, superentitiesCount );

      matrix.setRowCapacities( superentitiesCounts );

      superentitiesCounts.setValue( 0 );



      // initialize superentities storage

      for( GlobalIndexType superentityIndex = 0;

           superentityIndex < mesh.template getEntitiesCount< SuperdimensionTag::value >();

           superentityIndex++ )

      for( GlobalIndexType superentityIndex = 0; superentityIndex < superentitiesCount; superentityIndex++ )

      {

         for( LocalIndexType i = 0;

              i < mesh.template getSubentitiesCount< SuperdimensionTag::value, SubdimensionTag::value >( superentityIndex );

              i++ )

         {

            const GlobalIndexType subentityIndex = mesh.template getSubentityIndex< SuperdimensionTag::value, SubdimensionTag::value >( superentityIndex, i );

            superentityStorageNetwork.getValues( subentityIndex ).setValue( superentitiesCounts[ subentityIndex ]++, superentityIndex );

            auto row = matrix.getRow( subentityIndex );

            row.setElement( superentitiesCounts[ subentityIndex ]++, superentityIndex, true );

         }

      }
@@ -205,22 +213,24 @@ class EntityInitializerLayer< MeshConfig, 
   using SuperentityTraitsType      = typename MeshTraits< MeshConfig >::template EntityTraits< SuperdimensionTag::value >;

   using SuperentityTopology        = typename SuperentityTraitsType::EntityTopology;

   using SubentitySeedsCreatorType  = SubentitySeedsCreator< MeshConfig, SuperdimensionTag, SubdimensionTag >;

   using SuperentityStorageNetwork  = typename MeshTraits< MeshConfig >::template SuperentityTraits< SubentityTopology, SuperdimensionTag::value >::StorageNetworkType;

   using ValuesAllocationVectorType = typename SuperentityStorageNetwork::ValuesAllocationVectorType;

   using SuperentityMatrixType      = typename MeshTraits< MeshConfig >::SuperentityMatrixType;



public:

   static void initSuperentities( InitializerType& meshInitializer, MeshType& mesh )

   {

      //std::cout << "   Initiating superentities with dimension " << SuperdimensionTag::value << " for subentities with dimension " << SubdimensionTag::value << " ... " << std::endl;



      // counter for superentities of each subentity

      const GlobalIndexType subentitiesCount = mesh.template getEntitiesCount< SubdimensionTag::value >();

      ValuesAllocationVectorType superentitiesCounts( subentitiesCount );

      const GlobalIndexType superentitiesCount = mesh.template getEntitiesCount< SuperdimensionTag::value >();

      if( SubdimensionTag::value > 0 )

         meshInitializer.template initSubentityMatrix< SuperdimensionTag::value, SubdimensionTag::value >( superentitiesCount, subentitiesCount );



      // counter for superentities of each subentity

      auto& superentitiesCounts = meshInitializer.template getSuperentitiesCountsArray< SubdimensionTag::value, SuperdimensionTag::value >();

      superentitiesCounts.setSize( subentitiesCount );

      superentitiesCounts.setValue( 0 );



      for( GlobalIndexType superentityIndex = 0;

           superentityIndex < mesh.template getEntitiesCount< SuperdimensionTag::value >();

           superentityIndex++ )

      for( GlobalIndexType superentityIndex = 0; superentityIndex < superentitiesCount; superentityIndex++ )

      {

         auto subentitySeeds = SubentitySeedsCreatorType::create( meshInitializer.template getSubvertices< SuperdimensionTag::value >( superentityIndex ) );

         auto& subentityOrientationsArray = meshInitializer.template subentityOrientationsArray< SuperdimensionTag::value, SubdimensionTag::value >( superentityIndex );
@@ -234,21 +244,21 @@ public: 
      }



      // allocate superentities storage

      SuperentityStorageNetwork& superentityStorageNetwork = meshInitializer.template meshSuperentityStorageNetwork< SubdimensionTag::value, SuperdimensionTag::value >();

      superentityStorageNetwork.allocate( superentitiesCounts );

      SuperentityMatrixType& matrix = meshInitializer.template getSuperentitiesMatrix< SubdimensionTag::value, SuperdimensionTag::value >();

      matrix.setDimensions( subentitiesCount, superentitiesCount );

      matrix.setRowCapacities( superentitiesCounts );

      superentitiesCounts.setValue( 0 );



      // initialize superentities storage

      for( GlobalIndexType superentityIndex = 0;

           superentityIndex < mesh.template getEntitiesCount< SuperdimensionTag::value >();

           superentityIndex++ )

      for( GlobalIndexType superentityIndex = 0; superentityIndex < superentitiesCount; superentityIndex++ )

      {

         for( LocalIndexType i = 0;

              i < mesh.template getSubentitiesCount< SuperdimensionTag::value, SubdimensionTag::value >( superentityIndex );

              i++ )

         {

            const GlobalIndexType subentityIndex = mesh.template getSubentityIndex< SuperdimensionTag::value, SubdimensionTag::value >( superentityIndex, i );

            superentityStorageNetwork.getValues( subentityIndex ).setValue( superentitiesCounts[ subentityIndex ]++, superentityIndex );

            auto row = matrix.getRow( subentityIndex );

            row.setElement( superentitiesCounts[ subentityIndex ]++, superentityIndex, true );

         }

      }
@@ -292,6 +302,12 @@ public: 
   static void initSuperentities( InitializerType& meshInitializer, MeshType& mesh )

   {

      //std::cout << "   Initiating superentities with dimension " << SuperdimensionTag::value << " for subentities with dimension " << SubdimensionTag::value << " ... " << std::endl;



      const GlobalIndexType subentitiesCount = mesh.template getEntitiesCount< SubdimensionTag::value >();

      const GlobalIndexType superentitiesCount = mesh.template getEntitiesCount< SuperdimensionTag::value >();

      if( SubdimensionTag::value > 0 )

         meshInitializer.template initSubentityMatrix< SuperdimensionTag::value, SubdimensionTag::value >( superentitiesCount, subentitiesCount );



      for( GlobalIndexType superentityIndex = 0;

           superentityIndex < mesh.template getEntitiesCount< SuperdimensionTag::value >();

           superentityIndex++ )
@@ -349,6 +365,12 @@ public: 
   static void initSuperentities( InitializerType& meshInitializer, MeshType& mesh )

   {

      //std::cout << "   Initiating superentities with dimension " << SuperdimensionTag::value << " for subentities with dimension " << SubdimensionTag::value << " ... " << std::endl;



      const GlobalIndexType subentitiesCount = mesh.template getEntitiesCount< SubdimensionTag::value >();

      const GlobalIndexType superentitiesCount = mesh.template getEntitiesCount< SuperdimensionTag::value >();

      if( SubdimensionTag::value > 0 )

         meshInitializer.template initSubentityMatrix< SuperdimensionTag::value, SubdimensionTag::value >( superentitiesCount, subentitiesCount );



      for( GlobalIndexType superentityIndex = 0;

           superentityIndex < mesh.template getEntitiesCount< SuperdimensionTag::value >();

           superentityIndex++ )
@@ -399,22 +421,22 @@ class EntityInitializerLayer< MeshConfig, 
   using SuperentityTraitsType      = typename MeshTraits< MeshConfig >::template EntityTraits< SuperdimensionTag::value >;

   using SuperentityTopology        = typename SuperentityTraitsType::EntityTopology;

   using SubentitySeedsCreatorType  = SubentitySeedsCreator< MeshConfig, SuperdimensionTag, SubdimensionTag >;

   using SuperentityStorageNetwork  = typename MeshTraits< MeshConfig >::template SuperentityTraits< SubentityTopology, SuperdimensionTag::value >::StorageNetworkType;

   using ValuesAllocationVectorType = typename SuperentityStorageNetwork::ValuesAllocationVectorType;

   using SuperentityMatrixType      = typename MeshTraits< MeshConfig >::SuperentityMatrixType;



public:

   static void initSuperentities( InitializerType& meshInitializer, MeshType& mesh )

   {

      //std::cout << "   Initiating superentities with dimension " << SuperdimensionTag::value << " for subentities with dimension " << SubdimensionTag::value << " ... " << std::endl;



      // counter for superentities of each subentity

      const GlobalIndexType subentitiesCount = mesh.template getEntitiesCount< SubdimensionTag::value >();

      ValuesAllocationVectorType superentitiesCounts( subentitiesCount );

      const GlobalIndexType superentitiesCount = mesh.template getEntitiesCount< SuperdimensionTag::value >();



      // counter for superentities of each subentity

      auto& superentitiesCounts = meshInitializer.template getSuperentitiesCountsArray< SubdimensionTag::value, SuperdimensionTag::value >();

      superentitiesCounts.setSize( subentitiesCount );

      superentitiesCounts.setValue( 0 );



      for( GlobalIndexType superentityIndex = 0;

           superentityIndex < mesh.template getEntitiesCount< SuperdimensionTag::value >();

           superentityIndex++ )

      for( GlobalIndexType superentityIndex = 0; superentityIndex < superentitiesCount; superentityIndex++ )

      {

         auto subentitySeeds = SubentitySeedsCreatorType::create( meshInitializer.template getSubvertices< SuperdimensionTag::value >( superentityIndex ) );

         for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )
@@ -425,20 +447,20 @@ public: 
      }



      // allocate superentities storage

      SuperentityStorageNetwork& superentityStorageNetwork = meshInitializer.template meshSuperentityStorageNetwork< SubdimensionTag::value, SuperdimensionTag::value >();

      superentityStorageNetwork.allocate( superentitiesCounts );

      SuperentityMatrixType& matrix = meshInitializer.template getSuperentitiesMatrix< SubdimensionTag::value, SuperdimensionTag::value >();

      matrix.setDimensions( subentitiesCount, superentitiesCount );

      matrix.setRowCapacities( superentitiesCounts );

      superentitiesCounts.setValue( 0 );



      // initialize superentities storage

      for( GlobalIndexType superentityIndex = 0;

           superentityIndex < mesh.template getEntitiesCount< SuperdimensionTag::value >();

           superentityIndex++ )

      for( GlobalIndexType superentityIndex = 0; superentityIndex < superentitiesCount; superentityIndex++ )

      {

         auto subentitySeeds = SubentitySeedsCreatorType::create( meshInitializer.template getSubvertices< SuperdimensionTag::value >( superentityIndex ) );

         for( LocalIndexType i = 0; i < subentitySeeds.getSize(); i++ )

         {

            const GlobalIndexType subentityIndex = meshInitializer.findEntitySeedIndex( subentitySeeds[ i ] );

            superentityStorageNetwork.getValues( subentityIndex ).setValue( superentitiesCounts[ subentityIndex ]++, superentityIndex );

            auto row = matrix.getRow( subentityIndex );

            row.setElement( superentitiesCounts[ subentityIndex ]++, superentityIndex, true );

         }

      }