diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 1f21a300d52..ef1e9f1d009 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -66,6 +66,7 @@ #include #include #include +#include #ifdef GEOS_USE_MPI #include @@ -820,7 +821,381 @@ scatterByBlock( vtkDataSet & mesh ) return result; } +/** + * @brief Separate 2D and 3D cells from a mesh + * + * Extracts cells by topological dimension into separate unstructured grids while + * maintaining mappings to original cell indices. 1D and 0D cells are ignored. + * + * @param[in] mesh Input mesh containing cells of mixed dimensions + * @param[out] cells3D Unstructured grid containing only 3D volumetric cells + * @param[out] cells2D Unstructured grid containing only 2D surface cells + * @param[out] cells3DToOriginal Index mapping from extracted 3D cells to original mesh + * @param[out] cells2DToOriginal Index mapping from extracted 2D cells to original mesh + * + * @note The output grids are deep copies and independent of the input mesh + */ +static void separateCellsByDimension( vtkDataSet & mesh, + vtkSmartPointer< vtkUnstructuredGrid > & cells3D, + vtkSmartPointer< vtkUnstructuredGrid > & cells2D, + array1d< vtkIdType > & cells3DToOriginal, + array1d< vtkIdType > & cells2DToOriginal ) +{ + GEOS_MARK_FUNCTION; + + vtkIdType const numCells = mesh.GetNumberOfCells(); + + vtkNew< vtkIdList > indices3D; + vtkNew< vtkIdList > indices2D; + + // Classify cells by dimension + for( vtkIdType i = 0; i < numCells; ++i ) + { + int const cellDim = mesh.GetCell( i )->GetCellDimension(); + if( cellDim == 3 ) + { + indices3D->InsertNextId( i ); + } + else if( cellDim == 2 ) + { + indices2D->InsertNextId( i ); + } + // Note: 1D edges and 0D vertices are intentionally ignored + } + + GEOS_LOG_RANK_0( GEOS_FMT( "Separating mesh: {} 3D cells, {} 2D cells (from {} total)", + indices3D->GetNumberOfIds(), + indices2D->GetNumberOfIds(), + numCells ) ); + + // Extract and deep-copy 3D cells + vtkNew< vtkExtractCells > extractor3D; + extractor3D->SetInputDataObject( &mesh ); + extractor3D->SetExtractAllCells( false ); + extractor3D->SetCellList( indices3D ); + extractor3D->Update(); + + cells3D = vtkSmartPointer< vtkUnstructuredGrid >::New(); + cells3D->DeepCopy( extractor3D->GetOutput() ); + + // Store index mapping for 3D cells + cells3DToOriginal.resize( indices3D->GetNumberOfIds() ); + for( vtkIdType i = 0; i < indices3D->GetNumberOfIds(); ++i ) + { + cells3DToOriginal[i] = indices3D->GetId( i ); + } + + // Extract and deep-copy 2D cells + vtkNew< vtkExtractCells > extractor2D; + extractor2D->SetInputDataObject( &mesh ); + extractor2D->SetExtractAllCells( false ); + extractor2D->SetCellList( indices2D ); + extractor2D->Update(); + + cells2D = vtkSmartPointer< vtkUnstructuredGrid >::New(); + cells2D->DeepCopy( extractor2D->GetOutput() ); + + // Store index mapping for 2D cells + cells2DToOriginal.resize( indices2D->GetNumberOfIds() ); + for( vtkIdType i = 0; i < indices2D->GetNumberOfIds(); ++i ) + { + cells2DToOriginal[i] = indices2D->GetId( i ); + } +} + +/** + * @brief Build mapping from 2D cells to their neighboring 3D cells using global cell IDs + * @param[in] mesh Original mesh containing both 2D and 3D cells with global IDs + * @param[in] cells2DToOriginal Mapping from local 2D indices to original mesh indices + * @param[in] cells3DToOriginal Mapping from local 3D indices to original mesh indices + * @return ArrayOfArrays mapping each 2D cell index to global IDs of neighboring 3D cells + */ +static ArrayOfArrays< vtkIdType, int64_t > +build2DTo3DNeighbors( vtkDataSet & mesh, + arrayView1d< vtkIdType const > cells2DToOriginal, + arrayView1d< vtkIdType const > cells3DToOriginal ) +{ + GEOS_MARK_FUNCTION; + + // Retrieve global cell ID array + vtkDataArray * globalCellIds = mesh.GetCellData()->GetGlobalIds(); + GEOS_ERROR_IF( globalCellIds == nullptr, + "Global cell IDs must be present in mesh for 2D-3D neighbor mapping" ); + + +// Build reverse lookup: original mesh index -> global cell ID (3D cells only) + stdUnorderedMap< vtkIdType, int64_t > original3DToGlobalId; + original3DToGlobalId.reserve( cells3DToOriginal.size() ); + + for( localIndex i = 0; i < cells3DToOriginal.size(); ++i ) + { + vtkIdType const origIdx = cells3DToOriginal[i]; + int64_t const globalId = static_cast< int64_t >( globalCellIds->GetTuple1( origIdx ) ); + original3DToGlobalId.emplace( origIdx, globalId ); // ← FIXED + } + + ArrayOfArrays< vtkIdType, int64_t > neighbors2Dto3D; + neighbors2Dto3D.reserve( cells2DToOriginal.size() ); + + // Topology statistics + localIndex numStandalone = 0; + localIndex numBoundary = 0; // 1 neighbor + localIndex numInternal = 0; // 2 neighbors + localIndex numJunction = 0; // >2 neighbors + + // Build neighbor list for each 2D cell + for( localIndex i = 0; i < cells2DToOriginal.size(); ++i ) + { + vtkIdType const origIdx2D = cells2DToOriginal[i]; + vtkCell * cell2D = mesh.GetCell( origIdx2D ); + vtkIdList * pointIds2D = cell2D->GetPointIds(); + + // Find all cells sharing ALL nodes with this 2D cell (exact face match) + vtkNew< vtkIdList > neighborCells; + mesh.GetCellNeighbors( origIdx2D, pointIds2D, neighborCells ); + + // Filter for 3D neighbors and retrieve their global IDs + array1d< int64_t > neighbor3DGlobalIds; + neighbor3DGlobalIds.reserve( neighborCells->GetNumberOfIds() ); + + for( vtkIdType n = 0; n < neighborCells->GetNumberOfIds(); ++n ) + { + vtkIdType const neighborIdx = neighborCells->GetId( n ); + + // Only consider 3D cells + if( mesh.GetCell( neighborIdx )->GetCellDimension() == 3 ) + { + auto it = original3DToGlobalId.find( neighborIdx ); + GEOS_ERROR_IF( it == original3DToGlobalId.end(), + GEOS_FMT( "3D neighbor at index {} not found in 3D cell mapping", neighborIdx ) ); + + neighbor3DGlobalIds.emplace_back( it->second ); + } + } + + // Update topology statistics + localIndex const numNeighbors = neighbor3DGlobalIds.size(); + switch( numNeighbors ) + { + case 0: numStandalone++; break; + case 1: numBoundary++; break; + case 2: numInternal++; break; + default: numJunction++; break; + } + + neighbors2Dto3D.appendArray( neighbor3DGlobalIds.begin(), neighbor3DGlobalIds.end() ); + } + + // Print diagnostic summary + GEOS_LOG_RANK_0( "\n2D-to-3D Neighbor Topology " ); + GEOS_LOG_RANK_0( GEOS_FMT( " Total 2D cells: {}", cells2DToOriginal.size() ) ); + GEOS_LOG_RANK_0( GEOS_FMT( " Standalone (0 neighbors): {}", numStandalone ) ); + GEOS_LOG_RANK_0( GEOS_FMT( " Boundary (1 neighbor): {}", numBoundary ) ); + GEOS_LOG_RANK_0( GEOS_FMT( " Internal (2 neighbors): {}", numInternal ) ); + GEOS_LOG_RANK_0( GEOS_FMT( " Junction (>2 neighbors):{}", numJunction ) ); + + // Standalone 2D cells indicate mesh topology errors + GEOS_ERROR_IF( numStandalone > 0, + GEOS_FMT( "{} orphaned 2D cells detected with no 3D neighbors. " + "These may be artifacts or detached surfaces. " + "Please clean the mesh or verify the geometry.", numStandalone ) ); + + return neighbors2Dto3D; +} + +/** + * @brief Assign partition ownership for 2D cells based on their 3D neighbor locations + * + * Uses a deterministic tie-breaking rule (minimum global ID) to assign 2D cells that + * have multiple 3D neighbors (e.g., internal fractures). + * + * @param[in] neighbors2Dto3D Neighbor map from build2DTo3DNeighbors() + * @param[in] globalIdTo3DPartition Complete map of 3D cell global ID -> rank assignment + * @param[in] comm MPI communicator (unused, kept for API consistency) + * @return Partition (rank) assignment for each 2D cell + */ +static array1d< int64_t > +assign2DCellsTo3DPartitions( ArrayOfArrays< vtkIdType, int64_t > const & neighbors2Dto3D, + stdUnorderedMap< int64_t, int64_t > const & globalIdTo3DPartition ) +{ + GEOS_MARK_FUNCTION; + + array1d< int64_t > partitions2D( neighbors2Dto3D.size() ); + + for( localIndex i = 0; i < neighbors2Dto3D.size(); ++i ) + { + auto neighbors = neighbors2Dto3D[i]; + + GEOS_ASSERT_MSG( neighbors.size() > 0, + "2D cell should have at least one neighbor (enforced by build2DTo3DNeighbors)" ); + + // Deterministic tie-breaking: choose neighbor with minimum global ID + int64_t minGlobalId = neighbors[0]; + for( localIndex n = 1; n < neighbors.size(); ++n ) + { + if( neighbors[n] < minGlobalId ) + { + minGlobalId = neighbors[n]; + } + } + + // Look up partition assignment from the complete 3D partition map + auto it = globalIdTo3DPartition.find( minGlobalId ); + GEOS_ERROR_IF( it == globalIdTo3DPartition.end(), + GEOS_FMT( "3D neighbor with global ID {} not found in partition map", minGlobalId ) ); + + partitions2D[i] = it->second; + } + + return partitions2D; +} + +/** + * @brief Merge 2D surface cells with redistributed 3D volume cells + * + * This function completes the 2D/3D redistribution workflow: + * 1. Builds a global 3D partition map (globalID -> rank) + * 2. Assigns 2D cells to ranks based on their 3D neighbor locations (rank 0 only) + * 3. Redistributes 2D cells to appropriate ranks + * 4. Merges local 2D and 3D cells on each rank into a unified mesh + * + * @param[in] redistributed3D Already partitioned 3D cells with global IDs + * @param[in] cells2D 2D boundary cells (still on rank 0 or original distribution) + * @param[in] neighbors2Dto3D Pre-computed 2D-to-3D neighbor mapping + * @param[in] redistributedFractures Already partitioned fracture meshes (pass-through) + * @param[in] comm MPI communicator + * @return Complete AllMeshes object with merged main mesh and fractures + */ +static AllMeshes +merge2D3DCellsAndRedistribute( vtkSmartPointer< vtkDataSet > redistributed3D, + vtkSmartPointer< vtkUnstructuredGrid > cells2D, + ArrayOfArrays< vtkIdType, int64_t > const & neighbors2Dto3D, + stdMap< string, vtkSmartPointer< vtkDataSet > > const & redistributedFractures, + MPI_Comm const comm ) +{ + GEOS_MARK_FUNCTION; + + int const rank = MpiWrapper::commRank( comm ); + int const numRanks = MpiWrapper::commSize( comm ); + + // Step 1: Build complete 3D partition map (all ranks participate) + vtkDataArray * redistributed3DGlobalIds = redistributed3D->GetCellData()->GetGlobalIds(); + GEOS_ERROR_IF( redistributed3DGlobalIds == nullptr, + "Global cell IDs required in redistributed 3D mesh for 2D cell assignment" ); + + vtkIdType const local3DCells = redistributed3D->GetNumberOfCells(); + + // Gather counts from all ranks + array1d< int > cellCounts; + MpiWrapper::allGather( static_cast< int >( local3DCells ), cellCounts, comm ); + + // Compute offsets + array1d< int > offsets( numRanks ); + int totalCells = 0; + for( int r = 0; r < numRanks; ++r ) + { + offsets[r] = totalCells; + totalCells += cellCounts[r]; + } + + // Gather local global IDs + array1d< int64_t > localGlobalIds( local3DCells ); + for( vtkIdType i = 0; i < local3DCells; ++i ) + { + localGlobalIds[i] = static_cast< int64_t >( redistributed3DGlobalIds->GetTuple1( i ) ); + } + + // All-gather global IDs to all ranks + array1d< int64_t > allGlobalIds( totalCells ); + MpiWrapper::allgatherv( localGlobalIds.data(), static_cast< int >( local3DCells ), + allGlobalIds.data(), cellCounts.data(), offsets.data(), + comm ); + +// Build complete partition map: global cell ID -> owning rank + stdUnorderedMap< int64_t, int64_t > complete3DPartitionMap; + complete3DPartitionMap.reserve( totalCells ); + + for( int r = 0; r < numRanks; ++r ) + { + for( int i = 0; i < cellCounts[r]; ++i ) + { + int64_t const globalId = allGlobalIds[offsets[r] + i]; + complete3DPartitionMap.emplace( globalId, r ); // ← FIXED + } + } + // Step 2: Rank 0 assigns 2D partitions and splits the mesh + vtkSmartPointer< vtkPartitionedDataSet > split2DCells = vtkSmartPointer< vtkPartitionedDataSet >::New(); + vtkIdType expected2DCells = 0; + + if( rank == 0 ) + { + vtkIdType const numCells2D = cells2D->GetNumberOfCells(); + expected2DCells = numCells2D; + + GEOS_LOG_RANK_0( GEOS_FMT( "Assigning {} 2D cells to partitions", numCells2D ) ); + + // Use the helper function instead of inlining + array1d< int64_t > partitions2D = assign2DCellsTo3DPartitions( + neighbors2Dto3D, + complete3DPartitionMap ); + + split2DCells = splitMeshByPartition( cells2D, numRanks, partitions2D.toViewConst() ); + } + else + { + // Non-root ranks: create empty partitioned dataset with correct structure + // Required for vtk::redistribute to work correctly + split2DCells->SetNumberOfPartitions( numRanks ); + for( int r = 0; r < numRanks; ++r ) + { + vtkNew< vtkUnstructuredGrid > emptyPart; + split2DCells->SetPartition( r, emptyPart ); + } + } + + // Step 3: Redistribute 2D cells to target ranks + vtkSmartPointer< vtkUnstructuredGrid > local2DCells = vtk::redistribute( *split2DCells, comm ); + + // Conservation check: verify no 2D cells were lost during redistribution + vtkIdType const total2DCells = MpiWrapper::sum( local2DCells->GetNumberOfCells(), comm ); + MpiWrapper::broadcast( expected2DCells, 0, comm ); + + GEOS_ERROR_IF( total2DCells != expected2DCells, + GEOS_FMT( "2D cell redistribution failed: expected {} cells, got {} cells", + expected2DCells, total2DCells ) ); + + + // Step 4: Merge local 3D and 2D cells on each rank + vtkSmartPointer< vtkUnstructuredGrid > mergedMesh = vtkSmartPointer< vtkUnstructuredGrid >::New(); + + if( local2DCells->GetNumberOfCells() > 0 ) + { + // Merge 3D and 2D cells using VTK append filter + vtkNew< vtkAppendFilter > appendFilter; + appendFilter->AddInputData( redistributed3D ); + appendFilter->AddInputData( local2DCells ); + appendFilter->MergePointsOn(); // Ensures shared nodes are not duplicated + appendFilter->Update(); + + mergedMesh->DeepCopy( appendFilter->GetOutput() ); + + //GEOS_LOG_RANK( GEOS_FMT( "Rank merged {} 3D cells + {} 2D cells = {} total", + // redistributed3D->GetNumberOfCells(), + // local2DCells->GetNumberOfCells(), + // mergedMesh->GetNumberOfCells() ) ); + } + else + { + mergedMesh->DeepCopy( redistributed3D ); + + //GEOS_LOG_RANK( GEOS_FMT( "Rank has {} 3D cells only (no 2D cells)", + // redistributed3D->GetNumberOfCells() ) ); + + } + + return AllMeshes( mergedMesh, redistributedFractures ); +} vtkSmartPointer< vtkUnstructuredGrid > threshold( vtkDataSet & mesh, @@ -1192,7 +1567,6 @@ ensureNoEmptyRank( vtkSmartPointer< vtkDataSet > mesh, return vtk::redistribute( *splitMesh, MPI_COMM_GEOS ); } - AllMeshes redistributeMeshes( integer const logLevel, vtkSmartPointer< vtkDataSet > loadedMesh, @@ -1223,51 +1597,134 @@ redistributeMeshes( integer const logLevel, } } - // Determine if redistribution is required - vtkIdType const minCellsOnAnyRank = MpiWrapper::min( mesh->GetNumberOfCells(), comm ); - if( minCellsOnAnyRank == 0 ) + + // Step 1: Separate 2D and 3D cells from main mesh + vtkSmartPointer< vtkUnstructuredGrid > cells3D, cells2D; + array1d< vtkIdType > cells3DToOriginal, cells2DToOriginal; + separateCellsByDimension( *mesh, cells3D, cells2D, cells3DToOriginal, cells2DToOriginal ); + + // Build 2D-to-3D neighbor mapping using global IDs + ArrayOfArrays< vtkIdType, int64_t > neighbors2Dto3D; + if( cells2D->GetNumberOfCells() > 0 ) { - // Redistribute the mesh over all ranks using simple octree partitions - mesh = redistributeByKdTree( *mesh ); + neighbors2Dto3D = build2DTo3DNeighbors( *mesh, + cells2DToOriginal.toViewConst(), + cells3DToOriginal.toViewConst() ); } + else + { + neighbors2Dto3D.resize( 0, 0 ); + } + + GEOS_LOG_RANK_0( GEOS_FMT( "Separated main mesh into {} 3D cells and {} 2D cells", cells3D->GetNumberOfCells(), cells2D->GetNumberOfCells() )); + + + // Step 2: Redistribute the 3D cells (+ fractures if using a presplit mesh from meshDoctor) + // Determine if redistribution is required (check 3D cells) + vtkIdType const minCellsOnAnyRank = MpiWrapper::min( cells3D->GetNumberOfCells(), comm ); - // Check if a rank does not have a cell after the redistribution - // If this is the case, we need a fix otherwise the next redistribution will fail - // We expect this function to only be called in some pathological cases - if( MpiWrapper::min( mesh->GetNumberOfCells(), comm ) == 0 ) + AllMeshes result3DAndFractures; + + if( minCellsOnAnyRank == 0 ) { - mesh = ensureNoEmptyRank( mesh, comm ); + // Redistribute the 3D mesh over all ranks using simple octree partitions + vtkSmartPointer< vtkDataSet > redistributed3D = redistributeByKdTree( *cells3D ); + + // Check if a rank does not have a cell after the redistribution + // If this is the case, we need a fix otherwise the next redistribution will fail + // We expect this function to only be called in some pathological cases + if( MpiWrapper::min( redistributed3D->GetNumberOfCells(), comm ) == 0 ) + { + redistributed3D = ensureNoEmptyRank( redistributed3D, comm ); + } + + result3DAndFractures.setMainMesh( redistributed3D ); + result3DAndFractures.setFaceBlocks( namesToFractures ); + } + else + { + result3DAndFractures.setMainMesh( cells3D ); + result3DAndFractures.setFaceBlocks( namesToFractures ); } - AllMeshes result; - // Redistribute the mesh again using higher-quality graph partitioner + // Redistribute again using higher-quality graph partitioner if( !structuredIndexAttributeName.empty() ) { - AllMeshes input( mesh, namesToFractures ); - result = redistributeByAreaGraphAndLayer( input, - method, - structuredIndexAttributeName, - comm, - numPartZ, - partitionRefinement - 1 ); + result3DAndFractures = redistributeByAreaGraphAndLayer( result3DAndFractures, + method, + structuredIndexAttributeName, + comm, + numPartZ, + partitionRefinement - 1 ); } else if( partitionRefinement > 0 ) { - AllMeshes input( mesh, namesToFractures ); - result = redistributeByCellGraph( input, method, comm, partitionRefinement - 1 ); + result3DAndFractures = redistributeByCellGraph( result3DAndFractures, method, comm, partitionRefinement - 1 ); } - else + + // Step 3: Merge 2D cells back with redistributed 3D cells + AllMeshes finalResult = merge2D3DCellsAndRedistribute( + result3DAndFractures.getMainMesh(), + cells2D, + neighbors2Dto3D, + result3DAndFractures.getFaceBlocks(), + comm + ); + + // Step 4: Diagnostics { - result.setMainMesh( mesh ); - result.setFaceBlocks( namesToFractures ); + int const rank = MpiWrapper::commRank( comm ); + int const numRanks = MpiWrapper::commSize( comm ); + + vtkIdType local2DCells = 0; + vtkIdType local3DCells = 0; + + vtkSmartPointer< vtkDataSet > finalMesh = finalResult.getMainMesh(); + for( vtkIdType i = 0; i < finalMesh->GetNumberOfCells(); ++i ) + { + int dim = finalMesh->GetCell( i )->GetCellDimension(); + if( dim == 2 ) + local2DCells++; + else if( dim == 3 ) + local3DCells++; + } + + array1d< vtkIdType > all2D, all3D; + MpiWrapper::allGather( local2DCells, all2D, comm ); + MpiWrapper::allGather( local3DCells, all3D, comm ); + + if( rank == 0 ) + { + GEOS_LOG_RANK_0( "\n------------------------------------------------" ); + GEOS_LOG_RANK_0( "| Rk | 3D Cells | 2D Cells | Total (Main) |" ); + GEOS_LOG_RANK_0( "------------------------------------------------" ); + + vtkIdType sum2D = 0, sum3D = 0; + for( int r = 0; r < numRanks; ++r ) + { + sum2D += all2D[r]; + sum3D += all3D[r]; + + GEOS_LOG_RANK_0( "| " << std::setw( 4 ) << r << " | " + << std::setw( 9 ) << all3D[r] << " | " + << std::setw( 9 ) << all2D[r] << " | " + << std::setw( 13 ) << (all3D[r] + all2D[r]) << " |" ); + } + + GEOS_LOG_RANK_0( "------------------------------------------------" ); + GEOS_LOG_RANK_0( "|Total | " << std::setw( 9 ) << sum3D << " | " + << std::setw( 9 ) << sum2D << " | " + << std::setw( 13 ) << (sum3D + sum2D) << " |" ); + GEOS_LOG_RANK_0( "------------------------------------------------" ); + } } - // Logging some information about the redistribution. + // Step 5: Final logging { string const pattern = "{}: {}"; stdVector< string > messages; - messages.push_back( GEOS_FMT( pattern, "Local mesh size", result.getMainMesh()->GetNumberOfCells() ) ); - for( auto const & [faceName, faceMesh]: result.getFaceBlocks() ) + messages.push_back( GEOS_FMT( pattern, "Local mesh size", finalResult.getMainMesh()->GetNumberOfCells() ) ); + for( auto const & [faceName, faceMesh]: finalResult.getFaceBlocks() ) { messages.push_back( GEOS_FMT( pattern, faceName, faceMesh->GetNumberOfCells() ) ); } @@ -1277,7 +1734,7 @@ redistributeMeshes( integer const logLevel, } } - return result; + return finalResult; } /**