fix global to local mapping

src/coreComponents/mesh/generators/CellBlock.hpp

@@ -196,6 +196,12 @@ class CellBlock : public CellBlockABC
   arrayView1d< globalIndex const > localToGlobalMapConstView() const
   { return m_localToGlobalMap.toViewConst(); }
 
+  unordered_map< globalIndex, localIndex > & globalToLocalMap()
+  { return m_globalToLocalMap; }
+
+  unordered_map< globalIndex, localIndex > const & globalToLocalMap() const
+  { return m_globalToLocalMap; }
+
   /**
    * @brief Resize the cell block to hold @p numElements
    * @param numElements The new number of elements.
@@ -242,6 +248,9 @@ class CellBlock : public CellBlockABC
   /// Contains the global index of each object.
   array1d< globalIndex > m_localToGlobalMap;
 
+  /// Map from object global index to the local index.
+  unordered_map< globalIndex, localIndex > m_globalToLocalMap;
+
   /// Name of the properties registered from an external mesh
   string_array m_externalPropertyNames;
 

src/coreComponents/mesh/generators/CellBlockManager.hpp

@@ -225,6 +225,15 @@ class CellBlockManager : public CellBlockManagerABC
    */
   void setGlobalLength( real64 globalLength ) { m_globalLength = globalLength; }
 
+  /**
+   * @brief Get cell block at index @p blockIndex.
+   * @param[in] blockIndex The cell block index.
+   * @return Const reference to the instance.
+   *
+   * @note Mainly useful for iteration purposes.
+   */
+  CellBlock const & getCellBlock( localIndex const blockIndex ) const;
+
 private:
 
   struct viewKeyStruct
@@ -253,15 +262,6 @@ class CellBlockManager : public CellBlockManagerABC
    */
   Group & getLineBlocks();
 
-  /**
-   * @brief Get cell block at index @p blockIndex.
-   * @param[in] blockIndex The cell block index.
-   * @return Const reference to the instance.
-   *
-   * @note Mainly useful for iteration purposes.
-   */
-  CellBlock const & getCellBlock( localIndex const blockIndex ) const;
-
   /**
    * @brief Returns the number of cells blocks
    * @return Number of cell blocks

src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp

@@ -102,6 +102,17 @@ class EmbeddedSurfaceBlock : public EmbeddedSurfaceBlockABC
    */
   void setEmbeddedSurfElemPermeability( array1d< real64 > && _perms );
 
+  unordered_map< globalIndex, localIndex > & globalToLocalMap()
+  { return m_globalToLocalMap; }
+
+  unordered_map< globalIndex, localIndex > const & globalToLocalMap() const
+  { return m_globalToLocalMap; }
+
+  void setGlobalToLocalMap( unordered_map< globalIndex, localIndex > && g2l )
+  {
+    m_globalToLocalMap = g2l;
+  }
+
 private:
 
   localIndex m_numEmbeddedSurfaces;
@@ -113,6 +124,10 @@ class EmbeddedSurfaceBlock : public EmbeddedSurfaceBlockABC
   ArrayOfArrays< real64 > m_embeddedSurfElemWidthVectors;
   array1d< real64 > m_embeddedSurfElemApertures;
   array1d< real64 > m_embeddedSurfElemPermeability;
+
+  /// Map from object global index to the local index.
+  unordered_map< globalIndex, localIndex > m_globalToLocalMap;
+
 };
 
 }

src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp

