diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_new_smoke.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_new_smoke.xml
new file mode 100644
index 0000000000..05679f4825
--- /dev/null
+++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_new_smoke.xml
@@ -0,0 +1,116 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml
new file mode 100644
index 0000000000..f91d3f1e8f
--- /dev/null
+++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml
@@ -0,0 +1,235 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml
new file mode 100644
index 0000000000..2f56fcd462
--- /dev/null
+++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml
@@ -0,0 +1,120 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/coreComponents/common/format/table/TableFormatter.cpp b/src/coreComponents/common/format/table/TableFormatter.cpp
index 7392834e06..3e6a032f8e 100644
--- a/src/coreComponents/common/format/table/TableFormatter.cpp
+++ b/src/coreComponents/common/format/table/TableFormatter.cpp
@@ -172,7 +172,6 @@ string TableTextFormatter::toString< TableData >( TableData const & tableData )
outputLayout( tableOutput, columns, msgTableError, sectionSeparatingLine );
outputSectionRows( columns, sectionSeparatingLine, tableOutput, nbValuesRows, TableLayout::Section::values );
- tableOutput << '\n';
return tableOutput.str();
}
diff --git a/src/coreComponents/common/format/table/unitTests/testTable.cpp b/src/coreComponents/common/format/table/unitTests/testTable.cpp
index 0258385b1a..893eb03ef7 100644
--- a/src/coreComponents/common/format/table/unitTests/testTable.cpp
+++ b/src/coreComponents/common/format/table/unitTests/testTable.cpp
@@ -55,7 +55,7 @@ TEST( testTable, tableEmptyRow )
"| value1 | [30.21543] | 3.0 | 54 | 0 |\n"
"| | | | | |\n"
"| Duis fringilla, ligula sed porta fringilla, ligula wisi commodo felis,ut adipiscing felis dui in enim. Suspendisse malesuada ultrices ante | [30.21543] | 30.45465142 | 787442 | 10 |\n"
- "-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n\n"
+ "-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n"
);
}
@@ -83,7 +83,7 @@ TEST( testTable, tableClassic )
"| value1 | [30.21543] | 3.0 | 54 | 0 |\n"
"| | | | | |\n"
"| value23 | [30.21543] | 30.45465142 | 787442 | 10 |\n"
- "-----------------------------------------------------------------------------------------------------------------------------------------------------------\n\n"
+ "-----------------------------------------------------------------------------------------------------------------------------------------------------------\n"
);
}
@@ -112,7 +112,7 @@ TEST( testTable, tableColumnParamClassic )
"-------------------------------------------------------------------------------------------\n"
"| value1 | | 3.0 | 3.0129877 | 2 | 1 |\n"
"| val1 | v | [3.045,42.02,89.25] | 3 | 10 | 3 |\n"
- "-------------------------------------------------------------------------------------------\n\n"
+ "-------------------------------------------------------------------------------------------\n"
);
}
@@ -140,7 +140,7 @@ TEST( testTable, tableHiddenColumn )
"----------------------------------------------------------------------------------------------------------------\n"
"| value1 | | 3.0 | 3.0129877 |\n"
"| val1 | v | [3.045,42.02,89.25] | 3 |\n"
- "----------------------------------------------------------------------------------------------------------------\n\n" );
+ "----------------------------------------------------------------------------------------------------------------\n" );
}
TEST( testTable, tableUniqueColumn )
@@ -162,7 +162,7 @@ TEST( testTable, tableUniqueColumn )
"---------------------------------------------------------------------------------------------------------------\n"
"| value1 |\n"
"| val1 |\n"
- "---------------------------------------------------------------------------------------------------------------\n\n" );
+ "---------------------------------------------------------------------------------------------------------------\n" );
}
TEST( testTable, tableEmptyTitle )
@@ -189,7 +189,7 @@ TEST( testTable, tableEmptyTitle )
"-------------------------------------------------------------------------------------------\n"
"| value1 | | 3.0 | 3.0129877 | 2 | 1 |\n"
"| val1 | v | [3.045,42.02,89.25] | 3 | 10 | 3 |\n"
- "-------------------------------------------------------------------------------------------\n\n"
+ "-------------------------------------------------------------------------------------------\n"
);
}
@@ -226,7 +226,7 @@ TEST( testTable, table2DTable )
"| Temperature = 300 | 0.03 | 0.02 |\n"
"| Temperature = 350 | 0.035 | 0.023333333333333334 |\n"
"| Temperature = 400 | 0.04 | 0.02666666666666667 |\n"
- "---------------------------------------------------------------------\n\n"
+ "---------------------------------------------------------------------\n"
);
}
@@ -266,7 +266,7 @@ TEST( testTable, table2DColumnMismatch )
"| Temperature = 300 | 0.03 | 0.02 |\n"
"| Temperature = 350 | 0.035 | |\n"
"| Temperature = 400 | 0.04 | 0.02666666666666667 |\n"
- "--------------------------------------------------------------------\n\n"
+ "--------------------------------------------------------------------\n"
);
}
}
@@ -321,7 +321,7 @@ TEST( testTable, tableSetMargin )
// "------------------------------------------------------------\n"
// "| value 1 |long value 1| 3.0034|3.0129877| | 1|\n"
// "| value 1 |long value 2| 100.45|4.0129877| 1| 2|\n"
-// "------------------------------------------------------------\n\n"
+// "------------------------------------------------------------\n"
// );
// }
}
diff --git a/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp b/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp
index 500d1da5ec..968bbfe3c2 100644
--- a/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp
+++ b/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp
@@ -22,7 +22,6 @@
#include "finiteVolume/FluxApproximationBase.hpp"
#include "finiteVolume/HybridMimeticDiscretization.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
namespace geos
diff --git a/src/coreComponents/mesh/MeshForLoopInterface.hpp b/src/coreComponents/mesh/MeshForLoopInterface.hpp
index b7c126b1e2..eb1e6040d8 100644
--- a/src/coreComponents/mesh/MeshForLoopInterface.hpp
+++ b/src/coreComponents/mesh/MeshForLoopInterface.hpp
@@ -69,11 +69,11 @@ minLocOverElemsInMesh( MeshLevel const & mesh, LAMBDA && lambda )
ElementRegionManager const & elemManager = mesh.getElemManager();
- for( localIndex er=0; er( [&]( localIndex const esr, CellElementSubRegion const & subRegion )
+ elemRegion.forElementSubRegionsIndex< ElementSubRegionBase >( [&]( localIndex const esr, ElementSubRegionBase const & subRegion )
{
localIndex const size = subRegion.size();
for( localIndex k = 0; k < size; ++k )
@@ -93,6 +93,68 @@ minLocOverElemsInMesh( MeshLevel const & mesh, LAMBDA && lambda )
return std::make_pair( minVal, std::make_tuple( minReg, minSubreg, minIndex ));
}
+/**
+ * @brief @return Return the minimum location/indices for a value condition specified by @p lambda.
+ * @tparam LAMBDA The type of the lambda function to be used to specify the minimum condition.
+ * @param region The region that will have all of its elements processed by this function.
+ * @param lambda A lambda function that returns as value that will be used in the minimum comparison.
+ */
+template< typename LAMBDA >
+auto
+minLocOverElemsInRegion( ElementRegionBase const & region, LAMBDA && lambda )
+{
+ using NUMBER = decltype( lambda( 0, 0 ) );
+
+ NUMBER minVal = std::numeric_limits< NUMBER >::max();
+ localIndex minSubreg = -1, minIndex = -1;
+
+ region.forElementSubRegionsIndex< ElementSubRegionBase >( [&]( localIndex const esr, ElementSubRegionBase const & subRegion )
+ {
+ localIndex const size = subRegion.size();
+ for( localIndex k = 0; k < size; ++k )
+ {
+ NUMBER const val = lambda( esr, k );
+ if( val < minVal )
+ {
+ minVal = val;
+ minSubreg = esr;
+ minIndex = k;
+ }
+ }
+ } );
+
+ return std::make_pair( minVal, std::make_tuple( minSubreg, minIndex ));
+}
+
+/**
+ * @brief @return Return the minimum location/indices for a value condition specified by @p lambda.
+ * @tparam LAMBDA The type of the lambda function to be used to specify the minimum condition.
+ * @param subRegion The subregion that will have all of its elements processed by this function.
+ * @param lambda A lambda function that returns as value that will be used in the minimum comparison.
+ */
+template< typename LAMBDA >
+auto
+minLocOverElemsInSubRegion( ElementSubRegionBase const & subRegion, LAMBDA && lambda )
+{
+ using NUMBER = decltype( lambda( 0 ) );
+
+ NUMBER minVal = std::numeric_limits< NUMBER >::max();
+ localIndex minIndex = -1;
+
+ localIndex const size = subRegion.size();
+ for( localIndex k = 0; k < size; ++k )
+ {
+ NUMBER const val = lambda( k );
+ if( val < minVal )
+ {
+ minVal = val;
+ minIndex = k;
+ }
+ }
+
+ return std::make_pair( minVal, minIndex );
+}
+
} // namespace geos
#endif // GEOS_MESH_MESHFORLOOPINTERFACE_HPP
diff --git a/src/coreComponents/mesh/Perforation.cpp b/src/coreComponents/mesh/Perforation.cpp
index 74ade6d426..77c5ea57a8 100644
--- a/src/coreComponents/mesh/Perforation.cpp
+++ b/src/coreComponents/mesh/Perforation.cpp
@@ -46,6 +46,11 @@ Perforation::Perforation( string const & name, Group * const parent )
setApplyDefaultValue( 0.0 ).
setInputFlag( InputFlags::OPTIONAL ).
setDescription( "Perforation skin factor" );
+
+ registerWrapper( viewKeyStruct::targetRegionString(), &m_targetRegionName ).
+ setRTTypeName( rtTypes::CustomTypes::groupNameRef ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Target region to connect the perforation" );
}
diff --git a/src/coreComponents/mesh/Perforation.hpp b/src/coreComponents/mesh/Perforation.hpp
index 5a3a3ae872..fb3636ef5d 100644
--- a/src/coreComponents/mesh/Perforation.hpp
+++ b/src/coreComponents/mesh/Perforation.hpp
@@ -102,6 +102,12 @@ class Perforation : public dataRepository::Group
*/
real64 getWellSkinFactor() const { return m_wellSkinFactor; }
+ /**
+ * @brief Get the target region for the perforation.
+ * @return the list of target regions
+ */
+ string const & getTargetRegion() const { return m_targetRegionName; }
+
///@}
/**
@@ -116,12 +122,8 @@ class Perforation : public dataRepository::Group
static constexpr char const *wellTransmissibilityString() { return "transmissibility"; }
/// @return String key for the well skin factor at this perforation
static constexpr char const *wellSkinFactorString() { return "skinFactor"; }
- /// ViewKey for the linear distance from well head
- dataRepository::ViewKey distanceFromHead = {distanceFromHeadString()};
- /// ViewKey for the well transmissibility at this perforation
- dataRepository::ViewKey wellTransmissibility = {wellTransmissibilityString()};
- /// ViewKey for the well transmissibility at this perforation
- dataRepository::ViewKey wellSkinFactor = { wellSkinFactorString() };
+ /// @return Target region for this perforation
+ static constexpr char const *targetRegionString() { return "targetRegion"; }
}
/// ViewKey struct for the Perforation class
viewKeysPerforation;
@@ -138,6 +140,9 @@ class Perforation : public dataRepository::Group
/// Well skin factor at this perforation
real64 m_wellSkinFactor;
+
+ /// Name of region the perforation will be connected to
+ string m_targetRegionName;
};
} // namespace geos
diff --git a/src/coreComponents/mesh/WellElementSubRegion.cpp b/src/coreComponents/mesh/WellElementSubRegion.cpp
index cd8d5cde24..f2eaf65cfb 100644
--- a/src/coreComponents/mesh/WellElementSubRegion.cpp
+++ b/src/coreComponents/mesh/WellElementSubRegion.cpp
@@ -120,7 +120,8 @@ void collectLocalAndBoundaryNodes( LineBlockABC const & lineBlock,
* @param[in] ei the index of the reservoir element
* @param[inout] nodes the nodes that have already been visited
*/
-void collectElementNodes( CellElementSubRegion const & subRegion,
+template< typename SUBREGION_TYPE >
+void collectElementNodes( SUBREGION_TYPE const & subRegion,
localIndex ei,
SortedArray< localIndex > & nodes )
{
@@ -137,6 +138,33 @@ void collectElementNodes( CellElementSubRegion const & subRegion,
}
}
+template< typename SUBREGION_TYPE >
+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
+ return false;
+}
+
+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();
+ return computationalGeometry::isPointInsidePolyhedron( referencePosition,
+ elemsToFaces[eiLocal],
+ facesToNodes,
+ elemCenter,
+ location );
+}
+
/**
* @brief Search the reservoir elements that can be accessed from the set "nodes".
Stop if a reservoir element containing the perforation is found.
@@ -221,7 +249,7 @@ bool visitNeighborElements( MeshLevel const & mesh,
erMatched = er;
esrMatched = esr;
eiMatched = eiLocal;
- giMatched = eiGlobal;
+ giMatched = eiGlobal;
matched = true;
break;
}
@@ -243,6 +271,111 @@ bool visitNeighborElements( MeshLevel const & mesh,
return matched;
}
+/**
+ * @brief Search the reservoir elements that can be accessed from the set "nodes".
+ Stop if a reservoir element containing the perforation is found.
+ If not, enlarge the set "nodes"
+ * @param[in] meshLevel the mesh object (single level only)
+ * @param[in] location the location of that we are trying to match with a reservoir element
+ * @param[inout] nodes the nodes that have already been visited
+ * @param[inout] elements the reservoir elements that have already been visited
+ * @param[in] targetRegionIndex the target region index for the reservoir element
+ * @param[in] targetSubRegionIndex the target subregion index for the reservoir element
+ * @param[inout] eiMatched the element index of the reservoir element that contains "location", if any
+ * @param[inout] giMatched the element global index of the reservoir element that contains "location", if any
+ */
+template< typename SUBREGION_TYPE >
+bool visitNeighborElements( MeshLevel const & mesh,
+ real64 const (&location)[3],
+ SortedArray< localIndex > & nodes,
+ SortedArray< globalIndex > & elements,
+ localIndex const & targetRegionIndex,
+ localIndex const & targetSubRegionIndex,
+ localIndex & eiMatched,
+ globalIndex & giMatched )
+{
+ ElementRegionManager const & elemManager = mesh.getElemManager();
+ NodeManager const & nodeManager = mesh.getNodeManager();
+ FaceManager const & faceManager = mesh.getFaceManager();
+
+ ArrayOfArraysView< localIndex const > const & toElementRegionList = nodeManager.elementRegionList();
+ ArrayOfArraysView< localIndex const > const & toElementSubRegionList = nodeManager.elementSubRegionList();
+ ArrayOfArraysView< localIndex const > const & toElementList = nodeManager.elementList();
+
+ arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const referencePosition =
+ nodeManager.referencePosition().toViewConst();
+
+ ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();
+
+ bool matched = false;
+
+ // In this function, we loop over all the nodes that we have collected so far.
+ // For each node, we grab all the reservoir elements that contain the node
+ // For reservoir element that has not been visited yet, we check if it contains "location"
+ // If the reservoir element contains "location" we can stop the search
+ // If the reservoir element does not contain "location", we collect its nodes
+
+ // we will enlarge the set of nodes in the loop below
+ // to do this we have to create a new set, "currNodes"
+ // that contains only the nodes that have already been visited
+ // the newly added nodes will be added to "nodes"
+ SortedArray< localIndex > currNodes = nodes;
+ giMatched = -1;
+ // for all the nodes already visited
+ for( localIndex currNode : currNodes )
+ {
+ // collect the elements that have not been visited yet
+ for( localIndex b = 0; b < toElementRegionList.sizeOfArray( currNode ); ++b )
+ {
+ localIndex const er = toElementRegionList[currNode][b];
+ localIndex const esr = toElementSubRegionList[currNode][b];
+ localIndex const eiLocal = toElementList[currNode][b];
+
+ if( er != targetRegionIndex || esr != targetSubRegionIndex )
+ continue;
+
+ 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 ) )
+ {
+ eiMatched = eiLocal;
+ giMatched = eiGlobal;
+ matched = true;
+ break;
+ }
+ // otherwise add the nodes of this element to the set of new nodes to visit
+ else
+ {
+ collectElementNodes( subRegion, eiLocal, nodes );
+ }
+ }
+ }
+
+ if( matched )
+ {
+ break;
+ }
+ }
+
+ // if not matched, insert the new nodes
+ return matched;
+}
+
/**
* @brief Search for the reservoir element that is the *closest* from the center of well element.
Note that this reservoir element does not necessarily contain the center of the well element.
@@ -279,7 +412,42 @@ void initializeLocalSearch( MeshLevel const & mesh,
erInit = std::get< 0 >( ret.second );
esrInit = std::get< 1 >( ret.second );
eiInit = std::get< 2 >( ret.second );
+}
+
+/**
+ * @brief Search for the reservoir element that is the *closest* from the center of well element.
+ Note that this reservoir element does not necessarily contain the center of the well element.
+ This "init" reservoir element will be used in SearchLocalElements to find the reservoir element that
+ contains the well element.
+ * @param[in] meshLevel the mesh object (single level only)
+ * @param[in] location the location of that we are trying to match with a reservoir element
+ * @param[in] targetRegionIndex the region index of the reservoir element from which we start the search
+ * @param[in] targetSubRegionIndex the subregion index of the reservoir element from which we start the search
+ * @param[inout] eiInit the element index of the reservoir element from which we start the search
+ */
+void initializeLocalSearch( MeshLevel const & mesh,
+ real64 const (&location)[3],
+ localIndex const & targetRegionIndex,
+ localIndex const & targetSubRegionIndex,
+ 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() );
+ // to initialize the local search for the reservoir element that contains "location",
+ // we find the reservoir element that minimizes the distance from "location" to the reservoir element center
+ auto ret = minLocOverElemsInSubRegion( subRegion, [&] ( localIndex const ei )
+ {
+ real64 v[3] = { location[0], location[1], location[2] };
+ LvArray::tensorOps::subtract< 3 >( v, resElemCenter[targetRegionIndex][targetSubRegionIndex][ei] );
+ auto dist = LvArray::tensorOps::l2Norm< 3 >( v );
+ return dist;
+ } );
+ // save the index of the reservoir element
+ // note that this reservoir element does not necessarily contains "location"
+ eiInit = ret.second;
}
/**
@@ -345,6 +513,91 @@ bool searchLocalElements( MeshLevel const & mesh,
return resElemFound;
}
+/**
+ * @brief Search for the reservoir element that contains the well element.
+ To do that, loop over the reservoir elements that are in the neighborhood of (erInit,esrInit,eiInit)
+ * @param[in] meshLevel the mesh object (single level only)
+ * @param[in] location the location of that we are trying to match with a reservoir element
+ * @param[in] targetRegionIndex the target region index for the reservoir element
+ * @param[inout] esrMatched the subregion index of the reservoir element that contains "location", if any
+ * @param[inout] eiMatched the element index of the reservoir element that contains "location", if any
+ * @param[inout] giMatched the element global index of the reservoir element that contains "location", if any
+ */
+bool searchLocalElements( MeshLevel const & mesh,
+ real64 const (&location)[3],
+ localIndex const & searchDepth,
+ localIndex const & targetRegionIndex,
+ localIndex & esrMatched,
+ localIndex & eiMatched,
+ globalIndex & giMatched )
+{
+ ElementRegionBase const & region = mesh.getElemManager().getRegion< ElementRegionBase >( targetRegionIndex );
+
+ bool resElemFound = false;
+ for( localIndex esr = 0; esr < region.numSubRegions(); ++esr )
+ {
+ ElementSubRegionBase const & subRegionBase = region.getSubRegion( esr );
+ region.applyLambdaToContainer< CellElementSubRegion, SurfaceElementSubRegion >( subRegionBase, [&]( auto const & subRegion )
+ {
+ GEOS_LOG( GEOS_FMT( " searching well connections with region/subregion: {}/{}", region.getName(), subRegion.getName() ) );
+
+ // first, we search for the reservoir element that is the *closest* from the center of well element
+ // note that this reservoir element does not necessarily contain the center of the well element
+ // this "init" reservoir element will be used later to find the reservoir element that
+ // contains the well element
+ localIndex eiInit = -1;
+ initializeLocalSearch( mesh, location, targetRegionIndex, esr, eiInit );
+
+ // loop over the reservoir elements that are in the neighborhood of (esrInit,eiInit)
+ // search locally, starting from the location of the previous perforation
+ // the assumption here is that perforations have been entered in order of depth
+
+ SortedArray< localIndex > nodes;
+ SortedArray< globalIndex > elements;
+
+ // here is how the search is done:
+ // 1 - We check if "location" is within the "init" reservoir element defined by (erInit,esrMatched,eiMatched)
+ // 2 - If yes, stop
+ // - If not, a) collect the nodes of the reservoir element defined by (erInit,esrMatched,eiMatched)
+ // b) use these nodes to grab the neighbors of (erInit,esrMatched,eiMatched)
+ // c) check if "location" is within the neighbors. If not, grab the neighbors of the neighbors, and so
+ // on...
+
+ // collect the nodes of the current element
+ // they will be used to access the neighbors and check if they contain the perforation
+ collectElementNodes( subRegion, eiInit, nodes );
+
+ // if no match is found, enlarge the neighborhood m_searchDepth'th times
+ for( localIndex d = 0; d < searchDepth; ++d )
+ {
+ localIndex nNodes = nodes.size();
+
+ // search the reservoir elements that can be accessed from the set "nodes"
+ // stop if a reservoir element containing the perforation is found
+ // if not, enlarge the set "nodes"
+
+ resElemFound =
+ visitNeighborElements< TYPEOFREF( subRegion ) >( mesh, location, nodes, elements, targetRegionIndex, esr, eiMatched, giMatched );
+
+ if( resElemFound || nNodes == nodes.size())
+ {
+ esrMatched = esr;
+ GEOS_LOG( GEOS_FMT( " found {}/{}/{}", region.getName(), subRegion.getName(), eiMatched ) );
+ // TODO learn how to exit forElementSubRegionsIndex
+ break;
+ }
+ }
+ } );
+
+ if( resElemFound )
+ {
+ break;
+ }
+ }
+
+ return resElemFound;
+}
+
}
void WellElementSubRegion::generate( MeshLevel & mesh,
@@ -448,7 +701,7 @@ void WellElementSubRegion::generate( MeshLevel & mesh,
void WellElementSubRegion::assignUnownedElementsInReservoir( MeshLevel & mesh,
LineBlockABC const & lineBlock,
- SortedArray< globalIndex > const & unownedElems,
+ SortedArray< globalIndex > const & unownedElems,
SortedArray< globalIndex > & localElems,
arrayView1d< integer > & elemStatusGlobal ) const
{
@@ -770,59 +1023,66 @@ void WellElementSubRegion::connectPerforationsToMeshElements( MeshLevel & mesh,
arrayView2d< real64 const > const perfCoordsGlobal = lineBlock.getPerfCoords();
arrayView1d< real64 const > const perfWellTransmissibilityGlobal = lineBlock.getPerfTransmissibility();
arrayView1d< real64 const > const perfWellSkinFactorGlobal = lineBlock.getPerfSkinFactor();
+ arrayView1d< string const > const perfTargetTegionGlobal = lineBlock.getPerfTargetRegion();
m_perforationData.resize( perfCoordsGlobal.size( 0 ) );
localIndex iperfLocal = 0;
arrayView2d< real64 > const perfLocation = m_perforationData.getLocation();
+ ElementRegionManager const & elemManager = mesh.getElemManager();
+
// loop over all the perforations
for( globalIndex iperfGlobal = 0; iperfGlobal < perfCoordsGlobal.size( 0 ); ++iperfGlobal )
{
real64 const location[3] = { perfCoordsGlobal[iperfGlobal][0],
perfCoordsGlobal[iperfGlobal][1],
perfCoordsGlobal[iperfGlobal][2] };
+ GEOS_LOG( GEOS_FMT( "{}: perforation {} location = ({}, {}, {})", lineBlock.getName(), iperfGlobal,
+ location[0], location[1], location[2] ) );
- localIndex erMatched = -1;
- localIndex esrMatched = -1;
- localIndex eiMatched = -1;
+ localIndex erStart = -1, erEnd = -1;
- localIndex erInit = -1;
- localIndex esrInit = -1;
- localIndex eiInit = -1;
- globalIndex giMatched = -1;
+ localIndex const targetRegionIndex = elemManager.getRegions().getIndex( perfTargetTegionGlobal[iperfGlobal] );
+ if( targetRegionIndex >= 0 )
+ {
+ erStart = targetRegionIndex;
+ erEnd = erStart + 1;
+ }
+ else // default is all regions
+ {
+ erStart = 0;
+ erEnd = elemManager.numRegions();
+ }
// for each perforation, we have to find the reservoir element that contains the perforation
+ for( localIndex er = erStart; er < erEnd; er++ )
+ {
+ // search for the reservoir element that contains the well element
+ localIndex esrMatched = -1;
+ localIndex eiMatched = -1;
+ globalIndex giMatched = -1;
+ bool resElemFound = searchLocalElements( mesh, location, m_searchDepth, er, esrMatched, eiMatched, giMatched );
+
+ // if the element was found
+ if( resElemFound )
+ {
+ // set the indices for the matched reservoir element
+ m_perforationData.getMeshElements().m_toElementRegion[iperfLocal] = er;
+ m_perforationData.getMeshElements().m_toElementSubRegion[iperfLocal] = esrMatched;
+ m_perforationData.getMeshElements().m_toElementIndex[iperfLocal] = eiMatched;
+ m_perforationData.getReservoirElementGlobalIndex()[iperfLocal] = giMatched;
- // Step 1: first, we search for the reservoir element that is the *closest* from the center of well element
- // note that this reservoir element does not necessarily contain the center of the well element
- // this "init" reservoir element will be used in SearchLocalElements to find the reservoir element that
- // contains the well element
- initializeLocalSearch( mesh, location,
- erInit, esrInit, eiInit );
+ // construct the local wellTransmissibility and location maps
+ m_perforationData.getWellTransmissibility()[iperfLocal] = perfWellTransmissibilityGlobal[iperfGlobal];
+ m_perforationData.getWellSkinFactor()[iperfLocal] = perfWellSkinFactorGlobal[iperfGlobal];
+ LvArray::tensorOps::copy< 3 >( perfLocation[iperfLocal], location );
- // Step 2: then, search for the reservoir element that contains the well element
- // to do that, we loop over the reservoir elements that are in the neighborhood of (erInit,esrInit,eiInit)
- bool resElemFound = searchLocalElements( mesh, location, m_searchDepth,
- erInit, esrInit, eiInit,
- erMatched, esrMatched, eiMatched, giMatched );
+ // increment the local to global map
+ m_perforationData.localToGlobalMap()[iperfLocal++] = iperfGlobal;
- // if the element was found
- if( resElemFound )
- {
- // set the indices for the matched reservoir element
- m_perforationData.getMeshElements().m_toElementRegion[iperfLocal] = erMatched;
- m_perforationData.getMeshElements().m_toElementSubRegion[iperfLocal] = esrMatched;
- m_perforationData.getMeshElements().m_toElementIndex[iperfLocal] = eiMatched;
- m_perforationData.getReservoirElementGlobalIndex()[iperfLocal] = giMatched;
-
- // construct the local wellTransmissibility and location maps
- m_perforationData.getWellTransmissibility()[iperfLocal] = perfWellTransmissibilityGlobal[iperfGlobal];
- m_perforationData.getWellSkinFactor()[iperfLocal] = perfWellSkinFactorGlobal[iperfGlobal];
- LvArray::tensorOps::copy< 3 >( perfLocation[iperfLocal], location );
-
- // increment the local to global map
- m_perforationData.localToGlobalMap()[iperfLocal++] = iperfGlobal;
+ break;
+ }
}
}
diff --git a/src/coreComponents/mesh/generators/LineBlock.hpp b/src/coreComponents/mesh/generators/LineBlock.hpp
index e2af5724ec..6d66051081 100644
--- a/src/coreComponents/mesh/generators/LineBlock.hpp
+++ b/src/coreComponents/mesh/generators/LineBlock.hpp
@@ -154,6 +154,7 @@ class LineBlock : public LineBlockABC
arrayView1d< real64 const > getPerfSkinFactor() const override final { return m_perfSkinFactor; }
+ arrayView1d< string const > getPerfTargetRegion() const override final { return m_perfTargetRegion; }
/**
* @brief Set the well skin factor at the perforations.
@@ -161,6 +162,12 @@ class LineBlock : public LineBlockABC
*/
void setPerfSkinFactor( arrayView1d< real64 const > perfSkinFactor ) { m_perfSkinFactor = perfSkinFactor; }
+ /**
+ * @brief Set the target region for the perforations.
+ * @param perfTargetRegion list of target regions for all the perforations on the well
+ */
+ void setPerfTargetRegion( arrayView1d< string const > perfTargetRegion ) { m_perfTargetRegion = perfTargetRegion; }
+
arrayView1d< globalIndex const > getPerfElemIndex() const override final { return m_perfElemId; }
/**
@@ -242,6 +249,9 @@ class LineBlock : public LineBlockABC
/// Well skin factor at the perforation
array1d< real64 > m_perfSkinFactor;
+ /// Target region for the perforation
+ array1d< string > m_perfTargetRegion;
+
/// Global index of the well element
array1d< globalIndex > m_perfElemId;
diff --git a/src/coreComponents/mesh/generators/LineBlockABC.hpp b/src/coreComponents/mesh/generators/LineBlockABC.hpp
index 2855789ed6..ad56b5b1ca 100644
--- a/src/coreComponents/mesh/generators/LineBlockABC.hpp
+++ b/src/coreComponents/mesh/generators/LineBlockABC.hpp
@@ -145,6 +145,12 @@ class LineBlockABC : public dataRepository::Group
*/
virtual arrayView1d< real64 const > getPerfSkinFactor() const = 0;
+ /**
+ * @brief Get the target region for the perforations.
+ * @return list of target regions for all the perforations on the well
+ */
+ virtual arrayView1d< string const > getPerfTargetRegion() const = 0;
+
/**
* @brief Get the global indices of the well elements connected to each perforation.
* @return list providing the global index of the connected well element for each perforation
diff --git a/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp b/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp
index 555bef9216..8e850479e6 100644
--- a/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp
+++ b/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp
@@ -89,6 +89,7 @@ void MeshGeneratorBase::attachWellInfo( CellBlockManager & cellBlockManager )
lb.setPerfCoords( wellGen.getPerfCoords() );
lb.setPerfTransmissibility( wellGen.getPerfTransmissibility() );
lb.setPerfSkinFactor( wellGen.getPerfSkinFactor() );
+ lb.setPerfTargetRegion( wellGen.getPerfTargetRegion() );
lb.setPerfElemIndex( wellGen.getPerfElemIndex() );
lb.setWellControlsName( wellGen.getWellControlsName() );
lb.setWellGeneratorName( wellGen.getName() );
diff --git a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp
index f64e0e9b7d..a49439667a 100644
--- a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp
+++ b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp
@@ -118,6 +118,7 @@ void WellGeneratorBase::generateWellGeometry( )
m_perfDistFromHead.resize( m_numPerforations );
m_perfTransmissibility.resize( m_numPerforations );
m_perfSkinFactor.resize( m_numPerforations );
+ m_perfTargetRegion.resize( m_numPerforations );
m_perfElemId.resize( m_numPerforations );
// construct a reverse map from the polyline nodes to the segments
@@ -348,6 +349,7 @@ void WellGeneratorBase::connectPerforationsToWellElements()
m_perfDistFromHead[iperf] = perf.getDistanceFromWellHead();
m_perfTransmissibility[iperf] = perf.getWellTransmissibility();
m_perfSkinFactor[iperf] = perf.getWellSkinFactor();
+ m_perfTargetRegion[iperf] = perf.getTargetRegion();
// search in all the elements of this well between head and bottom
globalIndex iwelemTop = 0;
diff --git a/src/coreComponents/mesh/generators/WellGeneratorBase.hpp b/src/coreComponents/mesh/generators/WellGeneratorBase.hpp
index fbe05157d6..db7505d8f0 100644
--- a/src/coreComponents/mesh/generators/WellGeneratorBase.hpp
+++ b/src/coreComponents/mesh/generators/WellGeneratorBase.hpp
@@ -185,6 +185,12 @@ class WellGeneratorBase : public MeshComponentBase
*/
arrayView1d< real64 const > getPerfSkinFactor() const { return m_perfSkinFactor; };
+ /**
+ * @brief Get the target region for a perforation.
+ * @return the target regions for a perforation
+ */
+ arrayView1d< string const > getPerfTargetRegion() const { return m_perfTargetRegion; };
+
/**
* @brief Get the global indices of the well elements connected to each perforation.
* @return list providing the global index of the connected well element for each perforation
@@ -393,6 +399,9 @@ class WellGeneratorBase : public MeshComponentBase
/// Physical location of the perforation wrt to well head
array1d< real64 > m_perfDistFromHead;
+ /// Target region for the perforation
+ array1d< string > m_perfTargetRegion;
+
};
}
#endif /* GEOS_MESH_GENERATORS_WELLGENERATORBASE_HPP_ */
diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp
index 9a76b95632..ada7fe4d2d 100644
--- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp
+++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp
@@ -30,7 +30,6 @@
#include "discretizationMethods/NumericalMethodsManager.hpp"
#include "mainInterface/ProblemManager.hpp"
#include "mesh/SurfaceElementRegion.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
#include "mesh/mpiCommunications/NeighborCommunicator.hpp"
#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" // needed to register pressure(_n)
#include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp"
diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp
index c9fe35f614..a9589a4424 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp
@@ -25,7 +25,6 @@
#include "finiteElement/Kinematics.h"
#include "finiteElement/FiniteElementDispatch.hpp"
#include "mesh/DomainPartition.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
#include "mesh/utilities/ComputationalGeometry.hpp"
namespace geos
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
index 70cb8b3c56..56fc98e9e6 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
@@ -33,6 +33,7 @@
#include "fieldSpecification/FieldSpecificationManager.hpp"
#include "fieldSpecification/TractionBoundaryCondition.hpp"
#include "finiteElement/FiniteElementDiscretizationManager.hpp"
+#include "finiteElement/FiniteElementDiscretization.hpp"
#include "LvArray/src/output.hpp"
#include "mainInterface/ProblemManager.hpp"
#include "mesh/DomainPartition.hpp"
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp
index 76d98cb514..ee5d8049f0 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp
@@ -24,7 +24,6 @@
#include "common/TimingMacros.hpp"
#include "kernels/SolidMechanicsLagrangianFEMKernels.hpp"
#include "kernels/StrainHelper.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
#include "mesh/mpiCommunications/CommunicationTools.hpp"
#include "mesh/mpiCommunications/MPI_iCommData.hpp"
#include "physicsSolvers/PhysicsSolverBase.hpp"
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
index a2b1d3ec2f..5eeb811ca8 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
@@ -31,6 +31,7 @@
#include "constitutive/ConstitutiveManager.hpp"
#include "constitutive/contact/FrictionBase.hpp"
#include "finiteElement/FiniteElementDiscretizationManager.hpp"
+#include "finiteElement/FiniteElementDiscretization.hpp"
#include "finiteElement/Kinematics.h"
#include "LvArray/src/output.hpp"
#include "mesh/DomainPartition.hpp"
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp
index 571e365f14..672ce06c89 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp
@@ -24,7 +24,6 @@
#include "common/TimingMacros.hpp"
#include "kernels/SolidMechanicsLagrangianFEMKernels.hpp"
#include "kernels/ExplicitMPM.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
#include "mesh/mpiCommunications/CommunicationTools.hpp"
#include "mesh/mpiCommunications/MPI_iCommData.hpp"
#include "physicsSolvers/PhysicsSolverBase.hpp"
diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp
index da4e285636..dd8fe2d42c 100644
--- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp
+++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp
@@ -191,7 +191,7 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups()
if( added )
{
- GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator, "Element " << cellIndex << " is fractured" );
+ GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator, GEOS_FMT( "{}: element {} is fractured", subRegion.getName(), cellIndex ) );
// Add the information to the CellElementSubRegion
subRegion.addFracturedElement( cellIndex, localNumberOfSurfaceElems );
diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
index 2e16f5bfee..f3711ab557 100644
--- a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
+++ b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
@@ -22,7 +22,7 @@
#include "mesh/mpiCommunications/CommunicationTools.hpp"
#include "mesh/mpiCommunications/NeighborCommunicator.hpp"
#include "mesh/mpiCommunications/SpatialPartition.hpp"
-#include "finiteElement/FiniteElementDiscretizationManager.hpp"
+#include "finiteElement/FiniteElementDiscretization.hpp"
#include "finiteVolume/FiniteVolumeManager.hpp"
#include "finiteVolume/FluxApproximationBase.hpp"
#include "mesh/SurfaceElementRegion.hpp"