feat: Add functions to connect well perforation to surface elements

inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml

@@ -29,7 +29,7 @@
 
     <ConstantPermeability
       name="rockPerm"
-      permeabilityComponents="{ 1.0e-18, 1.0e-18, 1.0e-18 }" />
+      permeabilityComponents="{ 1.0e-14, 1.0e-14, 1.0e-14 }" />
       <!-- Material inside the fault -->
     <CompressibleSolidParallelPlatesPermeability
       name="faultFilling"
@@ -188,14 +188,15 @@
         logLevel="1"
         wellRegionName="wellRegion2"
         wellControlsName="wellControls2"
-        polylineNodeCoords="{ { -5400, -5400, 0 },
-                              { -5400, -5400, -2500 } }"
+        polylineNodeCoords="{ { 1200, -700, 0 },
+                              { 1200, -700, -2500 } }"
         polylineSegmentConn="{ { 0, 1 } }"
         radius="0.1"
         numElementsPerSegment="1">
         <Perforation
           name="injector1_perf1"
-          distanceFromHead="2500"/>
+          targetRegion="Fault"
+          distanceFromHead="1700"/>
       </InternalWell>
     </VTKMesh>
   </Mesh>

inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml

@@ -39,6 +39,7 @@
       targetRegions="{ Region, Fault, wellRegion1, wellRegion2 }">
       <NonlinearSolverParameters
         newtonMaxIter="40"
+        newtonTol="1e-4"
         maxAllowedResidualNorm="1e+25"/>
       <LinearSolverParameters
         directParallel="0"/>

src/coreComponents/mesh/DomainPartition.cpp

@@ -242,10 +242,11 @@ void DomainPartition::setupCommunications( bool use_nonblocking )
         {
           NodeManager & nodeManager = meshLevel.getNodeManager();
           FaceManager & faceManager = meshLevel.getFaceManager();
+          ElementRegionManager & elemManager = meshLevel.getElemManager();
 
           CommunicationTools::getInstance().setupGhosts( meshLevel, m_neighbors, use_nonblocking );
-          faceManager.sortAllFaceNodes( nodeManager, meshLevel.getElemManager() );
-          faceManager.computeGeometry( nodeManager );
+          faceManager.sortAllFaceNodes( nodeManager, elemManager );
+          faceManager.computeGeometry( nodeManager, elemManager );
         }
         else if( !meshLevel.isShallowCopyOf( meshBody.getMeshLevels().getGroup< MeshLevel >( 0 )) )
         {

src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp

@@ -47,14 +47,6 @@ EmbeddedSurfaceSubRegion::EmbeddedSurfaceSubRegion( string const & name,
 {
   m_elementType = ElementType::Polygon;
 
-  registerWrapper( viewKeyStruct::elementCenterString(), &m_elementCenter ).
-    setDescription( "The center of each EmbeddedSurface element." );
-
-  registerWrapper( viewKeyStruct::elementVolumeString(), &m_elementVolume ).
-    setApplyDefaultValue( -1.0 ).
-    setPlotLevel( dataRepository::PlotLevel::LEVEL_0 ).
-    setDescription( "The volume of each EmbeddedSurface element." );
-
   registerWrapper( viewKeyStruct::connectivityIndexString(), &m_connectivityIndex ).
     setApplyDefaultValue( 1 ).
     setDescription( "Connectivity index of each EmbeddedSurface." );

src/coreComponents/mesh/FaceManager.cpp

@@ -203,7 +203,7 @@ void FaceManager::setGeometricalRelations( CellBlockManagerABC const & cellBlock
 
   if( isBaseMeshLevel )
   {
-    computeGeometry( nodeManager );
+    computeGeometry( nodeManager, elemRegionManager );
   }
 }
 
@@ -217,7 +217,8 @@ void FaceManager::setupRelatedObjectsInRelations( NodeManager const & nodeManage
   m_toElements.setElementRegionManager( elemRegionManager );
 }
 
-void FaceManager::computeGeometry( NodeManager const & nodeManager )
+void FaceManager::computeGeometry( NodeManager const & nodeManager,
+                                   ElementRegionManager const & elemManager )
 {
   arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X = nodeManager.referencePosition();
 
@@ -230,6 +231,11 @@ void FaceManager::computeGeometry( NodeManager const & nodeManager )
                                                                          m_faceNormal[ faceIndex ] );
 
   } );
+
+  elemManager.forElementSubRegions< CellElementSubRegion, FaceElementSubRegion >( [&] ( auto const & subRegion )
+  {
+    subRegion.calculateElementCenters( X );
+  } );
 }
 
 void FaceManager::setIsExternal()
