diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp
index 7109599cac4..36aa651d318 100644
--- a/src/coreComponents/mainInterface/ProblemManager.cpp
+++ b/src/coreComponents/mainInterface/ProblemManager.cpp
@@ -576,12 +576,9 @@ void ProblemManager::generateMesh()
 
       elemManager.generateMesh( cellBlockManager );
 
-      nodeManager.setGeometricalRelations( cellBlockManager );
+      nodeManager.setGeometricalRelations( cellBlockManager, elemManager );
       edgeManager.setGeometricalRelations( cellBlockManager );
-      faceManager.setGeometricalRelations( cellBlockManager, nodeManager );
-
-      nodeManager.buildRegionMaps( elemManager );
-      faceManager.buildRegionMaps( elemManager );
+      faceManager.setGeometricalRelations( cellBlockManager, elemManager, nodeManager );
 
       nodeManager.constructGlobalToLocalMap( cellBlockManager );
 
diff --git a/src/coreComponents/mesh/CellElementRegion.hpp b/src/coreComponents/mesh/CellElementRegion.hpp
index 8791857985c..ad10ac97d31 100644
--- a/src/coreComponents/mesh/CellElementRegion.hpp
+++ b/src/coreComponents/mesh/CellElementRegion.hpp
@@ -66,14 +66,14 @@ class CellElementRegion : public ElementRegionBase
    * @brief The key name for the FaceElementRegion in the object catalog.
    * @return A string containing the key name.
    */
-  static const string catalogName()
+  static string catalogName()
   { return "CellElementRegion"; }
 
   /**
    * @copydoc catalogName()
    */
-  virtual const string getCatalogName() const override final
-  { return CellElementRegion::catalogName(); }
+  virtual string getCatalogName() const override final
+  { return catalogName(); }
 
   ///@}
 