@@ -22,30 +22,51 @@
 namespace geos::vtk
 {
 
+unordered_map< globalIndex, localIndex > buildGlobalToLocal( vtkIdTypeArray const * edfmMeshCellGlobalIds )
+{
+  unordered_map< globalIndex, localIndex > g2l;
+
+  vtkIdType const numEdfms = edfmMeshCellGlobalIds ? edfmMeshCellGlobalIds->GetNumberOfTuples() : 0;
+
+  for( auto i = 0; i < numEdfms; ++i )
+  {
+    // Note that `numEdfms` is zero if `edfmMeshCellGlobalIds` is 0 too.
+    // This prevents from calling a member on a null instance.
+    g2l[edfmMeshCellGlobalIds->GetValue( i )] =  i;
+  }
+
+  return g2l;
+}
+
 void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName,
                                     vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh,
-                                    vtkSmartPointer< vtkDataSet > GEOS_UNUSED_PARAM( mesh ),
                                     CellBlockManager & cellBlockManager )
 {
+  EmbeddedSurfaceBlock & embeddedSurfBlock = cellBlockManager.registerEmbeddedSurfaceBlock( embeddedSurfaceBlockName );
+  vtkIdType const numEdfmFracs = embeddedSurfaceMesh->GetNumberOfCells();
+
+  vtkIdTypeArray const * edfmMeshCellGlobalIds = vtkIdTypeArray::FastDownCast( embeddedSurfaceMesh->GetCellData()->GetGlobalIds() );
+
+  GEOS_ERROR_IF(numEdfmFracs > 0 && !edfmMeshCellGlobalIds, "GlobalIds needs to be provided for " << embeddedSurfaceBlockName);
+
+  // build and set global to local mapping to use later for 'fracture_to_parent_matrix_cell_mapping'
+  embeddedSurfBlock.setGlobalToLocalMap( buildGlobalToLocal( edfmMeshCellGlobalIds ) );
 
-  // Ouassim change
-  vtkIdType numEdfmFracs = embeddedSurfaceMesh->GetNumberOfCells();
   //1.Get EDFM vertices
-  vtkUnstructuredGrid * grid = vtkUnstructuredGrid::SafeDownCast( embeddedSurfaceMesh );
+  vtkUnstructuredGrid * const grid = vtkUnstructuredGrid::SafeDownCast( embeddedSurfaceMesh );
   vtkPoints * const nodes = grid->GetPoints();
-  vtkIdType numNodes = nodes ? nodes->GetNumberOfPoints() : 0;
+  vtkIdType const numNodes = nodes ? nodes->GetNumberOfPoints() : 0;
 
   ArrayOfArrays< localIndex > elem2dToNodes( numEdfmFracs );
   for( vtkIdType i = 0; i < numEdfmFracs; ++i )
   {
-    auto val = grid->GetCell( i )->GetPointIds();
+    auto const val = grid->GetCell( i )->GetPointIds();
     elem2dToNodes.emplaceBack( i, val->GetId( 0 ));
     elem2dToNodes.emplaceBack( i, val->GetId( 1 ));
     elem2dToNodes.emplaceBack( i, val->GetId( 2 ));
     elem2dToNodes.emplaceBack( i, val->GetId( 3 ));
   }
 
-
   ArrayOfArrays< real64 > embeddedSurfaceElemNodes( LvArray::integerConversion< localIndex >( numNodes ) );
   for( vtkIdType i = 0; i < numNodes; ++i )
   {
@@ -59,10 +80,10 @@ void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName,
   //2. Get EDFM normal vectors
   ArrayOfArrays< real64 > embeddedSurfaceNormalVectors( numEdfmFracs );
   const char * const normal_str = "normal_vectors";
-  vtkDataArray * normals =grid->GetCellData()->GetArray( normal_str );
+  vtkDataArray * const normals = grid->GetCellData()->GetArray( normal_str );
   for( vtkIdType i = 0; i < numEdfmFracs; ++i )
   {
-    auto val = normals->GetTuple3( i );
+    auto const val = normals->GetTuple3( i );
     embeddedSurfaceNormalVectors.emplaceBack( i, val[0] );
     embeddedSurfaceNormalVectors.emplaceBack( i, val[1] );
     embeddedSurfaceNormalVectors.emplaceBack( i, val[2] );
@@ -71,10 +92,10 @@ void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName,
   //3. Get EDFM  tangential length vectors (horizontal)
   ArrayOfArrays< real64 > embeddedSurfaceTangentialLengthVectors( numEdfmFracs );
   const char * const tl_str = "tangential_length_vectors";
-  vtkDataArray * tangential_len =grid->GetCellData()->GetArray( tl_str );
+  vtkDataArray * const tangential_len = grid->GetCellData()->GetArray( tl_str );
   for( vtkIdType i = 0; i < numEdfmFracs; ++i )
   {
-    auto val = tangential_len->GetTuple3( i );
+    auto const val = tangential_len->GetTuple3( i );
     embeddedSurfaceTangentialLengthVectors.emplaceBack( i, val[0] );
     embeddedSurfaceTangentialLengthVectors.emplaceBack( i, val[1] );
     embeddedSurfaceTangentialLengthVectors.emplaceBack( i, val[2] );
@@ -83,10 +104,10 @@ void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName,
   //4. Get EDFM  tangential width vectors (vertical)
   ArrayOfArrays< real64 > embeddedSurfaceTangentialWidthVectors( numEdfmFracs );
   const char * const tw_str = "tangential_width_vectors";
-  vtkDataArray * tangential_width =grid->GetCellData()->GetArray( tw_str );
+  vtkDataArray * const tangential_width = grid->GetCellData()->GetArray( tw_str );
   for( vtkIdType i = 0; i < numEdfmFracs; ++i )
   {
-    auto val = tangential_width->GetTuple3( i );
+    auto const val = tangential_width->GetTuple3( i );
     embeddedSurfaceTangentialWidthVectors.emplaceBack( i, val[0] );
     embeddedSurfaceTangentialWidthVectors.emplaceBack( i, val[1] );
     embeddedSurfaceTangentialWidthVectors.emplaceBack( i, val[2] );
@@ -95,17 +116,17 @@ void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName,
   //4. Get EDFM aperture
   array1d< real64 > embeddedSurfaceAperture( numEdfmFracs );
   const char * const aperture_str = "aperture";
-  vtkDataArray * apertures =grid->GetCellData()->GetArray( aperture_str );
+  vtkDataArray * const apertures = grid->GetCellData()->GetArray( aperture_str );
   for( vtkIdType i = 0; i < numEdfmFracs; ++i )
   {
-    auto val =apertures->GetTuple1( i );
+    auto const val = apertures->GetTuple1( i );
     embeddedSurfaceAperture[i] = val;
   }
 
   //5. Get EDFM permeability
   array1d< real64 > embeddedSurfacePermeability( numEdfmFracs );
   const char * const perm_str = "permeability";
-  vtkDataArray * perms =grid->GetCellData()->GetArray( perm_str );
+  vtkDataArray * const perms = grid->GetCellData()->GetArray( perm_str );
   for( vtkIdType i = 0; i < numEdfmFracs; ++i )
   {
     auto val =perms->GetTuple1( i );
@@ -114,20 +135,26 @@ void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName,
 
   //6. Get edfm to matrix cell mapping
   const char * const fid_to_mid_mapping  = "fracture_to_parent_matrix_cell_mapping";
-  vtkDataArray * fid_mid =grid->GetCellData()->GetArray( fid_to_mid_mapping );
+  vtkDataArray * const fid_mid = grid->GetCellData()->GetArray( fid_to_mid_mapping );
+  localIndex const blockIndex = 0; // TODO what happens when there are multiple blocks?
+  CellBlock const & cellBlock = cellBlockManager.getCellBlock( blockIndex );
+  auto const & globalToLocalCell = cellBlock.globalToLocalMap();
+  auto const & globalToLocalEdfm = embeddedSurfBlock.globalToLocalMap();
 
   ArrayOfArrays< localIndex > toBlockIndex( numEdfmFracs );
   ArrayOfArrays< localIndex > toCellIndex( numEdfmFracs );
   for( vtkIdType i = 0; i < numEdfmFracs; ++i )
   {
-    auto fidmid = fid_mid->GetTuple2( i ); // TODO these indexes are global and needs to be converted to local
-    toBlockIndex.emplaceBack( fidmid[0], 0 );// cell block is set to 0 for now
-    toCellIndex.emplaceBack( fidmid[0], fidmid[1] );
+    // these indexes are global and needs to be converted to local
+    auto const fidmid = fid_mid->GetTuple2( i );
+    // convert global to local using the maps
+    localIndex const id0 = globalToLocalEdfm.at( fidmid[0] );
+    localIndex const id1 = globalToLocalCell.at( fidmid[1] );
+    toBlockIndex.emplaceBack( id0, blockIndex );// cell block is set to 0 for now
+    toCellIndex.emplaceBack( id0, id1 );
   }
   ToCellRelation< ArrayOfArrays< localIndex > > elem2dTo3d( std::move( toBlockIndex ), std::move( toCellIndex ) );
 
-  EmbeddedSurfaceBlock & embeddedSurfBlock = cellBlockManager.registerEmbeddedSurfaceBlock( embeddedSurfaceBlockName );
-
   embeddedSurfBlock.setNumEmbeddedSurfElem( numEdfmFracs );
   embeddedSurfBlock.setEmbeddedSurfElemNodes( std::move( embeddedSurfaceElemNodes ));
   embeddedSurfBlock.setEmbeddedSurfElemNormalVectors( std::move( embeddedSurfaceNormalVectors ));

src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp

@@ -29,12 +29,10 @@ namespace geos::vtk
  * @brief Attach the embedded surface block information to the cell block manager.
  * @param embeddedSurfaceBlockName[in] The name of the embedded surface block.
  * @param embeddedSurfaceMesh[in] The vtk mesh for the embedded surface block.
- * @param mesh[in] The vtk volumic mesh.
  * @param cellBlockManager[inout] The cell block manager that will receive the embedded surface block information.
  */
 void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName,
                                     vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh,
-                                    vtkSmartPointer< vtkDataSet > mesh,
                                     CellBlockManager & cellBlockManager );
 }
 

src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp

@@ -605,7 +605,6 @@ array1d< globalIndex > buildLocalToGlobal( vtkIdTypeArray const * faceMeshCellGl
   return l2g;
 }
 
-
 void importFractureNetwork( string const & faceBlockName,
                             vtkSmartPointer< vtkDataSet > faceMesh,
                             vtkSmartPointer< vtkDataSet > mesh,

src/coreComponents/mesh/generators/VTKMeshGenerator.cpp

@@ -228,7 +228,7 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager
   for( auto const & [name, mesh]: m_embeddedSurfaceBlockMeshes )
   {
     GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{} '{}': importing embedded fracture network {}...", catalogName(), getName(), name ) );
-    vtk::importEmbeddedFractureNetwork( name, mesh, m_vtkMesh, cellBlockManager );
+    vtk::importEmbeddedFractureNetwork( name, mesh, cellBlockManager );
   }
 
   GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{} '{}': done!", catalogName(), getName() ) );

src/coreComponents/mesh/generators/VTKUtilities.cpp

@@ -754,13 +754,15 @@ vtkSmartPointer< vtkDataSet > manageGlobalIds( vtkSmartPointer< vtkDataSet > mes
   {
     // Add global ids on the fly if needed
     int const me = hasGlobalIds( mesh );
+    GEOS_LOG_RANK( "me hasGlobalIds="<<me );
     int everyone;
     MpiWrapper::allReduce( &me, &everyone, 1, MPI_MAX, MPI_COMM_GEOS );
 
     if( everyone and not me )
     {
       mesh->GetPointData()->SetGlobalIds( vtkIdTypeArray::New() );
       mesh->GetCellData()->SetGlobalIds( vtkIdTypeArray::New() );
+      GEOS_LOG_RANK( "SetGlobalIds" );
     }
   }
 
@@ -930,6 +932,7 @@ redistributeMeshes( integer const logLevel,
 
   // Generate global IDs for vertices and cells, if needed
   vtkSmartPointer< vtkDataSet > mesh = manageGlobalIds( loadedMesh, useGlobalIds, !std::empty( fractures ) );
+  //vtkSmartPointer< vtkDataSet > mesh = manageGlobalIds( loadedMesh, useGlobalIds, !std::empty( fractures ) || !std::empty( edfms ) );
 
   if( MpiWrapper::commRank( comm ) != ( MpiWrapper::commSize( comm ) - 1 ) )
   {
@@ -1795,6 +1798,7 @@ void fillCellBlock( vtkDataSet & mesh,
   localIndex const numNodesPerElement = cellBlock.numNodesPerElement();
   arrayView2d< localIndex, cells::NODE_MAP_USD > const cellToVertex = cellBlock.getElemToNode();
   arrayView1d< globalIndex > const & localToGlobal = cellBlock.localToGlobalMap();
+  auto & globalToLocal = cellBlock.globalToLocalMap();
   vtkIdTypeArray const * const globalCellId = vtkIdTypeArray::FastDownCast( mesh.GetCellData()->GetGlobalIds() );
   GEOS_ERROR_IF( !cellIds.empty() && globalCellId == nullptr, "Global cell IDs have not been generated" );
 
@@ -1805,7 +1809,10 @@ void fillCellBlock( vtkDataSet & mesh,
     {
       cellToVertex[cellCount][v] = cell->GetPointId( nodeOrder[v] );
     }
-    localToGlobal[cellCount++] = globalCellId->GetValue( c );
+    globalIndex g = globalCellId->GetValue( c );
+    localToGlobal[cellCount] = g;
+    globalToLocal[g] = cellCount;
+    cellCount++;
   };
 
   // Writing connectivity and Local to Global