@@ -258,11 +264,6 @@ void FaceManager::sortAllFaceNodes( NodeManager const & nodeManager,
 
   ArrayOfArraysView< localIndex > const facesToNodes = nodeList().toView();
 
-  elemManager.forElementSubRegions< CellElementSubRegion, FaceElementSubRegion >( [&] ( auto const & subRegion )
-  {
-    subRegion.calculateElementCenters( X );
-  } );
-
   ElementRegionManager::ElementViewAccessor< arrayView2d< real64 const > > elemCenter =
     elemManager.constructArrayViewAccessor< real64, 2 >( ElementSubRegionBase::viewKeyStruct::elementCenterString() );
 

src/coreComponents/mesh/FaceManager.hpp

@@ -137,8 +137,10 @@ class FaceManager : public ObjectManagerBase
   /**
    * @brief Compute faces center, area and normal.
    * @param[in] nodeManager NodeManager associated with the current DomainPartition
+   * @param[in] elemManager element manager allowing access to the cell elements
    */
-  void computeGeometry( NodeManager const & nodeManager );
+  void computeGeometry( NodeManager const & nodeManager,
+                        ElementRegionManager const & elemManager );
 
   /**
    * @brief Builds the face-on-domain-boundary indicator.

src/coreComponents/mesh/PerforationData.cpp

@@ -233,23 +233,40 @@ void PerforationData::getReservoirElementDimensions( MeshLevel const & mesh,
 {
   ElementRegionManager const & elemManager = mesh.getElemManager();
   NodeManager const & nodeManager = mesh.getNodeManager();
-  CellElementRegion const & region = elemManager.getRegion< CellElementRegion >( er );
-  CellElementSubRegion const & subRegion = region.getSubRegion< CellElementSubRegion >( esr );
-
-  // compute the bounding box of the element
-  real64 boxDims[ 3 ];
-  computationalGeometry::getBoundingBox( ei,
-                                         subRegion.nodeList(),
-                                         nodeManager.referencePosition(),
-                                         boxDims );
-
-  // dx and dz from bounding box
-  dx = boxDims[ 0 ];
-  dy = boxDims[ 1 ];
-
-  // dz is computed as vol / (dx * dy)
-  dz  = subRegion.getElementVolume()[ei];
-  dz /= dx * dy;
+  ElementRegionBase const & region = elemManager.getRegion< ElementRegionBase >( er );
+  ElementSubRegionBase const & subRegionBase = region.getSubRegion< ElementSubRegionBase >( esr );
+
+  region.applyLambdaToContainer< CellElementSubRegion, SurfaceElementSubRegion >( subRegionBase, [&]( auto const & subRegion )
+  {
+    // compute the bounding box of the element
+    real64 boxDims[ 3 ];
+    computationalGeometry::getBoundingBox( ei,
+                                           subRegion.nodeList(),
+                                           nodeManager.referencePosition(),
+                                           boxDims );
+
+    // dx and dz from bounding box
+    dx = boxDims[ 0 ];
+    dy = boxDims[ 1 ];
+    dz = boxDims[ 2 ];
+
+    if( dx < 1e-10 )
+    {
+      dx  = subRegion.getElementVolume()[ei];
+      dx /= dy * dz;
+    }
+    else if( dy < 1e-10 )
+    {
+      dy  = subRegion.getElementVolume()[ei];
+      dy /= dx * dz;
+    }
+    else
+    {
+      // dz is computed as vol / (dx * dy)
+      dz  = subRegion.getElementVolume()[ei];
+      dz /= dx * dy;
+    }
+  } );
 }
 
 void PerforationData::connectToWellElements( LineBlockABC const & lineBlock,

src/coreComponents/mesh/WellElementSubRegion.cpp

@@ -21,6 +21,7 @@
 #include "common/MpiWrapper.hpp"
 #include "LvArray/src/output.hpp"
 
+#include <unordered_set>
 
 namespace geos
 {
@@ -143,7 +144,6 @@ bool isPointInsideElement( SUBREGION_TYPE const & GEOS_UNUSED_PARAM( subRegion )
                            arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & GEOS_UNUSED_PARAM( referencePosition ),
                            localIndex const & GEOS_UNUSED_PARAM( eiLocal ),
                            ArrayOfArraysView< localIndex const > const & GEOS_UNUSED_PARAM( facesToNodes ),
-                           real64 const (&GEOS_UNUSED_PARAM( elemCenter ))[3],
                            real64 const (&GEOS_UNUSED_PARAM( location ))[3] )
 {
   // only CellElementSubRegion is currently supported
@@ -154,17 +154,72 @@ bool isPointInsideElement( CellElementSubRegion const & subRegion,
                            arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & referencePosition,
                            localIndex const & eiLocal,
                            ArrayOfArraysView< localIndex const > const & facesToNodes,
-                           real64 const (&elemCenter)[3],
                            real64 const (&location)[3] )
 {
   arrayView2d< localIndex const > const elemsToFaces = subRegion.faceList();
+  arrayView2d< real64 const > const elemCenters = subRegion.getElementCenter();
+  real64 const elemCenter[3] = { elemCenters[eiLocal][0],
+                                 elemCenters[eiLocal][1],
+                                 elemCenters[eiLocal][2] };
   return computationalGeometry::isPointInsidePolyhedron( referencePosition,
                                                          elemsToFaces[eiLocal],
                                                          facesToNodes,
                                                          elemCenter,
                                                          location );
 }
 
+// Define a hash function
+template< typename POINT_TYPE >
+struct PointHash
+{
+  std::size_t operator()( POINT_TYPE const point ) const
+  {
+    std::size_t h1 = std::hash< real64 >()( point[0] );
+    std::size_t h2 = std::hash< real64 >()( point[1] );
+    std::size_t h3 = std::hash< real64 >()( point[2] );
+    return h1 ^ (h2 << 1) ^ (h3 << 2);
+  }
+};
+
+// Define equality operator
+template< typename POINT_TYPE >
+struct PointsEqual
+{
+  bool operator()( POINT_TYPE const & p1, POINT_TYPE const & p2 ) const
+  {
+    return (std::abs( p1[0] - p2[0] ) < 1e-10) && (std::abs( p1[1] - p2[1] ) < 1e-10) && (std::abs( p1[2] - p2[2] ) < 1e-10);
+  }
+};
+
+bool isPointInsideElement( SurfaceElementSubRegion const & subRegion,
+                           arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & referencePosition,
+                           localIndex const & eiLocal,
+                           ArrayOfArraysView< localIndex const > const & GEOS_UNUSED_PARAM( facesToNodes ),
+                           real64 const (&location)[3] )
+{
+  typedef std::array< real64, 3 > Point3d;
+
+  // collect element nodes
+  integer const nV = subRegion.numNodesPerElement( eiLocal );
+  SurfaceElementSubRegion::NodeMapType const & nodeList = subRegion.nodeList();
+  std::vector< Point3d > polygon( nV );
+  for( integer i = 0; i < nV; ++i )
+  {
+    for( integer j = 0; j < 3; ++j )
+    {
+      polygon[i][j] = referencePosition[nodeList( eiLocal, i )][j];
+    }
+  }
+
+  // remove duplicates
+  std::unordered_set< Point3d, PointHash< Point3d >, PointsEqual< Point3d > >
+  unique_points( polygon.begin(), polygon.end());
+  polygon.clear();
+  std::copy( unique_points.begin(), unique_points.end(), std::back_inserter( polygon ));
+
+  return computationalGeometry::isPointInPolygon3d( polygon, polygon.size(), location );
+}
+
 /**
  * @brief Search the reservoir elements that can be accessed from the set "nodes".
           Stop if a reservoir element containing the perforation is found.
@@ -336,22 +391,16 @@ bool visitNeighborElements( MeshLevel const & mesh,
 
       ElementRegionBase const & region = elemManager.getRegion< ElementRegionBase >( er );
       SUBREGION_TYPE const & subRegion = region.getSubRegion< SUBREGION_TYPE >( esr );
-      arrayView2d< real64 const > const elemCenters = subRegion.getElementCenter();
-
       globalIndex const eiGlobal = subRegion.localToGlobalMap()[eiLocal];
 
       // if this element has not been visited yet, save it
       if( !elements.contains( eiGlobal ))
       {
         elements.insert( eiGlobal );
 
-        real64 const elemCenter[3] = { elemCenters[eiLocal][0],
-                                       elemCenters[eiLocal][1],
-                                       elemCenters[eiLocal][2] };
-
         // perform the test to see if the point is in this reservoir element
         // if the point is in the resevoir element, save the indices and stop the search
-        if( isPointInsideElement( subRegion, referencePosition, eiLocal, facesToNodes, elemCenter, location ) )
+        if( isPointInsideElement( subRegion, referencePosition, eiLocal, facesToNodes, location ) )
         {
           eiMatched = eiLocal;
           giMatched = eiGlobal;
@@ -432,6 +481,7 @@ void initializeLocalSearch( MeshLevel const & mesh,
                             localIndex & eiInit )
 {
   ElementSubRegionBase const & subRegion = mesh.getElemManager().getRegion( targetRegionIndex ).getSubRegion( targetSubRegionIndex );
+
   ElementRegionManager::ElementViewAccessor< arrayView2d< real64 const > >
   resElemCenter = mesh.getElemManager().constructViewAccessor< array2d< real64 >,
                                                                arrayView2d< real64 const > >( ElementSubRegionBase::viewKeyStruct::elementCenterString() );