diff --git a/src/coreComponents/mesh/CellElementSubRegion.cpp b/src/coreComponents/mesh/CellElementSubRegion.cpp
index 264a1b9c490..7c84da5e6ea 100644
--- a/src/coreComponents/mesh/CellElementSubRegion.cpp
+++ b/src/coreComponents/mesh/CellElementSubRegion.cpp
@@ -330,7 +330,9 @@ void CellElementSubRegion::
 void CellElementSubRegion::calculateElementGeometricQuantities( NodeManager const & nodeManager,
                                                                 FaceManager const & GEOSX_UNUSED_PARAM( faceManager ) )
 {
-  arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X = nodeManager.referencePosition();
+  GEOSX_MARK_FUNCTION;
+
+  arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.referencePosition();
 
   forAll< parallelHostPolicy >( this->size(), [=] ( localIndex const k )
   {
diff --git a/src/coreComponents/mesh/CellElementSubRegion.hpp b/src/coreComponents/mesh/CellElementSubRegion.hpp
index 5cab955264a..5bee78b2328 100644
--- a/src/coreComponents/mesh/CellElementSubRegion.hpp
+++ b/src/coreComponents/mesh/CellElementSubRegion.hpp
@@ -50,14 +50,14 @@ class CellElementSubRegion : public ElementSubRegionBase
    * @brief Const getter for the catalog name.
    * @return the name of this type in the catalog
    */
-  static const string catalogName()
+  static string catalogName()
   { return "CellElementSubRegion"; }
 
   /**
    * @copydoc catalogName()
    */
-  virtual const string getCatalogName() const override final
-  { return CellElementSubRegion::catalogName(); }
+  virtual string getCatalogName() const override final
+  { return catalogName(); }
 
   /**
    * @name Constructor / Destructor
diff --git a/src/coreComponents/mesh/EdgeManager.cpp b/src/coreComponents/mesh/EdgeManager.cpp
index e21938cb622..1ccc462e762 100644
--- a/src/coreComponents/mesh/EdgeManager.cpp
+++ b/src/coreComponents/mesh/EdgeManager.cpp
@@ -94,21 +94,27 @@ void EdgeManager::buildEdges( localIndex const numNodes,
                               ArrayOfArraysView< localIndex const > const & faceToNodeMap,
                               ArrayOfArrays< localIndex > & faceToEdgeMap )
 {
+  ArrayOfArrays< localIndex > edgeToFace;
   localIndex const numEdges = buildEdgeMaps( numNodes,
                                              faceToNodeMap,
                                              faceToEdgeMap,
-                                             m_toFacesRelation,
+                                             edgeToFace,
                                              m_toNodesRelation );
 
+  m_toFacesRelation.base().assimilate< parallelHostPolicy >( std::move( edgeToFace ),
+                                                             LvArray::sortedArrayManipulation::UNSORTED_NO_DUPLICATES );
   resize( numEdges );
 }
 
 void EdgeManager::setGeometricalRelations( CellBlockManagerABC const & cellBlockManager )
 {
+  GEOSX_MARK_FUNCTION;
+
   resize( cellBlockManager.numEdges() );
 
   m_toNodesRelation.base() = cellBlockManager.getEdgeToNodes();
-  m_toFacesRelation.base() = cellBlockManager.getEdgeToFaces();
+  m_toFacesRelation.base().assimilate< parallelHostPolicy >( cellBlockManager.getEdgeToFaces(),
+                                                             LvArray::sortedArrayManipulation::UNSORTED_NO_DUPLICATES );
 }
 
 void EdgeManager::setupRelatedObjectsInRelations( NodeManager const & nodeManager,
@@ -120,31 +126,25 @@ void EdgeManager::setupRelatedObjectsInRelations( NodeManager const & nodeManage
 
 void EdgeManager::setDomainBoundaryObjects( FaceManager const & faceManager )
 {
-  // get the "isDomainBoundary" field from the faceManager. This should have
-  // been set already!
-  arrayView1d< integer const > const & isFaceOnDomainBoundary = faceManager.getDomainBoundaryIndicator();
+  // get the "isDomainBoundary" field from the faceManager. This should have been set already!
+  arrayView1d< integer const > const isFaceOnDomainBoundary = faceManager.getDomainBoundaryIndicator();
 
   // get the "isDomainBoundary" field from for *this, and set it to zero
-  arrayView1d< integer > const & isEdgeOnDomainBoundary = this->getDomainBoundaryIndicator();
+  arrayView1d< integer > const isEdgeOnDomainBoundary = this->getDomainBoundaryIndicator();
   isEdgeOnDomainBoundary.zero();
 
-  ArrayOfArraysView< localIndex const > const & faceToEdgeMap = faceManager.edgeList().toViewConst();
+  ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst();
 
-  // loop through all faces
-  for( localIndex kf = 0; kf < faceManager.size(); ++kf )
+  forAll< parallelHostPolicy >( faceManager.size(), [=]( localIndex const faceIndex )
   {
-    // check to see if the face is on a domain boundary
-    if( isFaceOnDomainBoundary[kf] == 1 )
+    if( isFaceOnDomainBoundary[faceIndex] == 1 )
     {
-      localIndex const numFaceEdges = faceToEdgeMap.sizeOfArray( kf );
-
-      // loop over all nodes connected to face, and set isNodeDomainBoundary
-      for( localIndex a = 0; a < numFaceEdges; ++a )
+      for( localIndex const edgeIndex : faceToEdgeMap[faceIndex] )
       {
-        isEdgeOnDomainBoundary( faceToEdgeMap( kf, a ) ) = 1;
+        isEdgeOnDomainBoundary[edgeIndex] = 1;
       }
     }
-  }
+  } );
 }
 
 bool EdgeManager::hasNode( const localIndex edgeIndex, const localIndex nodeIndex ) const
diff --git a/src/coreComponents/mesh/EdgeManager.hpp b/src/coreComponents/mesh/EdgeManager.hpp
index 70864f105c3..fb21d159d96 100644
--- a/src/coreComponents/mesh/EdgeManager.hpp
+++ b/src/coreComponents/mesh/EdgeManager.hpp
@@ -56,7 +56,7 @@ class EdgeManager : public ObjectManagerBase
   /**
    * @return the string representing the edge manager name in the catalog
    */
-  static const string catalogName()
+  static string catalogName()
   { return "EdgeManager"; }
 
 
@@ -64,8 +64,8 @@ class EdgeManager : public ObjectManagerBase
    * @brief Getter used to access the edge manager catalog name.
    * @return the edge manager catalog name
    */
-  virtual const string getCatalogName() const override
-  { return EdgeManager::catalogName(); }
+  virtual string getCatalogName() const override
+  { return catalogName(); }
 
   ///@}
 
@@ -103,9 +103,9 @@ class EdgeManager : public ObjectManagerBase
 
   /**
    * @brief Set the node of the domain boundary object.
-   * @param[in] faceManager The reference of the face manager.
+   * @param[in] faceIndex The reference of the face manager.
    */
-  void setDomainBoundaryObjects( FaceManager const & faceManager );
+  void setDomainBoundaryObjects( FaceManager const & faceIndex );
 
   /**
    * @brief Set external edges.
diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp
index 5c42257ba3f..6fe89fb1cac 100644
--- a/src/coreComponents/mesh/ElementRegionManager.cpp
+++ b/src/coreComponents/mesh/ElementRegionManager.cpp
@@ -59,11 +59,7 @@ void ElementRegionManager::setMaxGlobalIndex()
     m_localMaxGlobalIndex = std::max( m_localMaxGlobalIndex, subRegion.maxGlobalIndex() );
   } );
 
-  MpiWrapper::allReduce( &m_localMaxGlobalIndex,
-                         &m_maxGlobalIndex,
-                         1,
-                         MPI_MAX,
-                         MPI_COMM_GEOSX );
+  m_maxGlobalIndex = MpiWrapper::max( m_localMaxGlobalIndex, MPI_COMM_GEOSX );
 }
 
 
@@ -630,6 +626,33 @@ ElementRegionManager::unpackFracturedElements( buffer_unit_type const * & buffer
   return unpackedSize;
 }
 
+array2d< localIndex >
+ElementRegionManager::getCellBlockToSubRegionMap( CellBlockManagerABC const & cellBlockManager ) const
+{
+  array2d< localIndex > blockMap( cellBlockManager.getCellBlocks().numSubGroups(), 2 );
+  blockMap.setValues< serialPolicy >( -1 );
+
+  Group::subGroupMap const & cellBlocks = cellBlockManager.getCellBlocks().getSubGroups();
+
+  forElementSubRegionsComplete< CellElementSubRegion >( [blockMap = blockMap.toView(),
+                                                         &cellBlocks]( localIndex const er,
+                                                                       localIndex const esr,
+                                                                       ElementRegionBase const & region,
+                                                                       CellElementSubRegion const & subRegion )
+  {
+    localIndex const blockIndex = cellBlocks.getIndex( subRegion.getName() );
+    GEOSX_ERROR_IF( blockIndex == Group::subGroupMap::KeyIndex::invalid_index,
+                    GEOSX_FMT( "Cell block not found for subregion {}/{}", region.getName(), subRegion.getName() ) );
+    GEOSX_ERROR_IF( blockMap( blockIndex, 1 ) != -1,
+                    GEOSX_FMT( "Cell block {} mapped to more than one subregion", subRegion.getName() ) );
+
+    blockMap( blockIndex, 0 ) = er;
+    blockMap( blockIndex, 1 ) = esr;
+  } );
+
+  return blockMap;
+}
+
 
 REGISTER_CATALOG_ENTRY( ObjectManagerBase, ElementRegionManager, string const &, Group * const )
 }
diff --git a/src/coreComponents/mesh/ElementRegionManager.hpp b/src/coreComponents/mesh/ElementRegionManager.hpp
index a466f308be7..e4faad83f53 100644
--- a/src/coreComponents/mesh/ElementRegionManager.hpp
+++ b/src/coreComponents/mesh/ElementRegionManager.hpp
@@ -105,15 +105,15 @@ class ElementRegionManager : public ObjectManagerBase
    * @brief The function is to return the name of the ElementRegionManager in the object catalog
    * @return string that contains the catalog name used to register/lookup this class in  the object catalog
    */
-  static const string catalogName()
+  static string catalogName()
   { return "ZoneManager"; }
 
   /**
    * @brief Virtual access to catalogName()
    * @return string that contains the catalog name used to register/lookup this class in the object catalog
    */
-  virtual const string getCatalogName() const override final
-  { return ElementRegionManager::catalogName(); }
+  virtual string getCatalogName() const override final
+  { return catalogName(); }
 
   /**
    * @brief Constructor.
@@ -264,6 +264,14 @@ class ElementRegionManager : public ObjectManagerBase
     return this->getRegions().size();
   }
 
+  /**
+   * @brief Produce a map from cell block indices to element region and subregion indices
+   * @param cellBlockManager the CellBlocKManager
+   * @return a (numBlock x 2) array with each row corresponding to a cell block and containing
+   *         region (first entry) and subregion (second entry) indices, or -1 if block was not used.
+   */
+  array2d< localIndex > getCellBlockToSubRegionMap( CellBlockManagerABC const & cellBlockManager ) const;
+
   /**
    * @brief This function is used to launch kernel function over all the element regions with region type =
    * ElementRegionBase.
diff --git a/src/coreComponents/mesh/EmbeddedSurfaceNodeManager.hpp b/src/coreComponents/mesh/EmbeddedSurfaceNodeManager.hpp
index 73fcd34bbf3..877a0b8d3b3 100644
--- a/src/coreComponents/mesh/EmbeddedSurfaceNodeManager.hpp
+++ b/src/coreComponents/mesh/EmbeddedSurfaceNodeManager.hpp
@@ -102,8 +102,8 @@ class EmbeddedSurfaceNodeManager : public ObjectManagerBase
    * @brief Provide a virtual access to catalogName().
    * @return string that contains the EmbeddedSurfaceNodeManager catalog name
    */
-  const string getCatalogName() const override final
-  { return EmbeddedSurfaceNodeManager::catalogName(); }
+  string getCatalogName() const override final
+  { return catalogName(); }
 
   ///@}
 
diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp
index 126dc8c27f0..3e3c31e6b60 100644
--- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp
+++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp
@@ -79,16 +79,16 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion
    * @brief Get catalog name.
    * @return the catalog name
    */
-  static const string catalogName()
+  static string catalogName()
   { return "EmbeddedSurfaceSubRegion"; }
 
   /**
    * @brief Get catalog name.
    * @return the catalog name
    */
-  virtual const string getCatalogName() const override
+  virtual string getCatalogName() const override
   {
-    return EmbeddedSurfaceSubRegion::catalogName();
+    return catalogName();
   }
 
   ///@}
diff --git a/src/coreComponents/mesh/FaceElementSubRegion.hpp b/src/coreComponents/mesh/FaceElementSubRegion.hpp
index 5f53de07440..b5884ba6890 100644
--- a/src/coreComponents/mesh/FaceElementSubRegion.hpp
+++ b/src/coreComponents/mesh/FaceElementSubRegion.hpp
@@ -48,16 +48,16 @@ class FaceElementSubRegion : public SurfaceElementSubRegion
    * @brief Get catalog name.
    * @return the catalog name
    */
-  static const string catalogName()
+  static string catalogName()
   { return "FaceElementSubRegion"; }
 
   /**
    * @brief Get catalog name.
    * @return the catalog name
    */
-  virtual const string getCatalogName() const override
+  virtual string getCatalogName() const override
   {
-    return FaceElementSubRegion::catalogName();
+    return catalogName();
   }
 
   ///@}
diff --git a/src/coreComponents/mesh/FaceManager.cpp b/src/coreComponents/mesh/FaceManager.cpp
index 66eb142d9d7..dfc9d960c76 100644
--- a/src/coreComponents/mesh/FaceManager.cpp
+++ b/src/coreComponents/mesh/FaceManager.cpp
@@ -36,8 +36,8 @@ using namespace dataRepository;
 FaceManager::FaceManager( string const &, Group * const parent ):
   ObjectManagerBase( "FaceManager", parent )
 {
-  this->registerWrapper( viewKeyStruct::nodeListString(), &m_nodeList );
-  this->registerWrapper( viewKeyStruct::edgeListString(), &m_edgeList );
+  this->registerWrapper( viewKeyStruct::nodeListString(), &m_toNodesRelation );
+  this->registerWrapper( viewKeyStruct::edgeListString(), &m_toEdgesRelation );
 
   this->registerWrapper( viewKeyStruct::elementRegionListString(), &m_toElements.m_toElementRegion ).
     setApplyDefaultValue( -1 );
@@ -63,34 +63,20 @@ FaceManager::FaceManager( string const &, Group * const parent ):
 
 void FaceManager::resize( localIndex const newSize )
 {
-  m_nodeList.resize( newSize, 2 * nodeMapOverallocation() );
-  m_edgeList.resize( newSize, 2 * edgeMapOverallocation() );
+  m_toNodesRelation.resize( newSize, 2 * nodeMapOverallocation() );
+  m_toEdgesRelation.resize( newSize, 2 * edgeMapOverallocation() );
   ObjectManagerBase::resize( newSize );
 }
 
-void FaceManager::buildRegionMaps( ElementRegionManager const & elementRegionManager )
+void FaceManager::buildSets( NodeManager const & nodeManager )
 {
   GEOSX_MARK_FUNCTION;
 
-  auto const faceListGetter =
-    []( CellElementSubRegion const & subregion ) -> CellElementSubRegion::FaceMapType const &
-  { return subregion.faceList(); };
-
-  meshMapUtilities::buildToElementMaps( elementRegionManager,
-                                        size(),
-                                        m_toElements,
-                                        faceListGetter );
-}
-
-void FaceManager::buildSets( NodeManager const & nodeManager )
-{
   // First create the sets
   auto const & nodeSets = nodeManager.sets().wrappers();
-  for( localIndex i = 0; i < nodeSets.size(); ++i )
+  for( auto const & setWrapper : nodeSets )
   {
-    auto const & setWrapper = nodeSets[i];
-    string const & setName = setWrapper->getName();
-    createSet( setName );
+    createSet( setWrapper.second->getName() );
   }
 
   // Then loop over them in parallel and fill them in.
@@ -98,19 +84,22 @@ void FaceManager::buildSets( NodeManager const & nodeManager )
   {
     auto const & setWrapper = nodeSets[i];
     string const & setName = setWrapper->getName();
-    SortedArrayView< localIndex const > const & targetSet = nodeManager.sets().getReference< SortedArray< localIndex > >( setName ).toViewConst();
-    constructSetFromSetAndMap( targetSet, m_nodeList.toViewConst(), setName );
+    SortedArrayView< localIndex const > const targetSet =
+      nodeManager.sets().getReference< SortedArray< localIndex > >( setName ).toViewConst();
+    constructSetFromSetAndMap( targetSet, m_toNodesRelation.toViewConst(), setName );
   } );
 }
 
 void FaceManager::setDomainBoundaryObjects()
 {
-  arrayView1d< integer > const & isFaceOnDomainBoundary = getDomainBoundaryIndicator();
+  arrayView1d< integer > const isFaceOnDomainBoundary = getDomainBoundaryIndicator();
   isFaceOnDomainBoundary.zero();
 
-  forAll< parallelHostPolicy >( size(), [&]( localIndex const kf )
+  arrayView2d< localIndex const > const toElementRegion = m_toElements.m_toElementRegion.toViewConst();
+
+  forAll< parallelHostPolicy >( size(), [=]( localIndex const kf )
   {
-    if( m_toElements.m_toElementRegion[kf][1] == -1 )
+    if( toElementRegion( kf, 1 ) == -1 )
     {
       isFaceOnDomainBoundary( kf ) = 1;
     }
@@ -118,12 +107,21 @@ void FaceManager::setDomainBoundaryObjects()
 }
 
 void FaceManager::setGeometricalRelations( CellBlockManagerABC const & cellBlockManager,
+                                           ElementRegionManager const & elemRegionManager,
                                            NodeManager const & nodeManager )
 {
+  GEOSX_MARK_FUNCTION;
+
   resize( cellBlockManager.numFaces() );
 
-  m_nodeList.base() = cellBlockManager.getFaceToNodes();
-  m_edgeList.base() = cellBlockManager.getFaceToEdges();
+  m_toNodesRelation.base() = cellBlockManager.getFaceToNodes();
+  m_toEdgesRelation.base() = cellBlockManager.getFaceToEdges();
+
+  ToCellRelation< array2d< localIndex > > const toCellBlock = cellBlockManager.getFaceToElements();
+  array2d< localIndex > const blockToSubRegion = elemRegionManager.getCellBlockToSubRegionMap( cellBlockManager );
+  meshMapUtilities::transformCellBlockToRegionMap< parallelHostPolicy >( blockToSubRegion.toViewConst(),
+                                                                         toCellBlock,
+                                                                         m_toElements );
 
   computeGeometry( nodeManager );
 }
@@ -132,8 +130,8 @@ void FaceManager::setupRelatedObjectsInRelations( NodeManager const & nodeManage
                                                   EdgeManager const & edgeManager,
                                                   ElementRegionManager const & elementRegionManager )
 {
-  m_nodeList.setRelatedObject( nodeManager );
-  m_edgeList.setRelatedObject( edgeManager );
+  m_toNodesRelation.setRelatedObject( nodeManager );
+  m_toEdgesRelation.setRelatedObject( edgeManager );
 
   m_toElements.setElementRegionManager( elementRegionManager );
 }
@@ -145,7 +143,7 @@ void FaceManager::computeGeometry( NodeManager const & nodeManager )
   // loop over faces and calculate faceArea, faceNormal and faceCenter
   forAll< parallelHostPolicy >( this->size(), [&]( localIndex const faceIndex )
   {
-    m_faceArea[ faceIndex ] = computationalGeometry::centroid_3DPolygon( m_nodeList[ faceIndex ],
+    m_faceArea[ faceIndex ] = computationalGeometry::centroid_3DPolygon( m_toNodesRelation[ faceIndex ],
                                                                          X,
                                                                          m_faceCenter[ faceIndex ],
                                                                          m_faceNormal[ faceIndex ] );
@@ -167,18 +165,6 @@ void FaceManager::setIsExternal()
   }
 }
 
-localIndex FaceManager::getMaxFaceNodes() const
-{
-  localIndex maxSize = 0;
-  ArrayOfArraysView< localIndex const > const & faceToNodeMap = nodeList().toViewConst();
-  for( localIndex kf =0; kf < size(); ++kf )
-  {
-    maxSize = std::max( maxSize, faceToNodeMap.sizeOfArray( kf ) );
-  }
-
-  return maxSize;
-}
-
 void FaceManager::sortAllFaceNodes( NodeManager const & nodeManager,
                                     ElementRegionManager const & elemManager )
 {
@@ -187,55 +173,54 @@ void FaceManager::sortAllFaceNodes( NodeManager const & nodeManager,
   arrayView2d< localIndex const > const facesToElementRegions = elementRegionList();
   arrayView2d< localIndex const > const facesToElementSubRegions = elementSubRegionList();
   arrayView2d< localIndex const > const facesToElements = elementList();
-  arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X = nodeManager.referencePosition();
-
-  ArrayOfArraysView< localIndex > const & facesToNodes = nodeList().toView();
+  arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.referencePosition();
 
-  GEOSX_ERROR_IF( getMaxFaceNodes() >= MAX_FACE_NODES, "More nodes on a face than expected!" );
+  ArrayOfArraysView< localIndex > const facesToNodes = nodeList().toView();
 
   elemManager.forElementSubRegions< CellElementSubRegion >( [&] ( CellElementSubRegion const & subRegion )
-  { subRegion.calculateElementCenters( X ); } );
+  {
+    subRegion.calculateElementCenters( X );
+  } );
 
-  forAll< parallelHostPolicy >( size(), [&]( localIndex const faceIndex ) -> void
+  ElementRegionManager::ElementViewAccessor< arrayView2d< real64 const > > elemCenter =
+    elemManager.constructArrayViewAccessor< real64, 2 >( ElementSubRegionBase::viewKeyStruct::elementCenterString() );
+
+  forAll< parallelHostPolicy >( size(), [=, elemCenter = elemCenter.toNestedViewConst()]( localIndex const faceIndex )
   {
     // The face should be connected to at least one element.
-    if( facesToElements( faceIndex, 0 ) != -1 or facesToElements( faceIndex, 1 ) != -1 )
+    if( facesToElements( faceIndex, 0 ) < 0 && facesToElements( faceIndex, 1 ) < 0 )
     {
-      // Take the first defined face-to-(elt/region/sub region) to sorting direction.
-      localIndex const iElemLoc = facesToElements( faceIndex, 0 ) > -1 ? 0 : 1;
-
-      localIndex const er = facesToElementRegions[faceIndex][iElemLoc];
-      localIndex const esr = facesToElementSubRegions[faceIndex][iElemLoc];
-      localIndex const ei = facesToElements( faceIndex, iElemLoc );
-      if( er != -1 and esr != -1 and ei != -1 )
-      {
-        CellElementSubRegion const & subRegion = elemManager.getRegion( er ).getSubRegion< CellElementSubRegion >( esr );
-        arrayView2d< real64 const > const elemCenter = subRegion.getElementCenter();
-        localIndex const numFaceNodes = facesToNodes.sizeOfArray( faceIndex );
-        sortFaceNodes( X, elemCenter[ei], facesToNodes[faceIndex], numFaceNodes );
-      }
-      else
-      {
-        GEOSX_ERROR( "Face " << faceIndex << " is connected to at least one invalid region (" << er << "), sub-region (" << esr << ") or element (" << ei << ")." );
-      }
+      GEOSX_ERROR( "Face " << faceIndex << " is not connected to an element." );
     }
-    else
+
+    // Take the first defined face-to-(elt/region/sub region) to sorting direction.
+    localIndex const iElemLoc = facesToElements( faceIndex, 0 ) >= 0 ? 0 : 1;
+
+    localIndex const er = facesToElementRegions( faceIndex, iElemLoc );
+    localIndex const esr = facesToElementSubRegions( faceIndex, iElemLoc );
+    localIndex const ei = facesToElements( faceIndex, iElemLoc );
+
+    if( er < 0 || esr < 0 || ei < 0 )
     {
-      GEOSX_ERROR( "Face " << faceIndex << " does not seem connected to any element." );
+      GEOSX_ERROR( GEOSX_FMT( "Face {} is connected to an invalid element ({}/{}/{}).", faceIndex, er, esr, ei ) );
     }
+
+    sortFaceNodes( X, elemCenter[er][esr][ei], facesToNodes[faceIndex] );
   } );
 }
 
 void FaceManager::sortFaceNodes( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X,
                                  arraySlice1d< real64 const > const elementCenter,
-                                 localIndex * const faceNodes,
-                                 localIndex const numFaceNodes )
+                                 Span< localIndex > const faceNodes )
 {
+  localIndex const numFaceNodes = LvArray::integerConversion< localIndex >( faceNodes.size() );
+  GEOSX_ERROR_IF_GT_MSG( numFaceNodes, MAX_FACE_NODES, "Node per face limit exceeded" );
+
   localIndex const firstNodeIndex = faceNodes[0];
 
   // get face center (average vertex location)
   real64 fc[3] = { 0 };
-  for( localIndex n =0; n < numFaceNodes; ++n )
+  for( localIndex n = 0; n < numFaceNodes; ++n )
   {
     LvArray::tensorOps::add< 3 >( fc, X[faceNodes[n]] );
   }
@@ -280,7 +265,7 @@ void FaceManager::sortFaceNodes( arrayView2d< real64 const, nodes::REFERENCE_POS
     std::pair< real64, localIndex > thetaOrder[MAX_FACE_NODES];
 
     // Sort nodes counterclockwise around face center
-    for( localIndex n =0; n < numFaceNodes; ++n )
+    for( localIndex n = 0; n < numFaceNodes; ++n )
     {
       real64 v[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[faceNodes[n]] );
       LvArray::tensorOps::subtract< 3 >( v, fc );
@@ -290,7 +275,7 @@ void FaceManager::sortFaceNodes( arrayView2d< real64 const, nodes::REFERENCE_POS
     std::sort( thetaOrder, thetaOrder + numFaceNodes );
 
     // Reorder nodes on face
-    for( localIndex n =0; n < numFaceNodes; ++n )
+    for( localIndex n = 0; n < numFaceNodes; ++n )
     {
       faceNodes[n] = thetaOrder[n].second;
     }
@@ -298,7 +283,7 @@ void FaceManager::sortFaceNodes( arrayView2d< real64 const, nodes::REFERENCE_POS
     localIndex tempFaceNodes[MAX_FACE_NODES];
 
     localIndex firstIndexIndex = 0;
-    for( localIndex n =0; n < numFaceNodes; ++n )
+    for( localIndex n = 0; n < numFaceNodes; ++n )
     {
       tempFaceNodes[n] = thetaOrder[n].second;
       if( tempFaceNodes[n] == firstNodeIndex )
@@ -307,9 +292,9 @@ void FaceManager::sortFaceNodes( arrayView2d< real64 const, nodes::REFERENCE_POS
       }
     }
 
-    for( localIndex n=0; n < numFaceNodes; ++n )
+    for( localIndex n = 0; n < numFaceNodes; ++n )
     {
-      const localIndex index = (firstIndexIndex + n) % numFaceNodes;
+      localIndex const index = (firstIndexIndex + n) % numFaceNodes;
       faceNodes[n] = tempFaceNodes[index];
     }
   }
@@ -382,19 +367,19 @@ localIndex FaceManager::packUpDownMapsImpl( buffer_unit_type * & buffer,
   packedSize += bufferOps::Pack< DO_PACKING >( buffer, string( viewKeyStruct::nodeListString() ) );
 
   packedSize += bufferOps::Pack< DO_PACKING >( buffer,
-                                               m_nodeList.base().toViewConst(),
+                                               m_toNodesRelation.base().toViewConst(),
                                                m_unmappedGlobalIndicesInToNodes,
                                                packList,
                                                this->localToGlobalMap(),
-                                               m_nodeList.relatedObjectLocalToGlobal() );
+                                               m_toNodesRelation.relatedObjectLocalToGlobal() );
 
   packedSize += bufferOps::Pack< DO_PACKING >( buffer, string( viewKeyStruct::edgeListString() ) );
   packedSize += bufferOps::Pack< DO_PACKING >( buffer,
-                                               m_edgeList.base().toViewConst(),
+                                               m_toEdgesRelation.base().toViewConst(),
                                                m_unmappedGlobalIndicesInToEdges,
                                                packList,
                                                this->localToGlobalMap(),
-                                               m_edgeList.relatedObjectLocalToGlobal() );
+                                               m_toEdgesRelation.relatedObjectLocalToGlobal() );
 
   packedSize += bufferOps::Pack< DO_PACKING >( buffer, string( viewKeyStruct::elementListString() ) );
   packedSize += bufferOps::Pack< DO_PACKING >( buffer,
@@ -419,22 +404,22 @@ localIndex FaceManager::unpackUpDownMaps( buffer_unit_type const * & buffer,
   GEOSX_ERROR_IF_NE( nodeListString, viewKeyStruct::nodeListString() );
 
   unPackedSize += bufferOps::Unpack( buffer,
-                                     m_nodeList,
+                                     m_toNodesRelation,
                                      packList,
                                      m_unmappedGlobalIndicesInToNodes,
                                      this->globalToLocalMap(),
-                                     m_nodeList.relatedObjectGlobalToLocal() );
+                                     m_toNodesRelation.relatedObjectGlobalToLocal() );
 
   string edgeListString;
   unPackedSize += bufferOps::Unpack( buffer, edgeListString );
   GEOSX_ERROR_IF_NE( edgeListString, viewKeyStruct::edgeListString() );
 
   unPackedSize += bufferOps::Unpack( buffer,
-                                     m_edgeList,
+                                     m_toEdgesRelation,
                                      packList,
                                      m_unmappedGlobalIndicesInToEdges,
                                      this->globalToLocalMap(),
-                                     m_edgeList.relatedObjectGlobalToLocal() );
+                                     m_toEdgesRelation.relatedObjectGlobalToLocal() );
 
   string elementListString;
   unPackedSize += bufferOps::Unpack( buffer, elementListString );
@@ -451,19 +436,19 @@ localIndex FaceManager::unpackUpDownMaps( buffer_unit_type const * & buffer,
 
 void FaceManager::fixUpDownMaps( bool const clearIfUnmapped )
 {
-  ObjectManagerBase::fixUpDownMaps( m_nodeList,
+  ObjectManagerBase::fixUpDownMaps( m_toNodesRelation,
                                     m_unmappedGlobalIndicesInToNodes,
                                     clearIfUnmapped );
 
-  ObjectManagerBase::fixUpDownMaps( m_edgeList,
+  ObjectManagerBase::fixUpDownMaps( m_toEdgesRelation,
                                     m_unmappedGlobalIndicesInToEdges,
                                     clearIfUnmapped );
 }
 
 void FaceManager::compressRelationMaps()
 {
-  m_nodeList.compress();
-  m_edgeList.compress();
+  m_toNodesRelation.compress();
+  m_toEdgesRelation.compress();
 }
 
 void FaceManager::enforceStateFieldConsistencyPostTopologyChange( std::set< localIndex > const & targetIndices )
diff --git a/src/coreComponents/mesh/FaceManager.hpp b/src/coreComponents/mesh/FaceManager.hpp
index 97c471d4ef4..3a1586f225c 100644
--- a/src/coreComponents/mesh/FaceManager.hpp
+++ b/src/coreComponents/mesh/FaceManager.hpp
@@ -61,15 +61,15 @@ class FaceManager : public ObjectManagerBase
    * @brief Return the name of the FaceManager in the object catalog.
    * @return string that contains the catalog name of the FaceManager
    */
-  static const string catalogName()
+  static string catalogName()
   { return "FaceManager"; }
 
   /**
    * @brief Provide a virtual access to catalogName().
    * @return string that contains the catalog name of the FaceManager
    */
-  virtual const string getCatalogName() const override
-  { return FaceManager::catalogName(); }
+  virtual string getCatalogName() const override
+  { return catalogName(); }
   ///@}
 
   /**
@@ -105,28 +105,21 @@ class FaceManager : public ObjectManagerBase
   ///@}
 
   /**
-   * @brief Extend base class resize method resizing  m_nodeList, m_edgeList member containers.
+   * @brief Extend base class resize method resizing  m_toNodesRelation, m_toEdgesRelation member containers.
    * @details the \p newSize of this FaceManager is the number of faces it will contain
    * @param[in] newsize new size the FaceManager.
    */
   virtual void resize( localIndex const newsize ) override;
 
-  /**
-   * @brief Builds the faces to regions and faces to sub-regions mappings.
-   * @param [in] elementRegionManager the ElementRegionManager.
-   *
-   * @note Requires the sub-regions of the @p elementRegionManager to be fully defined.
-   * As well as the faces to elements mappings of the @p FaceManager.
-   */
-  void buildRegionMaps( ElementRegionManager const & elementRegionManager );
-
   /**
    * @brief Copies the nodes positions and the faces to (nodes|edges|elements) mappings from @p cellBlockManager.
    * Computes the faces center, area and normal too.
    * @param[in] cellBlockManager Provides the mappings.
+   * @param[in] elemRegionManager element region manager, needed to map blocks to subregion
    * @param[in] nodeManager Provides the nodes positions.
    */
   void setGeometricalRelations( CellBlockManagerABC const & cellBlockManager,
+                                ElementRegionManager const & elemRegionManager,
                                 NodeManager const & nodeManager );
 
   /**
@@ -159,13 +152,7 @@ class FaceManager : public ObjectManagerBase
   void buildSets( NodeManager const & nodeManager );
 
   /**
-   * @brief Return the number of nodes of the faces with the greatest number of nodes.
-   * @return the maximum number of nodes a face have
-   */
-  localIndex getMaxFaceNodes() const;
-
-  /**
-   * @brief Sort all faces nodes counterclockwisely in m_nodeList.
+   * @brief Sort all faces nodes counterclockwisely in m_toNodesRelation.
    * @param[in] nodeManager node manager allowing access to face nodes coordinates
    * @param[in] elemManager element manager allowing access to the cell elements
    */
@@ -177,12 +164,10 @@ class FaceManager : public ObjectManagerBase
    * @param[in] X array view of mesh nodes coordinates
    * @param[in] elementCenter coordinate of the element center
    * @param[in,out] faceNodes reordered local label list of nodes
-   * @param[in] numFaceNodes number of nodes for the face
    */
-  void sortFaceNodes( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X,
-                      arraySlice1d< real64 const > const elementCenter,
-                      localIndex * const faceNodes,
-                      localIndex const numFaceNodes );
+  static void sortFaceNodes( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X,
+                             arraySlice1d< real64 const > const elementCenter,
+                             Span< localIndex > const faceNodes );
 
   /**
    * @brief Flag faces on boundary or external to the DomainPartition.
@@ -348,25 +333,25 @@ class FaceManager : public ObjectManagerBase
    * @brief Get a mutable accessor to a map containing the list of each nodes for each faces.
    * @return non-const reference to a map containing the list of each nodes for each faces
    */
-  NodeMapType & nodeList()                    { return m_nodeList; }
+  NodeMapType & nodeList()                    { return m_toNodesRelation; }
 
   /**
    * @brief Get an immutable accessor to a map containing the list of each nodes for each faces
    * @return const reference to a map containing the list of each nodes for each faces
    */
-  NodeMapType const & nodeList() const { return m_nodeList; }
+  NodeMapType const & nodeList() const { return m_toNodesRelation; }
 
   /**
    * @brief Get a mutable accessor to a map containing the list of each edges for each faces
    * @return non-const reference to a map containing the list of each edges for each faces
    */
-  EdgeMapType & edgeList()       { return m_edgeList; }
+  EdgeMapType & edgeList()       { return m_toEdgesRelation; }
 
   /**
    * @brief Get an immutable accessor to a map containing the list of each edges for each faces.
    * @return const reference to a map containing the list of each edges for each faces
    */
-  EdgeMapType const & edgeList() const { return m_edgeList; }
+  EdgeMapType const & edgeList() const { return m_toEdgesRelation; }
 
   /**
    * @brief Get a mutable accessor to the faces-to-ElementRegion relation.
@@ -455,10 +440,10 @@ class FaceManager : public ObjectManagerBase
                                  arrayView1d< localIndex const > const & packList ) const;
 
   /// face keyed map containing face-to-node relation
-  NodeMapType m_nodeList;
+  NodeMapType m_toNodesRelation;
 
   /// face keyed map containing face-to-edge relation
-  EdgeMapType m_edgeList;
+  EdgeMapType m_toEdgesRelation;
 
   /// face keyed map containing face-to-element relation
   ElemMapType m_toElements;
diff --git a/src/coreComponents/mesh/NodeManager.cpp b/src/coreComponents/mesh/NodeManager.cpp
index 6d58941afab..f20dac03295 100644
--- a/src/coreComponents/mesh/NodeManager.cpp
+++ b/src/coreComponents/mesh/NodeManager.cpp
@@ -52,10 +52,6 @@ NodeManager::NodeManager( string const & name,
 }
 
 
-NodeManager::~NodeManager()
-{}
-
-
 void NodeManager::resize( localIndex const newSize )
 {
   m_toFacesRelation.resize( newSize, 2 * getFaceMapOverallocation() );
@@ -73,25 +69,11 @@ void NodeManager::constructGlobalToLocalMap( CellBlockManagerABC const & cellBlo
   ObjectManagerBase::constructGlobalToLocalMap();
 }
 
-
-void NodeManager::buildRegionMaps( ElementRegionManager const & elementRegionManager )
-{
-  GEOSX_MARK_FUNCTION;
-
-  auto const nodeListGetter =
-    []( CellElementSubRegion const & subregion ) -> CellElementSubRegion::NodeMapType const &
-  { return subregion.nodeList(); };
-
-  meshMapUtilities::buildToElementMaps( elementRegionManager,
-                                        size(),
-                                        m_toElements,
-                                        nodeListGetter,
-                                        getElemMapOverAllocation() );
-}
-
 void NodeManager::buildSets( CellBlockManagerABC const & cellBlockManager,
                              GeometricObjectManager const & geometries )
 {
+  GEOSX_MARK_FUNCTION;
+
   // Let's first copy the sets from the cell block manager.
   for( const auto & nameArray: cellBlockManager.getNodeSets() )
   {
@@ -103,8 +85,7 @@ void NodeManager::buildSets( CellBlockManagerABC const & cellBlockManager,
   arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = this->referencePosition();
   localIndex const numNodes = this->size();
 
-  geometries.forSubGroups< SimpleGeometricObjectBase >(
-    [&]( SimpleGeometricObjectBase const & object ) -> void
+  geometries.forSubGroups< SimpleGeometricObjectBase >( [&]( SimpleGeometricObjectBase const & object )
   {
     string const & name = object.getName();
     SortedArray< localIndex > & targetSet = m_sets.registerWrapper< SortedArray< localIndex > >( name ).reference();
@@ -121,33 +102,43 @@ void NodeManager::buildSets( CellBlockManagerABC const & cellBlockManager,
 
 void NodeManager::setDomainBoundaryObjects( FaceManager const & faceManager )
 {
-  arrayView1d< integer const > const & isFaceOnDomainBoundary = faceManager.getDomainBoundaryIndicator();
-  arrayView1d< integer > const & isNodeOnDomainBoundary = getDomainBoundaryIndicator();
+  arrayView1d< integer const > const isFaceOnDomainBoundary = faceManager.getDomainBoundaryIndicator();
+  arrayView1d< integer > const isNodeOnDomainBoundary = getDomainBoundaryIndicator();
   isNodeOnDomainBoundary.zero();
 
   ArrayOfArraysView< localIndex const > const faceToNodes = faceManager.nodeList().toViewConst();
 
-  forAll< parallelHostPolicy >( faceManager.size(), [&]( localIndex const k )
+  forAll< parallelHostPolicy >( faceManager.size(), [=]( localIndex const faceIndex )
   {
-    if( isFaceOnDomainBoundary[k] == 1 )
+    if( isFaceOnDomainBoundary[faceIndex] == 1 )
     {
-      localIndex const numNodes = faceToNodes.sizeOfArray( k );
-      for( localIndex a = 0; a < numNodes; ++a )
+      for( localIndex const nodeIndex : faceToNodes[faceIndex] )
       {
-        isNodeOnDomainBoundary[faceToNodes( k, a )] = 1;
+        isNodeOnDomainBoundary[nodeIndex] = 1;
       }
     }
   } );
 }
 
-void NodeManager::setGeometricalRelations( CellBlockManagerABC const & cellBlockManager )
+void NodeManager::setGeometricalRelations( CellBlockManagerABC const & cellBlockManager,
+                                           ElementRegionManager const & elemRegionManager )
 {
+  GEOSX_MARK_FUNCTION;
+
   resize( cellBlockManager.numNodes() );
 
-  m_referencePosition = cellBlockManager.getNodesPositions();
+  m_referencePosition = cellBlockManager.getNodePositions();
+
+  m_toEdgesRelation.base().assimilate< parallelHostPolicy >( cellBlockManager.getNodeToEdges(),
+                                                             LvArray::sortedArrayManipulation::UNSORTED_NO_DUPLICATES );
+  m_toFacesRelation.base().assimilate< parallelHostPolicy >( cellBlockManager.getNodeToFaces(),
+                                                             LvArray::sortedArrayManipulation::UNSORTED_NO_DUPLICATES );
 
-  m_toEdgesRelation.base() = cellBlockManager.getNodeToEdges();
-  m_toFacesRelation.base() = cellBlockManager.getNodeToFaces();
+  ToCellRelation< ArrayOfArrays< localIndex > > const toCellBlock = cellBlockManager.getNodeToElements();
+  array2d< localIndex > const blockToSubRegion = elemRegionManager.getCellBlockToSubRegionMap( cellBlockManager );
+  meshMapUtilities::transformCellBlockToRegionMap< parallelHostPolicy >( blockToSubRegion.toViewConst(),
+                                                                         toCellBlock,
+                                                                         m_toElements );
 }
 
 void NodeManager::setupRelatedObjectsInRelations( EdgeManager const & edgeManager,
diff --git a/src/coreComponents/mesh/NodeManager.hpp b/src/coreComponents/mesh/NodeManager.hpp
index c07ed17ba8a..1e91871399d 100644
--- a/src/coreComponents/mesh/NodeManager.hpp
+++ b/src/coreComponents/mesh/NodeManager.hpp
@@ -82,9 +82,9 @@ class NodeManager : public ObjectManagerBase
   static constexpr localIndex getElemMapOverAllocation()
   { return CellBlockManagerABC::elemMapExtraSpacePerNode(); }
 
-/**
- * @name Constructors/destructor
- */
+  /**
+   * @name Constructors/destructor
+   */
   ///@{
 
   /**
@@ -95,28 +95,6 @@ class NodeManager : public ObjectManagerBase
   NodeManager( string const & name,
                dataRepository::Group * const parent );
 
-  /**
-   * @brief The default NodeManager destructor.
-   */
-  ~NodeManager() override;
-
-  /// @cond DO_NOT_DOCUMENT
-  /**
-   * @brief deleted constructor
-   */
-  NodeManager() = delete;
-
-  /**
-   * @brief deleted copy constructor
-   */
-  NodeManager( const NodeManager & init ) = delete;
-
-  /**
-   * @brief deleted assignement operator
-   */
-  NodeManager & operator=( const NodeManager & ) = delete;
-  /// @endcond
-
   ///@}
 
   /**
@@ -142,20 +120,11 @@ class NodeManager : public ObjectManagerBase
    * @brief Provide a virtual access to catalogName().
    * @return string that contains the NodeManager catalog name
    */
-  const string getCatalogName() const override final
-  { return NodeManager::catalogName(); }
+  string getCatalogName() const override final
+  { return catalogName(); }
 
   ///@}
 
-  /**
-   * @brief Builds the nodes to regions and nodes to sub-regions mappings.
-   * @param [in] elementRegionManager the ElementRegionManager.
-   *
-   * @note Requires the sub-regions of the @p elementRegionManager to be fully defined.
-   * As well as the node to elements mappings of the @p NodeManager.
-   */
-  void buildRegionMaps( ElementRegionManager const & elementRegionManager );
-
   /**
    * @brief Copies the local to global mapping from @p cellBlockManager and invert to create the global to local mapping.
    * @param cellBlockManager Provides the local to global mapping.
@@ -172,16 +141,18 @@ class NodeManager : public ObjectManagerBase
 
   /**
    * @brief Builds the node-on-domain-boundary indicator.
-   * @param[in] faceManager The computation is based on the face-on-domain-boundary indicator.
+   * @param[in] faceIndex The computation is based on the face-on-domain-boundary indicator.
    * @see ObjectManagerBase::getDomainBoundaryIndicator()
    */
-  void setDomainBoundaryObjects( FaceManager const & faceManager );
+  void setDomainBoundaryObjects( FaceManager const & faceIndex );
 
   /**
    * @brief Copies the nodes positions and the nodes to (edges|faces|elements) mappings from @p cellBlockManager.
    * @param[in] cellBlockManager Will provide the mappings.
+   * @param[in] elemRegionManager element region manager, needed to map blocks to subregion
    */
-  void setGeometricalRelations( CellBlockManagerABC const & cellBlockManager );
+  void setGeometricalRelations( CellBlockManagerABC const & cellBlockManager,
+                                ElementRegionManager const & elemRegionManager );
 
   /**
    * @brief Link the current manager to other managers.
diff --git a/src/coreComponents/mesh/ObjectManagerBase.cpp b/src/coreComponents/mesh/ObjectManagerBase.cpp
index eb69b4050a2..5feeeb9a99f 100644
--- a/src/coreComponents/mesh/ObjectManagerBase.cpp
+++ b/src/coreComponents/mesh/ObjectManagerBase.cpp
@@ -141,30 +141,28 @@ void ObjectManagerBase::constructSetFromSetAndMap( SortedArrayView< localIndex c
                                                    ArrayOfArraysView< localIndex const > const & map,
                                                    const string & setName )
 {
-  SortedArray< localIndex > & newset = m_sets.getReference< SortedArray< localIndex > >( setName );
-  newset.clear();
+  SortedArray< localIndex > & newSet = m_sets.getReference< SortedArray< localIndex > >( setName );
+  newSet.clear();
 
   localIndex const numObjects = size();
-  GEOSX_ERROR_IF( map.size() != numObjects, "Size mismatch. " << map.size() << " != " << numObjects );
+  GEOSX_ERROR_IF_NE_MSG( map.size(), numObjects, "Map size does not match number of objects." );
 
   if( setName == "all" )
   {
-    newset.reserve( numObjects );
-
-    for( localIndex ka=0; ka<numObjects; ++ka )
+    newSet.reserve( numObjects );
+    for( localIndex ka = 0; ka < numObjects; ++ka )
     {
-      newset.insert( ka );
+      newSet.insert( ka );
     }
   }
   else
   {
-    for( localIndex ka=0; ka<numObjects; ++ka )
+    for( localIndex ka = 0; ka < numObjects; ++ka )
     {
-      localIndex const * const values = map[ka];
-      localIndex const numValues = map.sizeOfArray( ka );
-      if( std::all_of( values, values + numValues, [&]( localIndex const i ) { return inputSet.contains( i ); } ) )
+      arraySlice1d< localIndex const > const values = map[ka];
+      if( std::all_of( values.begin(), values.end(), [&]( localIndex const i ) { return inputSet.contains( i ); } ) )
       {
-        newset.insert( ka );
+        newSet.insert( ka );
       }
     }
   }
diff --git a/src/coreComponents/mesh/ObjectManagerBase.hpp b/src/coreComponents/mesh/ObjectManagerBase.hpp
index 00338e6e3f9..b5aa09558dd 100644
--- a/src/coreComponents/mesh/ObjectManagerBase.hpp
+++ b/src/coreComponents/mesh/ObjectManagerBase.hpp
@@ -69,7 +69,7 @@ class ObjectManagerBase : public dataRepository::Group
    * @brief Get the name of the catalog.
    * @return The name.
    */
-  virtual const string getCatalogName() const = 0;
+  virtual string getCatalogName() const = 0;
   ///@}
 
   using dataRepository::Group::packSize;
diff --git a/src/coreComponents/mesh/PerforationData.hpp b/src/coreComponents/mesh/PerforationData.hpp
index 68204a501ca..e3d589c7d0b 100644
--- a/src/coreComponents/mesh/PerforationData.hpp
+++ b/src/coreComponents/mesh/PerforationData.hpp
@@ -99,7 +99,7 @@ class PerforationData : public ObjectManagerBase
   /**
    * @copydoc catalogName()
    */
-  virtual const string getCatalogName() const override { return catalogName(); }
+  virtual string getCatalogName() const override { return catalogName(); }
 
   ///@}
 
diff --git a/src/coreComponents/mesh/SurfaceElementRegion.hpp b/src/coreComponents/mesh/SurfaceElementRegion.hpp
index 7533c6d60f0..d7119f9b1c8 100644
--- a/src/coreComponents/mesh/SurfaceElementRegion.hpp
+++ b/src/coreComponents/mesh/SurfaceElementRegion.hpp
@@ -84,11 +84,11 @@ class SurfaceElementRegion : public ElementRegionBase
    * @brief Get the key name for the SurfaceElementRegion in the object catalog.
    * @return A string containing the key name.
    */
-  static const string catalogName()
+  static string catalogName()
   { return "SurfaceElementRegion"; }
 
-  virtual const string getCatalogName() const override final
-  { return SurfaceElementRegion::catalogName(); }
+  virtual string getCatalogName() const override final
+  { return catalogName(); }
 
   ///@}
 
diff --git a/src/coreComponents/mesh/SurfaceElementSubRegion.hpp b/src/coreComponents/mesh/SurfaceElementSubRegion.hpp
index caafabe50eb..58b1fc0e56b 100644
--- a/src/coreComponents/mesh/SurfaceElementSubRegion.hpp
+++ b/src/coreComponents/mesh/SurfaceElementSubRegion.hpp
@@ -53,16 +53,16 @@ class SurfaceElementSubRegion : public ElementSubRegionBase
    * @brief Get catalog name.
    * @return the catalog name
    */
-  static const string catalogName()
+  static string catalogName()
   { return "SurfaceElementSubRegion"; }
 
   /**
    * @brief Get catalog name.
    * @return the catalog name
    */
-  virtual const string getCatalogName() const override
+  virtual string getCatalogName() const override
   {
-    return SurfaceElementSubRegion::catalogName();
+    return catalogName();
   }
 
   ///@}
diff --git a/src/coreComponents/mesh/WellElementRegion.hpp b/src/coreComponents/mesh/WellElementRegion.hpp
index 213a7aaeb9f..bd1b28c7eeb 100644
--- a/src/coreComponents/mesh/WellElementRegion.hpp
+++ b/src/coreComponents/mesh/WellElementRegion.hpp
@@ -71,14 +71,14 @@ class WellElementRegion : public ElementRegionBase
    * @brief Get the catalog name.
    * @return the name of this class in the catalog
    */
-  static const string catalogName()
+  static string catalogName()
   { return "WellElementRegion"; }
 
   /**
    * @copydoc catalogName()
    */
-  virtual const string getCatalogName() const override final
-  { return WellElementRegion::catalogName(); }
+  virtual string getCatalogName() const override final
+  { return catalogName(); }
 
   ///@}
 
diff --git a/src/coreComponents/mesh/WellElementSubRegion.hpp b/src/coreComponents/mesh/WellElementSubRegion.hpp
index 688600112ae..01821247690 100644
--- a/src/coreComponents/mesh/WellElementSubRegion.hpp
+++ b/src/coreComponents/mesh/WellElementSubRegion.hpp
@@ -72,12 +72,12 @@ class WellElementSubRegion : public ElementSubRegionBase
    * @brief Get the catalog name.
    * @return the name of this class in the catalog
    */
-  static const string catalogName() { return "wellElementSubRegion"; }
+  static string catalogName() { return "wellElementSubRegion"; }
 
   /**
    * @copydoc catalogName()
    */
-  virtual const string getCatalogName() const override { return WellElementSubRegion::catalogName(); }
+  virtual string getCatalogName() const override { return catalogName(); }
 
   ///@}
 
diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp
index 8f85576b7a3..216f2c29e3b 100644
--- a/src/coreComponents/mesh/generators/CellBlockManager.cpp
+++ b/src/coreComponents/mesh/generators/CellBlockManager.cpp
@@ -45,53 +45,152 @@ Group * CellBlockManager::createChild( string const & GEOSX_UNUSED_PARAM( childK
   return nullptr;
 }
 
-ArrayOfArrays< localIndex > CellBlockManager::getNodeToElements() const
+/// Element identifier containing (block index, cell index).
+using CellBlockIndexPair = std::pair< localIndex, localIndex >;
+
+template< typename POLICY >
+void convertFromCellBlockPairMap( ArrayOfArraysView< CellBlockIndexPair const > const & srcMap,
+                                  ToCellRelation< ArrayOfArrays< localIndex > > & dstMap )
 {
-  // The function works in three steps.
-  // First finds the number of elements attached to each node.
-  // Then second step allocates the vectors, including extra allocations.
-  // Last, the output vector is filled.
+  ArrayOfArrays< localIndex > & toBlock = dstMap.toBlockIndex;
+  ArrayOfArrays< localIndex > & toCell = dstMap.toCellIndex;
+
+  localIndex const numObjects = srcMap.size();
+
+  toBlock.resizeFromOffsets( numObjects, srcMap.toViewConst().getOffsets() );
+  toCell.resizeFromOffsets( numObjects, srcMap.toViewConst().getOffsets() );
+
+  forAll< parallelHostPolicy >( numObjects, [toBlock = toBlock.toView(),
+                                             toCell = toCell.toView(),
+                                             srcMap]( localIndex const objIndex )
+  {
+    arraySlice1d< CellBlockIndexPair const > const cells = srcMap[ objIndex ];
+    for( CellBlockIndexPair const & e : cells )
+    {
+      toBlock.emplaceBack( objIndex, std::get< 0 >( e ) );
+      toCell.emplaceBack( objIndex, std::get< 1 >( e ) );
+    }
+  } );
+}
+
+template< typename POLICY >
+void convertFromCellBlockPairMap( ArrayOfArraysView< CellBlockIndexPair const > const & srcMap,
+                                  ToCellRelation< array2d< localIndex > > & dstMap )
+{
+  array2d< localIndex > & toBlock = dstMap.toBlockIndex;
+  array2d< localIndex > & toCell = dstMap.toCellIndex;
+
+  localIndex const numObjects = srcMap.size();
+  localIndex const maxNumElem = toCell.size( 1 );
+
+  GEOSX_ERROR_IF_NE( toBlock.size( 1 ), maxNumElem );
+  GEOSX_ERROR_IF_NE( toCell.size( 1 ), maxNumElem );
+
+  toBlock.resizeDimension< 0 >( numObjects );
+  toCell.resizeDimension< 0 >( numObjects );
+
+  // We allow a fixed-size map to represent a variable relationship, as long as
+  // the number of elements does not exceed the fixed size (set by the caller).
+  // In this case, a dummy value "-1" is used to represent non-present elements.
+  toBlock.setValues< POLICY >( -1 );
+  toCell.setValues< POLICY >( -1 );
+
+  forAll< parallelHostPolicy >( numObjects, [=, // needed to optionally capture maxNumElem in Debug
+                                             toBlock = toBlock.toView(),
+                                             toCell = toCell.toView()]( localIndex const objIndex )
+  {
+    arraySlice1d< CellBlockIndexPair const > const cells = srcMap[ objIndex ];
+    GEOSX_ASSERT_GE( maxNumElem, cells.size() );
+    localIndex count = 0;
+    for( CellBlockIndexPair const & e : cells )
+    {
+      toBlock( objIndex, count ) = std::get< 0 >( e );
+      toCell( objIndex, count ) = std::get< 1 >( e );
+      ++count;
+    }
+  } );
+}
 
-  // First step: how many elements for each node, stored in the elemsPerNode array.
-  array1d< localIndex > elemsPerNode( m_numNodes );
-  elemsPerNode.setValues< serialPolicy >( elemMapExtraSpacePerNode() );
+/**
+ * @brief Build to-cell map by inverting existing maps in cell blocks.
+ * @tparam BASEMAP underlying type of to-cell map
+ * @tparam FUNC type of @p cellToObjectGetter
+ * @param cellIndex number of objects (nodes, faces, etc.) for which maps are built
+ * @param toCells container for to-cell (block, index) maps
+ * @param cellToObjectGetter function used to extract maps from subregions
+ * @param overAlloc overallocation for the resulting maps
+ *                  (extra capacity per row, only meaningful for variable-secondary-size containers)
+ */
+template< typename BASEMAP, typename FUNC >
+void CellBlockManager::buildToCellMap( localIndex const numObjects,
+                                       ToCellRelation< BASEMAP > & toCells,
+                                       FUNC cellToObjectGetter,
+                                       localIndex const overAlloc ) const
+{
+  // Calculate the number of entries in each sub-array
+  array1d< localIndex > counts( numObjects );
+  counts.setValues< serialPolicy >( overAlloc );
 
   for( localIndex blockIndex = 0; blockIndex < numCellBlocks(); ++blockIndex )
   {
-    CellBlock const & cb = this->getCellBlock( blockIndex );
-    arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemToNode = cb.getElemToNode();
-    forAll< parallelHostPolicy >( cb.numElements(), [elemsPerNode = elemsPerNode.toView(),
-                                                     elemToNode, &cb] ( localIndex const cellIndex )
+    CellBlock const & cb = getCellBlock( blockIndex );
+
+    auto const cellToObject = cellToObjectGetter( cb );
+    forAll< parallelHostPolicy >( cb.size(), [counts = counts.toView(),
+                                              cellToObject]( localIndex const ei )
     {
-      localIndex const numNodesPerElement = cb.numNodesPerElement();
-      for( localIndex a = 0; a < numNodesPerElement; ++a )
+      auto const objects = cellToObject[ ei ];
+      // can't use range-based for loop when slice is not contiguous
+      for( localIndex i = 0; i < objects.size(); ++i )
       {
-        localIndex const nodeIndex = elemToNode( cellIndex, a );
-        RAJA::atomicInc< parallelHostAtomic >( &elemsPerNode[ nodeIndex ] );
+        RAJA::atomicInc< parallelHostAtomic >( &counts[ objects[i] ] );
       }
     } );
   }
+  ;
 
-  // Second, allocation.
-  ArrayOfArrays< localIndex > result;
-  result.resizeFromCapacities< parallelHostPolicy >( m_numNodes, elemsPerNode.data() );
+  // Allocate memory
+  ArrayOfArrays< CellBlockIndexPair > cellBlockPairList;
+  cellBlockPairList.resizeFromCapacities< parallelHostPolicy >( numObjects, counts.data() );
 
-  // Third, filling the result.
-  for( localIndex blockIndex = 0; blockIndex < numCellBlocks(); ++blockIndex )// Do not need index
+  // Populate map of tuples in parallel
+  for( localIndex blockIndex = 0; blockIndex < numCellBlocks(); ++blockIndex )
   {
-    CellBlock const & cb = this->getCellBlock( blockIndex );
-    arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemToNode = cb.getElemToNode();
-    forAll< parallelHostPolicy >( cb.numElements(), [result = result.toView(),
-                                                     elemToNode, &cb] ( localIndex const cellIndex )
+    CellBlock const & cb = getCellBlock( blockIndex );
+
+    auto const cellToObject = cellToObjectGetter( cb );
+    forAll< parallelHostPolicy >( cb.size(), [cellBlockPairList = cellBlockPairList.toView(),
+                                              cellToObject,
+                                              blockIndex]( localIndex const cellIndex )
     {
-      for( localIndex a = 0; a < cb.numNodesPerElement(); ++a )
+      auto const objects = cellToObject[ cellIndex ];
+      // can't use range-based for loop when slice is not contiguous
+      for( localIndex i = 0; i < objects.size(); ++i )
       {
-        localIndex const nodeIndex = elemToNode( cellIndex, a );
-        result.emplaceBackAtomic< parallelHostAtomic >( nodeIndex, cellIndex );
+        cellBlockPairList.emplaceBackAtomic< parallelHostAtomic >( objects[i], blockIndex, cellIndex );
       }
     } );
   }
+  ;
 
+  // Sort each element list to ensure unique race-condition-free map order
+  forAll< parallelHostPolicy >( numObjects, [cellBlockPairList = cellBlockPairList.toView()]( localIndex const objIndex )
+  {
+    arraySlice1d< CellBlockIndexPair > const cells = cellBlockPairList[ objIndex ];
+    LvArray::sortedArrayManipulation::makeSorted( cells.begin(), cells.end() );
+  } );
+
+  // Finally, split arrays-of-tuples into separate arrays
+  convertFromCellBlockPairMap< parallelHostPolicy >( cellBlockPairList.toViewConst(), toCells );
+}
+
+ToCellRelation< ArrayOfArrays< localIndex > > CellBlockManager::getNodeToElements() const
+{
+  ToCellRelation< ArrayOfArrays< localIndex > > result;
+  buildToCellMap( m_numNodes,
+                  result,
+                  []( CellBlock const & cb ) { return cb.getElemToNode(); },
+                  elemMapExtraSpacePerNode() );
   return result;
 }
 
@@ -164,27 +263,31 @@ struct NodesAndElementOfFace
  * @brief Fills the face to nodes map and face to element maps
  * @param [in] lowestNodeToFaces and array of size numNodes of arrays of NodesAndElementOfFace associated with each node.
  * @param [in] uniqueFaceOffsets an array containing the unique ID of the first face associated with each node.
- * @param [inout] faceToElements the face to element map.
+ * @param [inout] faceToCells the face to element map.
  * @param [inout] faceToNodes the face to node map.
  */
 void populateFaceMaps( Group const & cellBlocks,
                        ArrayOfArraysView< NodesAndElementOfFace const > const & lowestNodeToFaces,
                        arrayView1d< localIndex const > const & uniqueFaceOffsets,
                        ArrayOfArraysView< localIndex > const & faceToNodes,
-                       arrayView2d< localIndex > const & faceToElements )
+                       arrayView2d< localIndex > const & faceToCells,
+                       arrayView2d< localIndex > const & faceToBlocks )
 {
   localIndex const numNodes = lowestNodeToFaces.size();
   localIndex const numUniqueFaces = uniqueFaceOffsets.back();
   GEOSX_ERROR_IF_NE( uniqueFaceOffsets.size() - 1, numNodes );
   GEOSX_ERROR_IF_NE( faceToNodes.size(), numUniqueFaces );
-  GEOSX_ERROR_IF_NE( faceToElements.size( 0 ), numUniqueFaces );
-  GEOSX_ERROR_IF_NE( faceToElements.size( 1 ), 2 );
+  GEOSX_ERROR_IF_NE( faceToCells.size( 0 ), numUniqueFaces );
+  GEOSX_ERROR_IF_NE( faceToCells.size( 1 ), 2 );
+  GEOSX_ERROR_IF_NE( faceToBlocks.size( 0 ), numUniqueFaces );
+  GEOSX_ERROR_IF_NE( faceToBlocks.size( 1 ), 2 );
 
   // loop over all the nodes.
   forAll< parallelHostPolicy >( numNodes, [uniqueFaceOffsets,
                                            lowestNodeToFaces,
                                            faceToNodes,
-                                           faceToElements,
+                                           faceToCells,
+                                           faceToBlocks,
                                            &cellBlocks]( localIndex const nodeIndex )
   {
     localIndex nodesInFace[ CellBlockManager::maxNodesPerFace() ];
@@ -202,13 +305,19 @@ void populateFaceMaps( Group const & cellBlocks,
         faceToNodes.emplaceBack( curFaceID, nodesInFace[i] );
       }
 
-      faceToElements( curFaceID, 0 ) = f0.cellIndex;
-      faceToElements( curFaceID, 1 ) = -1;
+      faceToCells( curFaceID, 0 ) = f0.cellIndex;
+      faceToBlocks( curFaceID, 0 ) = f0.blockIndex;
 
       if( ++first != last )
       {
         NodesAndElementOfFace const & f1 = *first++;
-        faceToElements( curFaceID, 1 ) = f1.cellIndex;
+        faceToCells( curFaceID, 1 ) = f1.cellIndex;
+        faceToBlocks( curFaceID, 1 ) = f1.blockIndex;
+      }
+      else
+      {
+        faceToCells( curFaceID, 1 ) = -1;
+        faceToBlocks( curFaceID, 1 ) = -1;
       }
       GEOSX_ASSERT( first == last ); // Should not be more than 2 faces
 
@@ -223,13 +332,14 @@ void populateFaceMaps( Group const & cellBlocks,
  * @param [in] lowestNodeToFaces and array of size numNodes of arrays of NodesAndElementOfFace associated with each node.
  * @param [in] uniqueFaceOffsets an containing the unique face IDs for each node in lowestNodeToFaces.
  * @param [out] faceToNodeMap the map from faces to nodes. This function resizes the array appropriately.
- * @param [out] faceToElemMap the map from faces to elements. This function resizes the array appropriately.
+ * @param [out] faceToCellMap the map from faces to elements. This function resizes the array appropriately.
  */
 void resizeFaceMaps( ArrayOfArraysView< NodesAndElementOfFace const > const & lowestNodeToFaces,
                      arrayView1d< localIndex const > const & uniqueFaceOffsets,
                      ArrayOfArrays< localIndex > & faceToNodeMap,
                      ArrayOfArrays< localIndex > & faceToEdgesMap,
-                     array2d< localIndex > & faceToElemMap )
+                     array2d< localIndex > & faceToCellMap,
+                     array2d< localIndex > & faceToBlockMap )
 {
   localIndex const numNodes = lowestNodeToFaces.size();
   localIndex const numUniqueFaces = uniqueFaceOffsets.back();
@@ -252,7 +362,8 @@ void resizeFaceMaps( ArrayOfArraysView< NodesAndElementOfFace const > const & lo
 
   // Each face may belong to _maximum_ 2 elements.
   // If it belongs to only one, we put `-1` for the undefined value.
-  faceToElemMap.resize( numUniqueFaces, 2 );
+  faceToCellMap.resize( numUniqueFaces, 2 );
+  faceToBlockMap.resize( numUniqueFaces, 2 );
 
   // TODO I'm really not sure.
   faceToEdgesMap.resize( numUniqueFaces, 2 * CellBlockManager::edgeMapExtraSpacePerFace() );
@@ -443,13 +554,15 @@ void CellBlockManager::buildFaceMaps()
                   uniqueFaceOffsets,
                   m_faceToNodes,
                   m_faceToEdges,
-                  m_faceToElements );
+                  m_faceToCells.toCellIndex,
+                  m_faceToCells.toBlockIndex );
 
   populateFaceMaps( getCellBlocks(),
                     lowestNodeToFaces.toViewConst(),
                     uniqueFaceOffsets,
                     m_faceToNodes.toView(),
-                    m_faceToElements );
+                    m_faceToCells.toCellIndex,
+                    m_faceToCells.toBlockIndex );
 
   fillElementToFacesOfCellBlocks( lowestNodeToFaces.toViewConst(),
                                   uniqueFaceOffsets,
@@ -458,12 +571,9 @@ void CellBlockManager::buildFaceMaps()
 
 void CellBlockManager::buildNodeToEdges()
 {
-  ArrayOfArrays< localIndex > nodeToEdges =
-    meshMapUtilities::transposeIndexMap< parallelHostPolicy >( m_edgeToNodes.toViewConst(),
-                                                               m_numNodes,
-                                                               edgeMapExtraSpacePerNode() );
-  m_nodeToEdges.assimilate< parallelHostPolicy >( std::move( nodeToEdges ),
-                                                  LvArray::sortedArrayManipulation::UNSORTED_NO_DUPLICATES );
+  m_nodeToEdges = meshMapUtilities::transposeIndexMap< parallelHostPolicy >( m_edgeToNodes.toViewConst(),
+                                                                             m_numNodes,
+                                                                             edgeMapExtraSpacePerNode() );
 }
 
 void CellBlockManager::buildMaps()
@@ -486,21 +596,16 @@ ArrayOfArrays< localIndex > CellBlockManager::getFaceToNodes() const
   return m_faceToNodes;
 }
 
-ArrayOfSets< localIndex > CellBlockManager::getNodeToFaces() const
+ArrayOfArrays< localIndex > CellBlockManager::getNodeToFaces() const
 {
-  ArrayOfArrays< localIndex > nodesToFaces =
-    meshMapUtilities::transposeIndexMap< parallelHostPolicy >( m_faceToNodes.toViewConst(),
-                                                               m_numNodes,
-                                                               faceMapExtraSpacePerNode() );
-  ArrayOfSets< localIndex > result;
-  result.assimilate< parallelHostPolicy >( std::move( nodesToFaces ),
-                                           LvArray::sortedArrayManipulation::UNSORTED_NO_DUPLICATES );
-  return result;
+  return meshMapUtilities::transposeIndexMap< parallelHostPolicy >( m_faceToNodes.toViewConst(),
+                                                                    m_numNodes,
+                                                                    faceMapExtraSpacePerNode() );
 }
 
-array2d< localIndex > CellBlockManager::getFaceToElements() const
+ToCellRelation< array2d< localIndex > > CellBlockManager::getFaceToElements() const
 {
-  return m_faceToElements;
+  return m_faceToCells;
 }
 
 const Group & CellBlockManager::getCellBlocks() const
@@ -533,7 +638,7 @@ localIndex CellBlockManager::numFaces() const
   return m_numFaces;
 }
 
-ArrayOfSets< localIndex > CellBlockManager::getEdgeToFaces() const
+ArrayOfArrays< localIndex > CellBlockManager::getEdgeToFaces() const
 {
   return m_edgeToFaces;
 }
@@ -548,9 +653,11 @@ ArrayOfArrays< localIndex > CellBlockManager::getFaceToEdges() const
   return m_faceToEdges;
 }
 
-ArrayOfSets< localIndex > CellBlockManager::getNodeToEdges() const
+ArrayOfArrays< localIndex > CellBlockManager::getNodeToEdges() const
 {
-  return m_nodeToEdges;
+  return meshMapUtilities::transposeIndexMap< parallelHostPolicy >( m_edgeToNodes.toViewConst(),
+                                                                    m_numNodes,
+                                                                    edgeMapExtraSpacePerNode() );
 }
 
 localIndex CellBlockManager::numEdges() const
@@ -563,12 +670,12 @@ CellBlock & CellBlockManager::registerCellBlock( string const & name )
   return this->getCellBlocks().registerGroup< CellBlock >( name );
 }
 
-array2d< real64, nodes::REFERENCE_POSITION_PERM > CellBlockManager::getNodesPositions() const
+array2d< real64, nodes::REFERENCE_POSITION_PERM > CellBlockManager::getNodePositions() const
 {
   return m_nodesPositions;
 }
 
-arrayView2d< real64, nodes::REFERENCE_POSITION_USD > CellBlockManager::getNodesPositions()
+arrayView2d< real64, nodes::REFERENCE_POSITION_USD > CellBlockManager::getNodePositions()
 {
   return m_nodesPositions.toView();
 }
diff --git a/src/coreComponents/mesh/generators/CellBlockManager.hpp b/src/coreComponents/mesh/generators/CellBlockManager.hpp
index e2b514697e4..29ee81982e8 100644
--- a/src/coreComponents/mesh/generators/CellBlockManager.hpp
+++ b/src/coreComponents/mesh/generators/CellBlockManager.hpp
@@ -49,7 +49,7 @@ class CellBlockManager : public CellBlockManagerABC
   static constexpr int maxNodesPerFace()
   { return 64; }
 
-  array2d< real64, nodes::REFERENCE_POSITION_PERM > getNodesPositions() const override;
+  array2d< real64, nodes::REFERENCE_POSITION_PERM > getNodePositions() const override;
 
   /**
    * @brief Returns a view to the vector holding the nodes coordinates
@@ -57,23 +57,23 @@ class CellBlockManager : public CellBlockManagerABC
    *
    * @note This is meant to be used as a values setter.
    */
-  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > getNodesPositions();
+  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > getNodePositions();
 
-  ArrayOfSets< localIndex > getNodeToEdges() const override;
+  ArrayOfArrays< localIndex > getNodeToEdges() const override;
 
-  ArrayOfSets< localIndex > getNodeToFaces() const override;
+  ArrayOfArrays< localIndex > getNodeToFaces() const override;
 
-  ArrayOfArrays< localIndex > getNodeToElements() const override;
+  ToCellRelation< ArrayOfArrays< localIndex > > getNodeToElements() const override;
 
   array2d< localIndex > getEdgeToNodes() const override;
 
-  ArrayOfSets< localIndex > getEdgeToFaces() const override;
+  ArrayOfArrays< localIndex > getEdgeToFaces() const override;
 
   ArrayOfArrays< localIndex > getFaceToNodes() const override;
 
   ArrayOfArrays< localIndex > getFaceToEdges() const override;
 
-  array2d< localIndex > getFaceToElements() const override;
+  ToCellRelation< array2d< localIndex > > getFaceToElements() const override;
 
   array1d< globalIndex > getNodeLocalToGlobal() const override;
 
@@ -195,14 +195,20 @@ class CellBlockManager : public CellBlockManagerABC
    */
   void buildFaceMaps();
 
+  template< typename BASEMAP, typename FUNC >
+  void buildToCellMap( localIndex const cellIndex,
+                       ToCellRelation< BASEMAP > & toCells,
+                       FUNC cellToObjectGetter,
+                       localIndex const overAlloc = 0 ) const;
+
   array2d< real64, nodes::REFERENCE_POSITION_PERM > m_nodesPositions;
 
-  ArrayOfSets< localIndex > m_nodeToEdges;
-  ArrayOfSets< localIndex > m_edgeToFaces;
+  ArrayOfArrays< localIndex > m_nodeToEdges;
+  ArrayOfArrays< localIndex > m_edgeToFaces;
   array2d< localIndex > m_edgeToNodes;
   ArrayOfArrays< localIndex >  m_faceToNodes;
   ArrayOfArrays< localIndex > m_faceToEdges;
-  array2d< localIndex >  m_faceToElements;
+  ToCellRelation< array2d< localIndex > > m_faceToCells;
 
   array1d< globalIndex >  m_nodeLocalToGlobal;
 
diff --git a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp
index 236e080afec..979d741bcc8 100644
--- a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp
+++ b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp
@@ -22,6 +22,17 @@
 namespace geosx
 {
 
+/**
+ * @brief Container for maps from a mesh object (node, edge or face) to cells.
+ * @tparam BASETYPE underlying map type
+ */
+template< typename BASETYPE >
+struct ToCellRelation
+{
+  BASETYPE toBlockIndex; ///< Map containing a list of cell block indices for each object
+  BASETYPE toCellIndex;  ///< Map containing cell indices, same shape as above
+};
+
 /**
  * @brief Abstract base class for CellBlockManager.
  */
@@ -120,28 +131,25 @@ class CellBlockManagerABC : public dataRepository::Group
    * @brief Returns the node coordinates in a (numNodes, 3) 2d array.
    * @return A const view to the array.
    */
-  virtual array2d< real64, nodes::REFERENCE_POSITION_PERM > getNodesPositions() const = 0;
+  virtual array2d< real64, nodes::REFERENCE_POSITION_PERM > getNodePositions() const = 0;
 
   /**
    * @brief Returns the node to edges mapping.
    * @return The one to many relationship.
    */
-  virtual ArrayOfSets< localIndex > getNodeToEdges() const = 0;
+  virtual ArrayOfArrays< localIndex > getNodeToEdges() const = 0;
 
   /**
    * @brief Returns the face to nodes mappings.
    * @return The one to many relationship.
    */
-  virtual ArrayOfSets< localIndex > getNodeToFaces() const = 0;
+  virtual ArrayOfArrays< localIndex > getNodeToFaces() const = 0;
 
   /**
    * @brief Returns the node to elements mapping.
    * @return A one to many relationship.
-   *
-   * @note The mapping is computed on the fly and returned. It is not stored in the instance.
    */
-  virtual ArrayOfArrays< localIndex > getNodeToElements() const = 0;
-
+  virtual ToCellRelation< ArrayOfArrays< localIndex > > getNodeToElements() const = 0;
   /**
    * @brief Returns the edge to nodes mapping.
    * @return A 1 to 2 relationship. The result is meant to have size (numEdges, 2).
@@ -152,7 +160,7 @@ class CellBlockManagerABC : public dataRepository::Group
    * @brief Returns the edge to faces mapping.
    * @return A one to many relationship.
    */
-  virtual ArrayOfSets< localIndex > getEdgeToFaces() const = 0;
+  virtual ArrayOfArrays< localIndex > getEdgeToFaces() const = 0;
 
   /**
    * @brief Returns the face to nodes mapping.
@@ -167,12 +175,12 @@ class CellBlockManagerABC : public dataRepository::Group
   virtual ArrayOfArrays< localIndex > getFaceToEdges() const = 0;
 
   /**
-   * @brief Returns the face to elements mappings.
+   * @brief Returns the face to elements mapping.
    * @return A 1 to 2 relationship. The result is meant to have size (numFaces, 2).
    *
    * In case the face only belongs to one single element, the second value of the table is -1.
    */
-  virtual array2d< localIndex > getFaceToElements() const = 0;
+  virtual ToCellRelation< array2d< localIndex > > getFaceToElements() const = 0;
 
   /**
    * @brief The node to global mapping for nodes.
diff --git a/src/coreComponents/mesh/generators/CellBlockUtilities.cpp b/src/coreComponents/mesh/generators/CellBlockUtilities.cpp
index 94b45a59489..e598fc2e821 100644
--- a/src/coreComponents/mesh/generators/CellBlockUtilities.cpp
+++ b/src/coreComponents/mesh/generators/CellBlockUtilities.cpp
@@ -397,7 +397,7 @@ createEdgesByLowestNode( localIndex const numNodes,
 void populateEdgeMaps( ArrayOfArraysView< EdgeBuilder const > const & edgesByLowestNode,
                        arrayView1d< localIndex const > const & uniqueEdgeOffsets,
                        ArrayOfArraysView< localIndex > const & faceToEdgeMap,
-                       ArrayOfSetsView< localIndex > const & edgeToFaceMap,
+                       ArrayOfArraysView< localIndex > const & edgeToFaceMap,
                        arrayView2d< localIndex > const & edgeToNodeMap )
 {
   localIndex const numNodes = edgesByLowestNode.size();
@@ -422,7 +422,7 @@ void populateEdgeMaps( ArrayOfArraysView< EdgeBuilder const > const & edgesByLow
       {
         EdgeBuilder const & edge = *first;
         faceToEdgeMap( edge.faceIndex, edge.edgeNumber ) = curEdgeID;
-        edgeToFaceMap.insertIntoSet( curEdgeID, edge.faceIndex );
+        edgeToFaceMap.emplaceBack( curEdgeID, edge.faceIndex );
         ++first;
       }
 
@@ -443,7 +443,7 @@ void resizeEdgeMaps( ArrayOfArraysView< EdgeBuilder const > const & edgesByLowes
                      arrayView1d< localIndex const > const & uniqueEdgeOffsets,
                      ArrayOfArraysView< localIndex const > const & faceToNodeMap,
                      ArrayOfArrays< localIndex > & faceToEdgeMap,
-                     ArrayOfSets< localIndex > & edgeToFaceMap,
+                     ArrayOfArrays< localIndex > & edgeToFaceMap,
                      array2d< localIndex > & edgeToNodeMap )
 {
   localIndex const numFaces = faceToNodeMap.size();
@@ -492,7 +492,7 @@ void resizeEdgeMaps( ArrayOfArraysView< EdgeBuilder const > const & edgesByLowes
 localIndex buildEdgeMaps( localIndex const numNodes,
                           ArrayOfArraysView< localIndex const > const & faceToNodeMap,
                           ArrayOfArrays< localIndex > & faceToEdgeMap,
-                          ArrayOfSets< localIndex > & edgeToFaceMap,
+                          ArrayOfArrays< localIndex > & edgeToFaceMap,
                           array2d< localIndex > & edgeToNodeMap )
 {
   GEOSX_MARK_FUNCTION;
diff --git a/src/coreComponents/mesh/generators/CellBlockUtilities.hpp b/src/coreComponents/mesh/generators/CellBlockUtilities.hpp
index b530efaae97..d519e59c415 100644
--- a/src/coreComponents/mesh/generators/CellBlockUtilities.hpp
+++ b/src/coreComponents/mesh/generators/CellBlockUtilities.hpp
@@ -35,7 +35,7 @@ namespace geosx
 localIndex buildEdgeMaps( localIndex numNodes,
                           ArrayOfArraysView< localIndex const > const & faceToNodeMap,
                           ArrayOfArrays< localIndex > & faceToEdgeMap,
-                          ArrayOfSets< localIndex > & edgeToFaceMap,
+                          ArrayOfArrays< localIndex > & edgeToFaceMap,
                           array2d< localIndex > & edgeToNodeMap );
 
 /**
diff --git a/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp b/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp
index 3b743e19ca2..12217696af2 100644
--- a/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp
+++ b/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp
@@ -731,7 +731,7 @@ void InternalMeshGenerator::generateMesh( DomainPartition & domain )
 
   cellBlockManager.setNumNodes( numNodes );
 
-  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > X = cellBlockManager.getNodesPositions();
+  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > X = cellBlockManager.getNodePositions();
 
   arrayView1d< globalIndex > const nodeLocalToGlobal = cellBlockManager.getNodeLocalToGlobal();
 
diff --git a/src/coreComponents/mesh/generators/PAMELAMeshGenerator.cpp b/src/coreComponents/mesh/generators/PAMELAMeshGenerator.cpp
index d74f319e66e..11b7ac1d73d 100644
--- a/src/coreComponents/mesh/generators/PAMELAMeshGenerator.cpp
+++ b/src/coreComponents/mesh/generators/PAMELAMeshGenerator.cpp
@@ -162,7 +162,7 @@ real64 importNodes( PAMELA::Mesh & srcMesh, // PAMELA is not const-correct,
                     CellBlockManager & cellBlockManager )
 {
   cellBlockManager.setNumNodes( srcMesh.get_PointCollection()->size_all() );
-  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > X = cellBlockManager.getNodesPositions();
+  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > X = cellBlockManager.getNodePositions();
 
   arrayView1d< globalIndex > const nodeLocalToGlobal = cellBlockManager.getNodeLocalToGlobal();
 
diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp
index 1e5b1072c69..697c33801e5 100644
--- a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp
+++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp
@@ -259,7 +259,7 @@ real64 writeMeshNodes( vtkUnstructuredGrid & mesh,
 
   // Writing the points
   arrayView1d< globalIndex > const & nodeLocalToGlobal = cellBlockManager.getNodeLocalToGlobal();
-  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > const & X = cellBlockManager.getNodesPositions();
+  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > const & X = cellBlockManager.getNodePositions();
 
   constexpr real64 minReal = LvArray::NumericLimits< real64 >::min;
   constexpr real64 maxReal = LvArray::NumericLimits< real64 >::max;
diff --git a/src/coreComponents/mesh/utilities/MeshMapUtilities.hpp b/src/coreComponents/mesh/utilities/MeshMapUtilities.hpp
index 64c2ed58b52..e357acba0d4 100644
--- a/src/coreComponents/mesh/utilities/MeshMapUtilities.hpp
+++ b/src/coreComponents/mesh/utilities/MeshMapUtilities.hpp
@@ -116,158 +116,97 @@ transposeIndexMap( VIEW_TYPE const & srcMap,
   return dstMap;
 }
 
-namespace internal
-{
-
-/// Element identifier containing (region index, subregion index, element index).
-using ElementTuple = std::tuple< localIndex, localIndex, localIndex >;
-
+/**
+ * @brief Convert ToCellRelation into ToElementRelation.
+ * @tparam POLICY execution policy
+ * @param blockToSubRegion a map from cell blocks to region/subregion pairs
+ * @param srcMap source map (object-to-cells)
+ * @param dstMap target map (object-to-elements)
+ */
 template< typename POLICY >
-void convertFromElementTupleMap( ArrayOfArraysView< ElementTuple const > const & elemList,
-                                 OrderedVariableToManyElementRelation & toElement )
+void transformCellBlockToRegionMap( arrayView2d< localIndex const > const & blockToSubRegion,
+                                    ToCellRelation< ArrayOfArrays< localIndex > > const & srcMap,
+                                    ToElementRelation< ArrayOfArrays< localIndex > > & dstMap )
 {
-  ArrayOfArrays< localIndex > & toElementRegionList = toElement.m_toElementRegion;
-  ArrayOfArrays< localIndex > & toElementSubRegionList = toElement.m_toElementSubRegion;
-  ArrayOfArrays< localIndex > & toElementList = toElement.m_toElementIndex;
-
-  localIndex const numObjects = elemList.size();
-
-  toElementRegionList.resizeFromOffsets( numObjects, elemList.toViewConst().getOffsets() );
-  toElementSubRegionList.resizeFromOffsets( numObjects, elemList.toViewConst().getOffsets() );
-  toElementList.resizeFromOffsets( numObjects, elemList.toViewConst().getOffsets() );
-
-  forAll< parallelHostPolicy >( numObjects, [elemList = elemList.toViewConst(),
-                                             elemRegion = toElementRegionList.toView(),
-                                             elemSubRegion = toElementSubRegionList.toView(),
-                                             elemIndex = toElementList.toView()]( localIndex const objIndex )
+  GEOSX_ASSERT_EQ( blockToSubRegion.size(), 2 );
+  localIndex const numObjects = srcMap.toCellIndex.size();
+
+  localIndex const * offsets = srcMap.toCellIndex.toViewConst().getOffsets();
+  dstMap.m_toElementRegion.resizeFromOffsets( numObjects, offsets );
+  dstMap.m_toElementSubRegion.resizeFromOffsets( numObjects, offsets );
+  dstMap.m_toElementIndex.resizeFromOffsets( numObjects, offsets );
+
+  forAll< POLICY >( numObjects, [toCell = srcMap.toCellIndex.toViewConst(),
+                                 toBlock = srcMap.toBlockIndex.toViewConst(),
+                                 toRegion = dstMap.m_toElementRegion.toView(),
+                                 toSubRegion = dstMap.m_toElementSubRegion.toView(),
+                                 toElement = dstMap.m_toElementIndex.toView(),
+                                 blockToSubRegion]( localIndex const objIndex )
   {
-    arraySlice1d< ElementTuple const > const elems = elemList[ objIndex ];
-    for( ElementTuple const & e : elems )
+    arraySlice1d< localIndex const > const cells = toCell[objIndex];
+    for( localIndex i = 0; i < cells.size(); ++i )
     {
-      elemRegion.emplaceBack( objIndex, std::get< 0 >( e ) );
-      elemSubRegion.emplaceBack( objIndex, std::get< 1 >( e ) );
-      elemIndex.emplaceBack( objIndex, std::get< 2 >( e ) );
-    }
-  } );
-}
-
-template< typename POLICY >
-void convertFromElementTupleMap( ArrayOfArraysView< ElementTuple const > const & elemList,
-                                 FixedToManyElementRelation & toElement )
-{
-  array2d< localIndex > & toElementRegionList = toElement.m_toElementRegion;
-  array2d< localIndex > & toElementSubRegionList = toElement.m_toElementSubRegion;
-  array2d< localIndex > & toElementList = toElement.m_toElementIndex;
-
-  localIndex const numObjects = elemList.size();
-  localIndex const maxNumElem = toElementList.size( 1 );
-
-  GEOSX_ERROR_IF_NE( toElementRegionList.size( 1 ), maxNumElem );
-  GEOSX_ERROR_IF_NE( toElementSubRegionList.size( 1 ), maxNumElem );
-
-  toElementRegionList.resizeDimension< 0 >( numObjects );
-  toElementSubRegionList.resizeDimension< 0 >( numObjects );
-  toElementList.resizeDimension< 0 >( numObjects );
-
-  // We allow a fixed-size map to represent a variable relationship, as long as
-  // the number of elements does not exceed the fixed size (set by the caller).
-  // In this case, a dummy value "-1" is used to represent non-present elements.
-  toElementRegionList.setValues< POLICY >( -1 );
-  toElementSubRegionList.setValues< POLICY >( -1 );
-  toElementList.setValues< POLICY >( -1 );
+      localIndex const blockIndex = toBlock( objIndex, i );
+      localIndex const er = blockToSubRegion( blockIndex, 0 );
+      localIndex const esr = blockToSubRegion( blockIndex, 1 );
 
-  forAll< parallelHostPolicy >( numObjects, [=, // needed to optionally capture maxNumElem in Debug
-                                             elemList = elemList.toViewConst(),
-                                             elemRegion = toElementRegionList.toView(),
-                                             elemSubRegion = toElementSubRegionList.toView(),
-                                             elemIndex = toElementList.toView()]( localIndex const objIndex )
-  {
-    arraySlice1d< ElementTuple const > const elems = elemList[ objIndex ];
-    GEOSX_ASSERT_GE( maxNumElem, elems.size() );
-    localIndex count = 0;
-    for( ElementTuple const & e : elems )
-    {
-      elemRegion( objIndex, count ) = std::get< 0 >( e );
-      elemSubRegion( objIndex, count ) = std::get< 1 >( e );
-      elemIndex( objIndex, count ) = std::get< 2 >( e );
-      ++count;
+      // Check needed because some blocks may remain unused
+      if( er >= 0 && esr >= 0 )
+      {
+        toRegion.emplaceBack( objIndex, er );
+        toSubRegion.emplaceBack( objIndex, esr );
+        toElement.emplaceBack( objIndex, cells[i] );
+      }
     }
   } );
 }
 
-}
-
 /**
- * @brief Build to-element map by inverting existing maps in element subregions.
- * @tparam BASEMAP underlying type of to-element map
- * @tparam FUNC type of @p cellToObjectGetter
- * @param elementRegionManager element region manager
- * @param numObjects number of objects (nodes, faces, etc.) for which maps are built
- * @param toElements container for to-element (region, subregion, index) maps
- * @param cellToObjectGetter function used to extract maps from subregions
- * @param overAlloc overallocation for the resulting maps
- *                  (extra capacity per row, only meaningful for variable-secondary-size containers)
+ * @brief Convert ToCellRelation into ToElementRelation.
+ * @tparam POLICY execution policy
+ * @param blockToSubRegion a map from cell blocks to region/subregion pairs
+ * @param srcMap source map (object-to-cells)
+ * @param dstMap target map (object-to-elements)
  */
-template< typename BASEMAP, typename FUNC >
-void buildToElementMaps( ElementRegionManager const & elementRegionManager,
-                         localIndex const numObjects,
-                         ToElementRelation< BASEMAP > & toElements,
-                         FUNC cellToObjectGetter,
-                         localIndex const overAlloc = 0 )
+template< typename POLICY, typename PERM >
+void transformCellBlockToRegionMap( arrayView2d< localIndex const > const & blockToSubRegion,
+                                    ToCellRelation< array2d< localIndex, PERM > > const & srcMap,
+                                    ToElementRelation< array2d< localIndex, PERM > > & dstMap )
 {
-  // Calculate the number of entries in each sub-array
-  array1d< localIndex > elemCounts( numObjects );
-  elemCounts.setValues< serialPolicy >( overAlloc );
-  elementRegionManager.forElementSubRegionsComplete< CellElementSubRegion >( [&]( localIndex const,
-                                                                                  localIndex const,
-                                                                                  ElementRegionBase const &,
-                                                                                  CellElementSubRegion const & subRegion )
+  GEOSX_ASSERT_EQ( blockToSubRegion.size(), 2 );
+  localIndex const numObjects = srcMap.toCellIndex.size( 0 );
+  localIndex const maxCellsPerObject = srcMap.toCellIndex.size( 1 );
+
+  dstMap.m_toElementRegion.resize( numObjects, maxCellsPerObject );
+  dstMap.m_toElementSubRegion.resize( numObjects, maxCellsPerObject );
+  dstMap.m_toElementIndex.resize( numObjects, maxCellsPerObject );
+
+  forAll< POLICY >( numObjects, [toCell = srcMap.toCellIndex.toViewConst(),
+                                 toBlock = srcMap.toBlockIndex.toViewConst(),
+                                 toRegion = dstMap.m_toElementRegion.toView(),
+                                 toSubRegion = dstMap.m_toElementSubRegion.toView(),
+                                 toElement = dstMap.m_toElementIndex.toView(),
+                                 blockToSubRegion,
+                                 maxCellsPerObject]( localIndex const objIndex )
   {
-    auto const elemToObject = cellToObjectGetter( subRegion );
-    forAll< parallelHostPolicy >( subRegion.size(), [elemCounts = elemCounts.toView(),
-                                                     elemToObject]( localIndex const ei )
+    arraySlice1d< localIndex const > const cells = toCell[objIndex];
+    localIndex cellCount = 0;
+    for( localIndex i = 0; i < maxCellsPerObject && cells[i] >= 0; ++i )
     {
-      auto const objects = elemToObject[ ei ];
-      // can't use range-based for loop when slice is not contiguous
-      for( localIndex i = 0; i < objects.size(); ++i )
-      {
-        RAJA::atomicInc< parallelHostAtomic >( &elemCounts[ objects[i] ] );
-      }
-    } );
-  } );
+      localIndex const blockIndex = toBlock( objIndex, i );
+      localIndex const er = blockToSubRegion( blockIndex, 0 );
+      localIndex const esr = blockToSubRegion( blockIndex, 1 );
 
-  // Allocate memory
-  ArrayOfArrays< internal::ElementTuple > elemList;
-  elemList.resizeFromCapacities< parallelHostPolicy >( numObjects, elemCounts.data() );
-
-  // Populate map of tuples in parallel
-  elementRegionManager.forElementSubRegionsComplete< CellElementSubRegion >( [&]( localIndex const er,
-                                                                                  localIndex const esr,
-                                                                                  ElementRegionBase const &,
-                                                                                  CellElementSubRegion const & subRegion )
-  {
-    auto const elemToObject = cellToObjectGetter( subRegion );
-    forAll< parallelHostPolicy >( subRegion.size(), [elemList = elemList.toView(),
-                                                     elemToObject, er, esr]( localIndex const ei )
-    {
-      auto const objects = elemToObject[ ei ];
-      // can't use range-based for loop when slice is not contiguous
-      for( localIndex i = 0; i < objects.size(); ++i )
+      // Check needed because some blocks may remain unused
+      if( er >= 0 && esr >= 0 )
       {
-        elemList.emplaceBackAtomic< parallelHostAtomic >( objects[i], er, esr, ei );
+        toRegion( objIndex, cellCount ) = er;
+        toSubRegion( objIndex, cellCount ) = esr;
+        toElement( objIndex, cellCount ) = cells[i];
+        ++cellCount;
       }
-    } );
-  } );
-
-  // Sort each element list to ensure unique race-condition-free map order
-  forAll< parallelHostPolicy >( numObjects, [elemList = elemList.toView()]( localIndex const objIndex )
-  {
-    arraySlice1d< internal::ElementTuple > const elems = elemList[ objIndex ];
-    LvArray::sortedArrayManipulation::makeSorted( elems.begin(), elems.end() );
+    }
   } );
-
-  // Finally, split arrays-of-tuples into separate arrays
-  internal::convertFromElementTupleMap< parallelHostPolicy >( elemList.toViewConst(), toElements );
 }
 
 } // namespace meshMapUtilities
diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
index 57647d9e6a5..8db330a6249 100644
--- a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
+++ b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
@@ -2097,7 +2097,7 @@ void SurfaceGenerator::performFracture( const localIndex nodeID,
                                                         getSubRegion< CellElementSubRegion >( faceToSubRegionMap[iFace][0] );
             arrayView2d< real64 const > const subRegionElemCenter = elementSubRegion.getElementCenter();
 
-            faceManager.sortFaceNodes( X, subRegionElemCenter[ elementIndex ], faceToNodeMap[ iFace ], faceToNodeMap.sizeOfArray( iFace ) );
+            FaceManager::sortFaceNodes( X, subRegionElemCenter[ elementIndex ], faceToNodeMap[ iFace ] );
 
             //Face normal need to be updated here
             real64 fCenter[ 3 ];
diff --git a/src/coreComponents/unitTests/meshTests/testPAMELAImport.cpp b/src/coreComponents/unitTests/meshTests/testPAMELAImport.cpp
index 1d330540dda..5236b603d55 100644
--- a/src/coreComponents/unitTests/meshTests/testPAMELAImport.cpp
+++ b/src/coreComponents/unitTests/meshTests/testPAMELAImport.cpp
@@ -65,7 +65,7 @@ void TestMeshImport( string const & inputStringMesh,
   elemManager.postProcessInputRecursive();
 
   CellBlockManager & cellBlockManager = meshBody.getGroup< CellBlockManager >( keys::cellManager );
-  nodeManager.setGeometricalRelations( cellBlockManager );
+  nodeManager.setGeometricalRelations( cellBlockManager, elemManager );
 
   elemManager.generateMesh( cellBlockManager );
 
diff --git a/src/coreComponents/unitTests/meshTests/testVTKImport.cpp b/src/coreComponents/unitTests/meshTests/testVTKImport.cpp
index eca9ed15eda..c53eeb9c7d1 100644
--- a/src/coreComponents/unitTests/meshTests/testVTKImport.cpp
+++ b/src/coreComponents/unitTests/meshTests/testVTKImport.cpp
@@ -92,7 +92,7 @@ TEST( VTKImport, cube )
 
     // The information in the tables is not filled yet. We can check the consistency of the sizes.
     ASSERT_EQ( cellBlockManager.getNodeToFaces().size(), expectedNumNodes );
-    ASSERT_EQ( cellBlockManager.getNodeToElements().size(), expectedNumNodes );
+    ASSERT_EQ( cellBlockManager.getNodeToElements().toCellIndex.size(), expectedNumNodes );
 
     // We have all the 4 x 4  x 4 = 64 nodes in the "all" set.
     SortedArray< localIndex > const & allNodes = cellBlockManager.getNodeSets().at( "all" );