fix: 2D/3D cell co-location in parallel mesh redistribution#3967
Merged
castelletto1 merged 17 commits intodevelopfrom Mar 27, 2026
Merged
fix: 2D/3D cell co-location in parallel mesh redistribution#3967castelletto1 merged 17 commits intodevelopfrom
castelletto1 merged 17 commits intodevelopfrom
Conversation
bd713
added
type: bug
bd713
changed the title
added 2 commits
February 9, 2026 11:25
Member
|
@bd713 why not just merge the 2d and 3d meshes into a single mesh object? |
added 2 commits
February 17, 2026 12:03
castelletto1
added
flag: requires rebaseline
| return result; | ||
| } | ||
|
|
||
| /** |
Collaborator
There was a problem hiding this comment.
Suggested change
| /** | |
| /** | |
| * @brief Extract cells by index and create index mapping to original mesh | |
| * | |
| * @param[in] mesh Source mesh | |
| * @param[in] indices Cell indices to extract | |
| * @param[out] output Deep copy of extracted cells | |
| * @param[out] mapping Mapping from output cell indices to original mesh indices | |
| */ | |
| static void extractCellsByIndices( vtkDataSet & mesh, | |
| vtkIdList * indices, | |
| vtkSmartPointer< vtkUnstructuredGrid > & output, | |
| array1d< vtkIdType > & mapping ) | |
| { | |
| vtkNew< vtkExtractCells > extractor; | |
| extractor->SetInputDataObject( &mesh ); | |
| extractor->SetExtractAllCells( false ); | |
| extractor->SetCellList( indices ); | |
| extractor->Update(); | |
| output = vtkSmartPointer< vtkUnstructuredGrid >::New(); | |
| output->DeepCopy( extractor->GetOutput() ); | |
| mapping.resize( indices->GetNumberOfIds() ); | |
| for( vtkIdType i = 0; i < indices->GetNumberOfIds(); ++i ) | |
| { | |
| mapping[i] = indices->GetId( i ); | |
| } | |
| } | |
| /** |
Adding an helper function to avoid doing the same thing for 3D and 2D cells.
Is there a way to avoid using DeepCopy?
Contributor
Author
There was a problem hiding this comment.
Added helper function.
Switched to shallowCopy and also removed the deep copies that separated 2D and 3D cells. Now rely on indices.
Comment on lines
+883
to
+915
| // 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 ); | ||
| } |
Collaborator
There was a problem hiding this comment.
Suggested change
| // 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 ); | |
| } | |
| // Extract and deep-copy 3D cells | |
| extractCellsByIndices( mesh, indices3D, cells3D, cells3DToOriginal ); | |
| // Extract and deep-copy 2D cells | |
| extractCellsByIndices( mesh, indices2D, cells2D, cells2DToOriginal ); |
| "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) |
Collaborator
There was a problem hiding this comment.
Suggested change
| // Build reverse lookup: original mesh index -> global cell ID (3D cells only) | |
| // Build reverse lookup: original mesh index -> global cell ID (3D cells only) |
Comment on lines
+979
to
+985
| // 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 ) ); | ||
|
|
Collaborator
There was a problem hiding this comment.
Suggested change
| // 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 ) ); | |
| // Check if neighbor is a 3D cell and get its global ID | |
| auto it = original3DToGlobalId.find( neighborIdx ); | |
| if( it != original3DToGlobalId.end() ) // Is 3D cell | |
| { | |
| neighbor3DGlobalIds.emplace_back( it->second ); | |
| } | |
| // Non-3D neighbors (2D/1D/0D) are silently skipped |
Comment on lines
+1106
to
+1120
| // 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 ) ); | ||
| } |
Collaborator
There was a problem hiding this comment.
Since global IDs are int64_t, I would expect offsets to be as well int64_t and not int.
| comm ); | ||
|
|
||
| // Build complete partition map: global cell ID -> owning rank | ||
| stdUnorderedMap< int64_t, int64_t > complete3DPartitionMap; |
Collaborator
There was a problem hiding this comment.
Suggested change
| stdUnorderedMap< int64_t, int64_t > complete3DPartitionMap; | |
| stdUnorderedMap< int64_t, int > complete3DPartitionMap; |
Rank is int
| for( int i = 0; i < cellCounts[r]; ++i ) | ||
| { | ||
| int64_t const globalId = allGlobalIds[offsets[r] + i]; | ||
| complete3DPartitionMap.emplace( globalId, r ); // ← FIXED |
Collaborator
There was a problem hiding this comment.
Suggested change
| complete3DPartitionMap.emplace( globalId, r ); // ← FIXED | |
| complete3DPartitionMap.emplace( globalId, r ); |
Comment on lines
+1122
to
+1126
| // 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 ); |
Collaborator
There was a problem hiding this comment.
This is going to be very expensive for large meshes.
Contributor
Author
There was a problem hiding this comment.
Limited gather operation to only the relevant indices
Comment on lines
+1197
to
+1200
| //GEOS_LOG_RANK( GEOS_FMT( "Rank merged {} 3D cells + {} 2D cells = {} total", | ||
| // redistributed3D->GetNumberOfCells(), | ||
| // local2DCells->GetNumberOfCells(), | ||
| // mergedMesh->GetNumberOfCells() ) ); |
Collaborator
There was a problem hiding this comment.
Suggested change
| //GEOS_LOG_RANK( GEOS_FMT( "Rank merged {} 3D cells + {} 2D cells = {} total", | |
| // redistributed3D->GetNumberOfCells(), | |
| // local2DCells->GetNumberOfCells(), | |
| // mergedMesh->GetNumberOfCells() ) ); |
Remove
Comment on lines
+1032
to
+1033
| assign2DCellsTo3DPartitions( ArrayOfArrays< vtkIdType, int64_t > const & neighbors2Dto3D, | ||
| stdUnorderedMap< int64_t, int64_t > const & globalIdTo3DPartition ) |
Collaborator
There was a problem hiding this comment.
Suggested change
| assign2DCellsTo3DPartitions( ArrayOfArrays< vtkIdType, int64_t > const & neighbors2Dto3D, | |
| stdUnorderedMap< int64_t, int64_t > const & globalIdTo3DPartition ) | |
| assign2DCellsTo3DPartitions( ArrayOfArrays< vtkIdType, int64_t > const & neighbors2Dto3D, | |
| stdUnorderedMap< int64_t, int > const & globalIdTo3DPartition ) |
Comment on lines
+1047
to
+1054
| int64_t minGlobalId = neighbors[0]; | ||
| for( localIndex n = 1; n < neighbors.size(); ++n ) | ||
| { | ||
| if( neighbors[n] < minGlobalId ) | ||
| { | ||
| minGlobalId = neighbors[n]; | ||
| } | ||
| } |
Collaborator
There was a problem hiding this comment.
Suggested change
| int64_t minGlobalId = neighbors[0]; | |
| for( localIndex n = 1; n < neighbors.size(); ++n ) | |
| { | |
| if( neighbors[n] < minGlobalId ) | |
| { | |
| minGlobalId = neighbors[n]; | |
| } | |
| } | |
| int64_t minGlobalId = *std::min_element( neighbors.begin(), neighbors.end() ); |
| * | ||
| * @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) |
Collaborator
There was a problem hiding this comment.
Suggested change
| * @param[in] comm MPI communicator (unused, kept for API consistency) |
Comment on lines
+1043
to
+1044
| GEOS_ASSERT_MSG( neighbors.size() > 0, | ||
| "2D cell should have at least one neighbor (enforced by build2DTo3DNeighbors)" ); |
Collaborator
There was a problem hiding this comment.
ASSERT or ERROR? Is the check needed here? build2DTo3DNeighbors should catch this earlier.
castelletto1
approved these changes
Mar 27, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Problem
The current
redistributeMeshes()partitions the entire mesh (both 2D and 3D cells) all at once using graph partitioning.However, there is no guarantee that a 2D surface element will be assigned to a rank containing its adjacent 3D face.
This can lead to configurations where 2D cells are orphaned on ranks with no access to their 3D neighbor faces.
See screenshot below for a real life illustration.
A side effect is that the nodes of these orphaned 2D surfaces are seen as owned by that rank and will prevent them from being properly ghosted in the 3D mesh.
Solution
In this PR, the redistribution workflow is refactored in 3 steps to ensure co-location of 2D and 3D cells:
NB: When using MeshDoctor splitting, the fracture elements suffer the exact same problem. These will be addressed in a separate PR.
Rebaseline
As the 2D elements are no longer directly included in the graph to be partitioned, the resulting partitions may differ from the baseline.
This notably affect the following tests:
In addition, we now reject 2D orphan cells that resulted from a bug in MeshDoctor (geosPythonPackages fix (PR #219)). This required changing the split mesh for: