Commit 999fd145 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Merge branch 'JK/mesh-fix' into 'develop' 

Distributed mesh fix

See merge request !74
parents 8ad20c2f 9c64790b

src/TNL/Algorithms/CudaReductionKernel.h

+2 −1
Original line number Diff line number Diff line
@@ -54,7 +54,8 @@ auto CudaReductionFunctorWrapper( Reduction&& reduction, Arg1&& arg1, Arg2&& arg 
// let's suppress the aforementioned warning...

#ifdef __NVCC__

#pragma push

#pragma diag_suppress 2979

#pragma diag_suppress 2979  // error number for nvcc 10.2

#pragma diag_suppress 3123  // error number for nvcc 11.1

#endif

   return std::forward<Reduction>(reduction)( std::forward<Arg1>(arg1), std::forward<Arg2>(arg2) );

#ifdef __NVCC__

src/TNL/Meshes/DistributedMeshes/DistributedMesh.h

+3 −0
Original line number Diff line number Diff line
@@ -197,11 +197,14 @@ public: 
                << "\tMesh dimension:\t" << getMeshDimension() << "\n"

                << "\tCell topology:\t" << getType( typename Cell::EntityTopology{} ) << "\n"

                << "\tCells count:\t" << cellsCount << "\n"

                << "\tFaces count:\t" << localMesh.template getEntitiesCount< Mesh::getMeshDimension() - 1 >() << "\n"

                << "\tvertices count:\t" << verticesCount << "\n"

                << "\tGhost levels:\t" << getGhostLevels() << "\n"

                << "\tGhost cells count:\t" << localMesh.template getGhostEntitiesCount< Mesh::getMeshDimension() >() << "\n"

                << "\tGhost faces count:\t" << localMesh.template getGhostEntitiesCount< Mesh::getMeshDimension() - 1 >() << "\n"

                << "\tGhost vertices count:\t" << localMesh.template getGhostEntitiesCount< 0 >() << "\n"

                << "\tBoundary cells count:\t" << localMesh.template getBoundaryIndices< Mesh::getMeshDimension() >().getSize() << "\n"

                << "\tBoundary faces count:\t" << localMesh.template getBoundaryIndices< Mesh::getMeshDimension() - 1 >().getSize() << "\n"

                << "\tBoundary vertices count:\t" << localMesh.template getBoundaryIndices< 0 >().getSize() << "\n";

            const GlobalIndexType globalPointIndices = getGlobalIndices< 0 >().getSize();

            const GlobalIndexType globalCellIndices = getGlobalIndices< Mesh::getMeshDimension() >().getSize();

src/TNL/Meshes/DistributedMeshes/distributeSubentities.h

+339 −131

File changed.

Preview size limit exceeded, changes collapsed.

src/Tools/CMakeLists.txt

+15 −16
Original line number Diff line number Diff line
@@ -7,6 +7,11 @@ ADD_EXECUTABLE(tnl-grid-setup tnl-grid-setup.cpp ) 
ADD_EXECUTABLE(tnl-grid-to-mesh tnl-grid-to-mesh.cpp )

ADD_EXECUTABLE(tnl-mesh-converter tnl-mesh-converter.cpp )

ADD_EXECUTABLE(tnl-game-of-life tnl-game-of-life.cpp )

if( BUILD_CUDA )

   CUDA_ADD_EXECUTABLE(tnl-test-distributed-mesh tnl-test-distributed-mesh.cu )

else()

   ADD_EXECUTABLE(tnl-test-distributed-mesh tnl-test-distributed-mesh.cpp )

endif()

ADD_EXECUTABLE(tnl-init tnl-init.cpp )

ADD_EXECUTABLE(tnl-view tnl-view.cpp )

ADD_EXECUTABLE(tnl-diff tnl-diff.cpp )
@@ -27,26 +32,19 @@ endif() 


find_package( ZLIB )

if( ZLIB_FOUND )

   target_compile_definitions(tnl-view PUBLIC "-DHAVE_ZLIB")

   target_include_directories(tnl-view PUBLIC ${ZLIB_INCLUDE_DIRS})

   target_link_libraries(tnl-view ${ZLIB_LIBRARIES})



   target_compile_definitions(tnl-mesh-converter PUBLIC "-DHAVE_ZLIB")

   target_include_directories(tnl-mesh-converter PUBLIC ${ZLIB_INCLUDE_DIRS})

   target_link_libraries(tnl-mesh-converter ${ZLIB_LIBRARIES})



   target_compile_definitions(tnl-game-of-life PUBLIC "-DHAVE_ZLIB")

   target_include_directories(tnl-game-of-life PUBLIC ${ZLIB_INCLUDE_DIRS})

   target_link_libraries(tnl-game-of-life ${ZLIB_LIBRARIES})

   foreach( target IN ITEMS tnl-view tnl-mesh-converter tnl-game-of-life tnl-test-distributed-mesh )

      target_compile_definitions(${target} PUBLIC "-DHAVE_ZLIB")

      target_include_directories(${target} PUBLIC ${ZLIB_INCLUDE_DIRS})

      target_link_libraries(${target} ${ZLIB_LIBRARIES})

   endforeach()

endif()



find_package( tinyxml2 QUIET )

if( tinyxml2_FOUND )

   target_compile_definitions(tnl-mesh-converter PUBLIC "-DHAVE_TINYXML2")

   target_link_libraries(tnl-mesh-converter tinyxml2::tinyxml2)



   target_compile_definitions(tnl-game-of-life PUBLIC "-DHAVE_TINYXML2")

   target_link_libraries(tnl-game-of-life tinyxml2::tinyxml2)

   foreach( target IN ITEMS tnl-mesh-converter tnl-game-of-life tnl-test-distributed-mesh )

      target_compile_definitions(${target} PUBLIC "-DHAVE_TINYXML2")

      target_link_libraries(${target} tinyxml2::tinyxml2)

   endforeach()

endif()



find_package( METIS QUIET )
@@ -80,6 +78,7 @@ INSTALL( TARGETS tnl-init 
                 tnl-grid-to-mesh

                 tnl-mesh-converter

                 tnl-game-of-life

                 tnl-test-distributed-mesh

                 tnl-dicom-reader

                 tnl-image-converter

                 tnl-lattice-init

src/Tools/tnl-decompose-mesh.cpp

+78 −28
Original line number Diff line number Diff line
@@ -14,9 +14,11 @@ 
#include <TNL/Meshes/TypeResolver/TypeResolver.h>

#include <TNL/Meshes/Writers/VTUWriter.h>

#include <TNL/Meshes/Writers/PVTUWriter.h>

#include <TNL/Meshes/MeshDetails/IndexPermutationApplier.h>



#include <metis.h>



#include <numeric>   // std::iota

#include <memory>

#include <vector>

#include <set>
@@ -392,6 +394,12 @@ struct DecomposeMesh 
      // set METIS options from parameters

      setMETISoptions(options, parameters);



      if( nparts == 1 ) {

         // k-way partitioning from Metis fails for nparts == 1 (segfault),

         // RB succeeds but produces nonsense

         part_array.setValue( 0 );

      }

      else {

         if( options[METIS_OPTION_PTYPE] == METIS_PTYPE_KWAY ) {

            std::cout << "Running METIS_PartGraphKway..." << std::endl;

            status = METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, vsize, adjwgt, &nparts, tpwgts, ubvec, options, &objval, part);
@@ -412,6 +420,7 @@ struct DecomposeMesh 
            default:

               throw std::runtime_error( "METIS_PartGraph failed with an unspecified error." );

         }

      }



      // deallocate auxiliary vectors

      connectivity.clear();
@@ -508,8 +517,8 @@ struct DecomposeMesh 
      points_counts.setValue( 0 );

      {

         // first assign points to subdomains - the subdomain with the highest number takes the point

         // (go over local cells, set subvertex owner = cell owner, higher rank will overwrite)

         IndexArray point_to_subdomain( pointsCount );

         point_to_subdomain.setValue( 0 );

         for( unsigned p = 0; p < nparts; p++ ) {

            for( Index local_idx = 0; local_idx < cells_counts[ p ]; local_idx++ ) {

               const Index global_idx = seed_to_cell_index[ cells_offsets[ p ] + local_idx ];
@@ -517,7 +526,6 @@ struct DecomposeMesh 
               const Index subvertices = cell.template getSubentitiesCount< 0 >();

               for( Index j = 0; j < subvertices; j++ ) {

                  const Index v = cell.template getSubentityIndex< 0 >( j );

                  if( (Index) p > point_to_subdomain[ v ] )

                  point_to_subdomain[ v ] = p;

               }

            }
@@ -623,7 +631,7 @@ struct DecomposeMesh 
         }



         // collect seed indices of ghost cells

         std::vector< Index > ghost_seed_indices;

         std::set< Index > ghost_seed_indices;

         for( unsigned gl = 0; gl < ghost_levels; gl++ ) {

            std::vector< Index > new_ghosts;

            for( auto global_idx : ghost_neighbors ) {
@@ -631,9 +639,15 @@ struct DecomposeMesh 
               const Index neighbors_end = dual_xadj.get()[ global_idx + 1 ];

               for( Index i = neighbors_start; i < neighbors_end; i++ ) {

                  const Index neighbor_idx = dual_adjncy.get()[ i ];

                  new_ghosts.push_back( neighbor_idx );

                  // skip neighbors on the local subdomain

                  if( part[ neighbor_idx ] == (int) p )

                     continue;

                  const Index neighbor_seed_idx = cell_to_seed_index[ neighbor_idx ];

                  ghost_seed_indices.push_back( neighbor_seed_idx );

                  // skip neighbors whose seed was already added

                  if( ghost_seed_indices.count( neighbor_seed_idx ) == 0 ) {

                     new_ghosts.push_back( neighbor_idx );

                     ghost_seed_indices.insert( neighbor_seed_idx );

                  }

               }

            }

            std::swap( ghost_neighbors, new_ghosts );
@@ -641,10 +655,7 @@ struct DecomposeMesh 
         ghost_neighbors.clear();

         ghost_neighbors.shrink_to_fit();



         // sort ghosts by their seed index (i.g. global index on the decomposed mesh)

         std::sort( ghost_seed_indices.begin(), ghost_seed_indices.end() );



         // add ghost cells

         // add ghost cells (the set is sorted by the seed index)

         for( auto idx : ghost_seed_indices ) {

            // the ghost_seed_indices array may contain duplicates and even local

            // cells, but add_cell takes care of uniqueness, so we don't have to
@@ -653,7 +664,6 @@ struct DecomposeMesh 
            add_cell( cell );

         }

         ghost_seed_indices.clear();

         ghost_seed_indices.shrink_to_fit();

         cell_global_indices.clear();



         // set points needed for the subdomain
@@ -685,14 +695,54 @@ struct DecomposeMesh 
         for( Index i = cells_counts[ p ]; i < cellSeeds.getSize(); i++ )

            cellGhosts[ i ] = (std::uint8_t) Meshes::VTK::CellGhostTypes::DUPLICATECELL;

         // point ghosts are more tricky because they were assigned to the subdomain with higher number

         Index pointsGhostCount = 0;

         for( Index i = 0; i < points.getSize(); i++ ) {

            const Index global_idx = pointsGlobalIndices[ i ];

            if( global_idx < points_offsets[ p ] || global_idx >= points_offsets[ p ] + points_counts[ p ] )

            if( global_idx < points_offsets[ p ] || global_idx >= points_offsets[ p ] + points_counts[ p ] ) {

               pointGhosts[ i ] = (std::uint8_t) Meshes::VTK::PointGhostTypes::DUPLICATEPOINT;

               pointsGhostCount++;

            }

            else

               pointGhosts[ i ] = 0;

         }



         // reorder ghost points to make sure that global indices are sorted

         {

            // prepare vector with an identity permutation

            std::vector< Index > permutation( points.getSize() );

            std::iota( permutation.begin(), permutation.end(), (Index) 0 );



            // sort the subarray corresponding to ghost entities by the global index

            std::stable_sort( permutation.begin() + points.getSize() - pointsGhostCount,

                              permutation.end(),

                              [&pointsGlobalIndices](auto& left, auto& right) {

               return pointsGlobalIndices[ left ] < pointsGlobalIndices[ right ];

            });



            // copy the permutation into TNL array

            typename Mesh::GlobalIndexArray perm( permutation );

            permutation.clear();

            permutation.shrink_to_fit();



            // apply the permutation

            using PermutationApplier = TNL::Meshes::IndexPermutationApplier< Mesh, 0 >;

            // - pointGhosts

            PermutationApplier::permuteArray( pointGhosts, perm );

            // - pointsGlobalIndices

            PermutationApplier::permuteArray( pointsGlobalIndices, perm );

            // - points

            PermutationApplier::permuteArray( points, perm );

            // - cellSeeds.setCornerID (inverse perm)

            std::vector< Index > iperm( points.getSize() );

            for( Index i = 0; i < perm.getSize(); i++ )

               iperm[ perm[ i ] ] = i;

            for( Index i = 0; i < cellSeeds.getSize(); i++ ) {

               auto& cornerIds = cellSeeds[ i ].getCornerIds();

               for( Index v = 0; v < cornerIds.getSize(); v++ )

                  cornerIds[ v ] = iperm[ cornerIds[ v ] ];

            }

         }



         // init mesh for the subdomain

         Mesh subdomain;

         subdomain.init( points, cellSeeds );