feat: Implement EDFM loading from VTK

src/coreComponents/mainInterface/ProblemManager.cpp

@@ -685,6 +685,9 @@ void ProblemManager::generateMesh()
   domain.setupCommunications( useNonblockingMPI );
   domain.outputPartitionInformation();
 
+  // Create Embedded fractures here.
+  generateEmbeddedFractures();
+
   domain.forMeshBodies( [&]( MeshBody & meshBody )
   {
     if( meshBody.hasGroup( keys::particleManager ) )
@@ -728,6 +731,47 @@ void ProblemManager::generateMesh()
 
 }
 
+void ProblemManager::generateEmbeddedFractures()
+{
+  DomainPartition & domain = getDomainPartition();
+  MeshManager & meshManager = this->getGroup< MeshManager >( groupKeys.meshManager );
+
+  meshManager.generateEmbeddedFractures( domain );
+
+  MeshBody & meshBody = domain.getMeshBody( 0 );
+  CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager();
+  Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks();
+  string const & faceBlockName = embeddedSurfaceRegion.getFaceBlockName();
+
+  if( embSurfBlocks.hasGroup( faceBlockName ))
+  {
+    EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( faceBlockName );
+
+    elemManager.forElementSubRegionsComplete< CellElementSubRegion >( [&]( localIndex const er,
+                                                                           localIndex const esr,
+                                                                           ElementRegionBase &,
+                                                                           CellElementSubRegion & subRegion )
+    {
+      FixedOneToManyRelation const & cellToEdges = subRegion.edgeList();
+
+      embeddedSurfaceSubRegion.copyFromEmbeddedSurfaceBlock( er,
+                                                             esr,
+                                                             nodeManager,
+                                                             embSurfNodeManager,
+                                                             edgeManager,
+                                                             cellToEdges,
+                                                             embSurf );
+
+      // Add all the fracture information to the CellElementSubRegion
+      for( localIndex edfmIndex=0; edfmIndex < embSurf.numEmbeddedSurfElem(); ++edfmIndex )
+      {
+        localIndex cellIndex = embSurf.getEmbeddedSurfElemTo3dElem().toCellIndex[edfmIndex][0];
+        subRegion.addFracturedElement( cellIndex, edfmIndex );
+        newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( edfmIndex );
+      }
+    } );   // end loop over subregions
+  }
+}
 
 void ProblemManager::importFields()
 {

src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp

@@ -27,11 +27,6 @@
 #include "BufferOps.hpp"
 #include "mesh/MeshFields.hpp"
 
-
-#include "mainInterface/GeosxState.hpp"
-#include "mainInterface/ProblemManager.hpp"
-#include "mesh/DomainPartition.hpp"
-
 namespace geos
 {
 using namespace dataRepository;
@@ -345,14 +340,13 @@ array1d< localIndex > EmbeddedSurfaceSubRegion::getEdfmNodeParentEdgeIndex( Arra
   return parentEdgeIndex;
 }
 
-bool EmbeddedSurfaceSubRegion::addAllEmbeddedSurfaces(
-  localIndex const regionIndex,
-  localIndex const subRegionIndex,
-  NodeManager const & nodeManager,
-  EmbeddedSurfaceNodeManager & embSurfNodeManager,
-  EdgeManager const & edgeManager,
-  FixedOneToManyRelation const & cellToEdges,
-  EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock )
+void EmbeddedSurfaceSubRegion::copyFromEmbeddedSurfaceBlock( localIndex const regionIndex,
+                                                             localIndex const subRegionIndex,
+                                                             NodeManager const & nodeManager,
+                                                             EmbeddedSurfaceNodeManager & embSurfNodeManager,
+                                                             EdgeManager const & edgeManager,
+                                                             FixedOneToManyRelation const & cellToEdges,
+                                                             EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock )
 {
 
   arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoord = nodeManager.referencePosition();
@@ -419,9 +413,6 @@ bool EmbeddedSurfaceSubRegion::addAllEmbeddedSurfaces(
     LvArray::tensorOps::copy< 3 >( m_tangentVector2[i], elemTangentialVectors2[i] );
     this->calculateElementGeometricQuantities( elemNodeCoords.toViewConst(), i );
   }
-
-  return true;
-
 }
 
 void EmbeddedSurfaceSubRegion::inheritGhostRank( array1d< array1d< arrayView1d< integer const > > > const & cellGhostRank )

src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp

@@ -290,6 +290,18 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion
 
   ///@}
 
+  /**
+   * @brief Fill the EmbeddedSurfaceSubRegion by copying those of the source CellBlock
+   * @param embeddedSurfaceBlock the CellBlEmbeddedSurfaceBlock which properties (connectivity info) will be copied.
+   */
+  void copyFromEmbeddedSurfaceBlock( localIndex const regionIndex,
+                                     localIndex const subRegionIndex,
+                                     NodeManager const & nodeManager,
+                                     EmbeddedSurfaceNodeManager & embSurfNodeManager,
+                                     EdgeManager const & edgeManager,
+                                     FixedOneToManyRelation const & cellToEdges,
+                                     EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock );
+
 private:
 
   /**

src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp

@@ -75,7 +75,7 @@ EmbeddedSurfaceGenerator::EmbeddedSurfaceGenerator( const string & name,
 
   registerWrapper( viewKeyStruct::targetObjectsNameString(), &m_targetObjectsName ).
     setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ).
-    setInputFlag( dataRepository::InputFlags::OPTIONAL ).
+    setInputFlag( dataRepository::InputFlags::REQUIRED ).
     setDescription( "List of geometric objects that will be used to initialized the embedded surfaces/fractures." );
 
   registerWrapper( viewKeyStruct::mpiCommOrderString(), &m_mpiCommOrder ).
@@ -130,121 +130,81 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups()
 
   NewObjectLists newObjects;
 
-  if( m_targetObjectsName.empty() ) // when targetObjects were not provided - try to import from vtk
+  // Loop over all the fracture planes
+  geometricObjManager.forGeometricObject< PlanarGeometricObject >( m_targetObjectsName, [&]( localIndex const,
+                                                                                             PlanarGeometricObject & fracture )
   {
-    MeshBody & meshBody = domain.getMeshBody( 0 );
-    CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager();
-    Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks();
-    string const & faceBlockName = embeddedSurfaceRegion.getFaceBlockName();
-    if( embSurfBlocks.hasGroup( faceBlockName ))
+    /* 1. Find out if an element is cut by the fracture or not.
+     * Loop over all the elements and for each one of them loop over the nodes and compute the
+     * dot product between the distance between the plane center and the node and the normal
+     * vector defining the plane. If two scalar products have different signs the plane cuts the
+     * cell. If a nodes gives a 0 dot product it has to be neglected or the method won't work.
+     */
+    real64 const planeCenter[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( fracture.getCenter() );
+    real64 const normalVector[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( fracture.getNormal() );
+    // Initialize variables
+    globalIndex nodeIndex;
+    integer isPositive, isNegative;
+    real64 distVec[ 3 ];
+
+    elemManager.forElementSubRegionsComplete< CellElementSubRegion >(
+      [&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion )
     {
-      EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( faceBlockName );
-
-      elemManager.forElementSubRegionsComplete< CellElementSubRegion >(
-        [&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion )
-      {
-        FixedOneToManyRelation const & cellToEdges = subRegion.edgeList();
+      arrayView2d< localIndex const, cells::NODE_MAP_USD > const cellToNodes = subRegion.nodeList();
+      FixedOneToManyRelation const & cellToEdges = subRegion.edgeList();
 
-        bool added = embeddedSurfaceSubRegion.addAllEmbeddedSurfaces( er,
-                                                                      esr,
-                                                                      nodeManager,
-                                                                      embSurfNodeManager,
-                                                                      edgeManager,
-                                                                      cellToEdges,
-                                                                      embSurf );
+      arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
 
-        if( added )
-        {
-          // Add all the fracture information to the CellElementSubRegion
-          for( localIndex edfmIndex=0; edfmIndex < embSurf.numEmbeddedSurfElem(); ++edfmIndex )
-          {
-            localIndex cellIndex = embSurf.getEmbeddedSurfElemTo3dElem().toCellIndex[edfmIndex][0];
-            subRegion.addFracturedElement( cellIndex, edfmIndex );
-            newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( edfmIndex );
-          }
-        }
-      } ); // end loop over subregions
-    }
-  }
-  else
-  {
-    // Loop over all the fracture planes
-    geometricObjManager.forGeometricObject< PlanarGeometricObject >( m_targetObjectsName, [&]( localIndex const,
-                                                                                               PlanarGeometricObject & fracture )
-    {
-      /* 1. Find out if an element is cut by the fracture or not.
-       * Loop over all the elements and for each one of them loop over the nodes and compute the
-       * dot product between the distance between the plane center and the node and the normal
-       * vector defining the plane. If two scalar products have different signs the plane cuts the
-       * cell. If a nodes gives a 0 dot product it has to be neglected or the method won't work.
-       */
-      real64 const planeCenter[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( fracture.getCenter() );
-      real64 const normalVector[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( fracture.getNormal() );
-      // Initialize variables
-      globalIndex nodeIndex;
-      integer isPositive, isNegative;
-      real64 distVec[ 3 ];
-
-      elemManager.forElementSubRegionsComplete< CellElementSubRegion >(
-        [&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion )
+      forAll< serialPolicy >( subRegion.size(), [ &, ghostRank,
+                                                  cellToNodes,
+                                                  nodesCoord ] ( localIndex const cellIndex )
       {
-        arrayView2d< localIndex const, cells::NODE_MAP_USD > const cellToNodes = subRegion.nodeList();
-        FixedOneToManyRelation const & cellToEdges = subRegion.edgeList();
-
-        arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
-
-        forAll< serialPolicy >( subRegion.size(), [ &, ghostRank,
-                                                    cellToNodes,
-                                                    nodesCoord ] ( localIndex const cellIndex )
+        if( ghostRank[cellIndex] < 0 )
         {
-          if( ghostRank[cellIndex] < 0 )
+          isPositive = 0;
+          isNegative = 0;
+          for( localIndex kn = 0; kn < subRegion.numNodesPerElement(); kn++ )
           {
-            isPositive = 0;
-            isNegative = 0;
-            for( localIndex kn = 0; kn < subRegion.numNodesPerElement(); kn++ )
+            nodeIndex = cellToNodes[cellIndex][kn];
+            LvArray::tensorOps::copy< 3 >( distVec, nodesCoord[nodeIndex] );
+            LvArray::tensorOps::subtract< 3 >( distVec, planeCenter );
+            // check if the dot product is zero
+            if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) > 0 )
             {
-              nodeIndex = cellToNodes[cellIndex][kn];
-              LvArray::tensorOps::copy< 3 >( distVec, nodesCoord[nodeIndex] );
-              LvArray::tensorOps::subtract< 3 >( distVec, planeCenter );
-              // check if the dot product is zero
-              if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) > 0 )
-              {
-                isPositive = 1;
-              }
-              else if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) < 0 )
-              {
-                isNegative = 1;
-              }
-            } // end loop over nodes
-
-            if( isPositive * isNegative == 1 )
+              isPositive = 1;
+            }
+            else if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) < 0 )
             {
-              bool added = embeddedSurfaceSubRegion.addNewEmbeddedSurface( cellIndex,
-                                                                           er,
-                                                                           esr,
-                                                                           nodeManager,
-                                                                           embSurfNodeManager,
-                                                                           edgeManager,
-                                                                           cellToEdges,
-                                                                           &fracture );
-
-              if( added )
-              {
-                GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator, "Element " << cellIndex << " is fractured" );
-
-                // Add the information to the CellElementSubRegion
-                subRegion.addFracturedElement( cellIndex, localNumberOfSurfaceElems );
-
-                newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( localNumberOfSurfaceElems );
-
-                localNumberOfSurfaceElems++;
-              }
+              isNegative = 1;
+            }
+          } // end loop over nodes
+          if( isPositive * isNegative == 1 )
+          {
+            bool added = embeddedSurfaceSubRegion.addNewEmbeddedSurface( cellIndex,
+                                                                         er,
+                                                                         esr,
+                                                                         nodeManager,
+                                                                         embSurfNodeManager,
+                                                                         edgeManager,
+                                                                         cellToEdges,
+                                                                         &fracture );
+
+            if( added )
+            {
+              GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator, "Element " << cellIndex << " is fractured" );
+
+              // Add the information to the CellElementSubRegion
+              subRegion.addFracturedElement( cellIndex, localNumberOfSurfaceElems );
+
+              newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( localNumberOfSurfaceElems );
+
+              localNumberOfSurfaceElems++;
             }
           }
-        } );// end loop over cells
-      } );// end loop over subregions
-    } );// end loop over planes
-  }
+        }
+      } );// end loop over cells
+    } );// end loop over subregions
+  } );// end loop over planes
 
   // Launch kernel to compute connectivity index of each fractured element.
   elemManager.forElementSubRegionsComplete< CellElementSubRegion >(