diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 9511aa9fbe3..4e065cceae2 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3843-13975-76d125d + baseline: integratedTests/baseline_integratedTests-pr3359-14024-3c9ebb4 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index eae23a6e3dc..2e3a8b70338 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3359 (2025-10-07) +===================== +Add functions to connect well perforation to surface elements. + PR #3843 (2025-10-02) ===================== Fix linear solver issues. diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_base.xml similarity index 95% rename from inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml rename to inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_base.xml index f91d3f1e8f1..c6bbab35273 100644 --- a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml +++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_base.xml @@ -29,7 +29,7 @@ + permeabilityComponents="{ 1.0e-14, 1.0e-14, 1.0e-14 }" /> + targetRegion="Fault" + distanceFromHead="1700"/> diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_smoke.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_fim_smoke.xml similarity index 93% rename from inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_smoke.xml rename to inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_fim_smoke.xml index 7d0ccf7fc71..4f25a4348d6 100644 --- a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_smoke.xml +++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_fim_smoke.xml @@ -1,6 +1,6 @@ - + + solverType="direct" + directParallel="1"/> + targetTotalRate="1e-3" /> + targetTotalRate="1e-3" /> diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_seq_smoke.xml similarity index 93% rename from inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml rename to inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_seq_smoke.xml index 39c5e7b1609..2b7244f3455 100644 --- a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml +++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_WellInFault_seq_smoke.xml @@ -1,6 +1,6 @@ - + + solverType="direct" + directParallel="1"/> @@ -60,7 +62,7 @@ control="BHP" referenceElevation="-2500" targetBHP="1e7" - targetTotalRate="1e10" /> + targetTotalRate="1e-3" /> + targetTotalRate="1e-3" /> diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp index f74a79d69a2..decbae981ef 100644 --- a/src/coreComponents/mainInterface/ProblemManager.cpp +++ b/src/coreComponents/mainInterface/ProblemManager.cpp @@ -863,11 +863,12 @@ void ProblemManager::generateMeshLevel( MeshLevel & meshLevel, elemRegionManager.forElementSubRegions< ElementSubRegionBase >( [&]( ElementSubRegionBase & subRegion ) { subRegion.setupRelatedObjectsInRelations( meshLevel ); + // TODO calling calculateElementGeometricQuantities for `FaceElementSubRegion` here is not very accurate: // `FaceElementSubRegion` has no node and therefore needs the nodes positions from the neighbor elements // in order to compute the geometric quantities. // And this point of the process, the ghosting has not been done and some elements of the `FaceElementSubRegion` - // can have no neighbor. Making impossible the computation, which is therfore postponed to after the ghosting. - if( isBaseMeshLevel && !dynamicCast< FaceElementSubRegion * >( &subRegion ) ) + // can have no neighbor. We still call it only to compute element centers to be used for well perforation computations. + if( isBaseMeshLevel ) // && !dynamicCast< FaceElementSubRegion * >( &subRegion ) ) { subRegion.calculateElementGeometricQuantities( nodeManager, faceManager ); } diff --git a/src/coreComponents/mesh/DomainPartition.cpp b/src/coreComponents/mesh/DomainPartition.cpp index 07a79316822..943c6366974 100644 --- a/src/coreComponents/mesh/DomainPartition.cpp +++ b/src/coreComponents/mesh/DomainPartition.cpp @@ -247,9 +247,10 @@ 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.sortAllFaceNodes( nodeManager, elemManager ); faceManager.computeGeometry( nodeManager ); } else if( !meshLevel.isShallowCopyOf( meshBody.getMeshLevels().getGroup< MeshLevel >( 0 )) ) diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index f8692a95664..0ba6bfc2b84 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -203,7 +203,11 @@ void ElementRegionManager::generateWells( CellBlockManagerABC const & cellBlockM // generate the local data (well elements, nodes, perforations) on this well // note: each MPI rank knows the global info on the entire well (constructed earlier in InternalWellGenerator) // so we only need node and element offsets to construct the local-to-global maps in each wellElemSubRegion - wellRegion.generateWell( meshLevel, lineBlock, nodeOffsetGlobal + wellNodeCount, elemOffsetGlobal + wellElemCount ); + wellRegion.generateWell( meshLevel, + lineBlock, + nodeOffsetGlobal + wellNodeCount, + elemOffsetGlobal + wellElemCount, + cellBlockManager.getGlobalLength() ); // increment counters with global number of nodes and elements wellElemCount += lineBlock.numElements(); diff --git a/src/coreComponents/mesh/ElementSubRegionBase.hpp b/src/coreComponents/mesh/ElementSubRegionBase.hpp index c61289839c8..21810289001 100644 --- a/src/coreComponents/mesh/ElementSubRegionBase.hpp +++ b/src/coreComponents/mesh/ElementSubRegionBase.hpp @@ -279,14 +279,45 @@ class ElementSubRegionBase : public ObjectManagerBase forAll< parallelHostPolicy >( size(), [=]( localIndex const k ) { - LvArray::tensorOps::copy< 3 >( elementCenters[k], X[e2n( k, 0 )] ); + // collect node coordinates for element k localIndex const numNodes = this->numNodesPerElement( k ); - for( localIndex a = 1; a < numNodes; ++a ) + std::vector< std::array< real64, 3 > > nodes; + nodes.reserve( numNodes ); + for( localIndex a = 0; a < numNodes; ++a ) { - LvArray::tensorOps::add< 3 >( elementCenters[k], X[e2n( k, a )] ); + localIndex const nA = e2n( k, a ); + nodes.push_back( { X( nA, 0 ), X( nA, 1 ), X( nA, 2 ) } ); } - LvArray::tensorOps::scale< 3 >( elementCenters[k], 1.0 / numNodes ); + // sort + unique with tolerance to remove near-duplicates (stable and robust for floating point) + auto cmpLex = []( auto const & A, auto const & B ) + { + return std::tie( A[0], A[1], A[2] ) < std::tie( B[0], B[1], B[2] ); + }; + const real64 tol = 1e-10; + auto almostEqual = [tol]( auto const & A, auto const & B ) + { + return ( std::abs( A[0] - B[0] ) < tol ) && + ( std::abs( A[1] - B[1] ) < tol ) && + ( std::abs( A[2] - B[2] ) < tol ); + }; + + std::sort( nodes.begin(), nodes.end(), cmpLex ); + nodes.erase( std::unique( nodes.begin(), nodes.end(), almostEqual ), nodes.end() ); + + // compute average (center) of unique nodes + std::array< real64, 3 > centre = { 0.0, 0.0, 0.0 }; + for( auto const & p : nodes ) + { + centre[0] += p[0]; + centre[1] += p[1]; + centre[2] += p[2]; + } + const real64 inv = 1.0 / static_cast< real64 >( nodes.size() ); + + elementCenters[k][0] = centre[0] * inv; + elementCenters[k][1] = centre[1] * inv; + elementCenters[k][2] = centre[2] * inv; } ); } }; diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index 875d23b289f..67ed74c2f87 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -48,14 +48,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." ); diff --git a/src/coreComponents/mesh/FaceElementSubRegion.cpp b/src/coreComponents/mesh/FaceElementSubRegion.cpp index b7ee007a270..041a9789751 100644 --- a/src/coreComponents/mesh/FaceElementSubRegion.cpp +++ b/src/coreComponents/mesh/FaceElementSubRegion.cpp @@ -230,7 +230,7 @@ void FaceElementSubRegion::calculateSingleElementGeometricQuantities( localIndex m_elementVolume[k] = m_elementAperture[k] * faceArea[m_toFacesRelation[k][0]]; } -void FaceElementSubRegion::calculateElementGeometricQuantities( NodeManager const & GEOS_UNUSED_PARAM( nodeManager ), +void FaceElementSubRegion::calculateElementGeometricQuantities( NodeManager const & nodeManager, FaceManager const & faceManager ) { arrayView1d< real64 const > const & faceArea = faceManager.faceArea(); @@ -239,6 +239,9 @@ void FaceElementSubRegion::calculateElementGeometricQuantities( NodeManager cons { calculateSingleElementGeometricQuantities( k, faceArea ); } ); + + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X = nodeManager.referencePosition(); + calculateElementCenters( X ); } ElementType FaceElementSubRegion::getElementType( localIndex ei ) const diff --git a/src/coreComponents/mesh/PerforationData.cpp b/src/coreComponents/mesh/PerforationData.cpp index 69d738abb95..fd0e4140e59 100644 --- a/src/coreComponents/mesh/PerforationData.cpp +++ b/src/coreComponents/mesh/PerforationData.cpp @@ -246,23 +246,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( isZero( dx ) ) + { + dx = subRegion.getElementVolume()[ei]; + dx /= dy * dz; + } + else if( isZero( dy ) ) + { + 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, diff --git a/src/coreComponents/mesh/WellElementRegion.cpp b/src/coreComponents/mesh/WellElementRegion.cpp index 695a0bb2377..e7b2f54e019 100644 --- a/src/coreComponents/mesh/WellElementRegion.cpp +++ b/src/coreComponents/mesh/WellElementRegion.cpp @@ -49,7 +49,8 @@ WellElementRegion::~WellElementRegion() void WellElementRegion::generateWell( MeshLevel & mesh, LineBlockABC const & lineBlock, globalIndex nodeOffsetGlobal, - globalIndex elemOffsetGlobal ) + globalIndex elemOffsetGlobal, + real64 const globalLength ) { // get the (unique) subregion WellElementSubRegion & @@ -64,8 +65,11 @@ void WellElementRegion::generateWell( MeshLevel & mesh, globalIndex const numElemsGlobal = lineBlock.numElements(); globalIndex const numPerforationsGlobal = lineBlock.numPerforations(); + // tolerance for geometrical calculations + real64 const geomTol = globalLength * 1e-10; + // 1) select the local perforations based on connectivity to the local reservoir elements - subRegion.connectPerforationsToMeshElements( mesh, lineBlock ); + subRegion.connectPerforationsToMeshElements( mesh, lineBlock, geomTol ); globalIndex const matchedPerforations = MpiWrapper::sum( perforationData->size() ); GEOS_THROW_IF( matchedPerforations != numPerforationsGlobal, @@ -105,7 +109,8 @@ void WellElementRegion::generateWell( MeshLevel & mesh, lineBlock, elemStatusGlobal, nodeOffsetGlobal, - elemOffsetGlobal ); + elemOffsetGlobal, + geomTol ); // 4) find out which rank is the owner of the top segment diff --git a/src/coreComponents/mesh/WellElementRegion.hpp b/src/coreComponents/mesh/WellElementRegion.hpp index 5a41cf19526..ee3e8be547b 100644 --- a/src/coreComponents/mesh/WellElementRegion.hpp +++ b/src/coreComponents/mesh/WellElementRegion.hpp @@ -130,11 +130,13 @@ class WellElementRegion : public ElementRegionBase * @param[in] lineBlock the LineBlockABC containing the global well topology * @param[in] nodeOffsetGlobal the offset of the first global well node ( = offset of last global mesh node + 1 ) * @param[in] elemOffsetGlobal the offset of the first global well element ( = offset of last global mesh elem + 1 ) + * @param[in] globalLength the global length of the mesh */ void generateWell( MeshLevel & mesh, LineBlockABC const & lineBlock, globalIndex nodeOffsetGlobal, - globalIndex elemOffsetGlobal ); + globalIndex elemOffsetGlobal, + real64 globalLength ); ///@} diff --git a/src/coreComponents/mesh/WellElementSubRegion.cpp b/src/coreComponents/mesh/WellElementSubRegion.cpp index 2442eb61240..1b61cc7c507 100644 --- a/src/coreComponents/mesh/WellElementSubRegion.cpp +++ b/src/coreComponents/mesh/WellElementSubRegion.cpp @@ -21,6 +21,7 @@ #include "common/MpiWrapper.hpp" #include "LvArray/src/output.hpp" +#include namespace geos { @@ -147,8 +148,8 @@ 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] ) + real64 const (&GEOS_UNUSED_PARAM( location ))[3], + real64 const GEOS_UNUSED_PARAM( geomTol ) ) { // only CellElementSubRegion is currently supported return false; @@ -158,15 +159,87 @@ 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] ) + real64 const (&location)[3], + real64 const geomTol ) { 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 ); + location, + geomTol ); +} + +// 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 +{ + PointsEqual( double tolerance ): m_tol( tolerance ) {} + + bool operator()( POINT_TYPE const & p1, POINT_TYPE const & p2 ) const + { + return (std::abs( p1[0] - p2[0] ) < m_tol) && (std::abs( p1[1] - p2[1] ) < m_tol) && (std::abs( p1[2] - p2[2] ) < m_tol); + } + +private: + real64 const m_tol; +}; + +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], + real64 const geomTol ) +{ + 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 + auto cmpLex = []( Point3d const & a, Point3d const & b ) + { + return std::tie( a[0], a[1], a[2] ) < std::tie( b[0], b[1], b[2] ); + }; + + auto almostEqual = [geomTol]( Point3d const & a, Point3d const & b ) + { + return (std::abs( a[0]-b[0] ) < geomTol) && (std::abs( a[1]-b[1] ) < geomTol) && (std::abs( a[2]-b[2] ) < geomTol); + }; + + std::sort( polygon.begin(), polygon.end(), cmpLex ); + auto it = std::unique( polygon.begin(), polygon.end(), almostEqual ); + polygon.erase( it, polygon.end() ); + + return computationalGeometry::isPointInPolygon3d( polygon, polygon.size(), location, geomTol ); } /** @@ -190,7 +263,8 @@ bool visitNeighborElements( MeshLevel const & mesh, localIndex const & targetRegionIndex, localIndex const & targetSubRegionIndex, localIndex & eiMatched, - globalIndex & giMatched ) + globalIndex & giMatched, + real64 const geomTol ) { ElementRegionManager const & elemManager = mesh.getElemManager(); NodeManager const & nodeManager = mesh.getNodeManager(); @@ -234,8 +308,6 @@ 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 @@ -243,13 +315,9 @@ bool visitNeighborElements( MeshLevel const & mesh, { 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, geomTol ) ) { eiMatched = eiLocal; giMatched = eiGlobal; @@ -326,7 +394,8 @@ bool searchLocalElements( MeshLevel const & mesh, localIndex const & targetRegionIndex, localIndex & esrMatched, localIndex & eiMatched, - globalIndex & giMatched ) + globalIndex & giMatched, + real64 const geomTol ) { ElementRegionBase const & region = mesh.getElemManager().getRegion< ElementRegionBase >( targetRegionIndex ); @@ -378,7 +447,15 @@ bool searchLocalElements( MeshLevel const & mesh, // if not, enlarge the set "nodes" resElemFound = - visitNeighborElements< TYPEOFREF( subRegion ) >( mesh, location, nodes, elements, targetRegionIndex, esr, eiMatched, giMatched ); + visitNeighborElements< TYPEOFREF( subRegion ) >( mesh, + location, + nodes, + elements, + targetRegionIndex, + esr, + eiMatched, + giMatched, + geomTol ); if( resElemFound || nNodes == nodes.size()) { @@ -407,7 +484,8 @@ void WellElementSubRegion::generate( MeshLevel & mesh, LineBlockABC const & lineBlock, arrayView1d< integer > & elemStatusGlobal, globalIndex nodeOffsetGlobal, - globalIndex elemOffsetGlobal ) + globalIndex elemOffsetGlobal, + real64 const geomTol ) { map< integer, SortedArray< globalIndex > > elemSetsByStatus; @@ -449,7 +527,8 @@ void WellElementSubRegion::generate( MeshLevel & mesh, lineBlock, unownedElems, localElems, - elemStatusGlobal ); + elemStatusGlobal, + geomTol ); // 1.b) Then we check that all the well elements have been assigned (and assigned once) // This is needed because if the center of the well element falls on the boundary of // a reservoir element, the assignment algorithm of 1.a) can assign the same well element @@ -508,7 +587,8 @@ void WellElementSubRegion::assignUnownedElementsInReservoir( MeshLevel & mesh, LineBlockABC const & lineBlock, SortedArray< globalIndex > const & unownedElems, SortedArray< globalIndex > & localElems, - arrayView1d< integer > & elemStatusGlobal ) const + arrayView1d< integer > & elemStatusGlobal, + real64 const geomTol ) const { ElementRegionManager const & elemManager = mesh.getElemManager(); // get the well and reservoir element coordinates @@ -530,7 +610,7 @@ void WellElementSubRegion::assignUnownedElementsInReservoir( MeshLevel & mesh, localIndex esrMatched = -1; localIndex eiMatched = -1; globalIndex giMatched = -1; - integer const resElemFound = searchLocalElements( mesh, location, m_searchDepth, er, esrMatched, eiMatched, giMatched ); + integer const resElemFound = searchLocalElements( mesh, location, m_searchDepth, er, esrMatched, eiMatched, giMatched, geomTol ); // if the element was found if( resElemFound ) @@ -811,7 +891,8 @@ void WellElementSubRegion::updateNodeManagerNodeToElementMap( MeshLevel & mesh ) } void WellElementSubRegion::connectPerforationsToMeshElements( MeshLevel & mesh, - LineBlockABC const & lineBlock ) + LineBlockABC const & lineBlock, + real64 const geomTol ) { arrayView2d< real64 const > const perfCoordsGlobal = lineBlock.getPerfCoords(); arrayView1d< real64 const > const perfWellTransmissibilityGlobal = lineBlock.getPerfTransmissibility(); @@ -857,7 +938,14 @@ void WellElementSubRegion::connectPerforationsToMeshElements( MeshLevel & mesh, localIndex esrMatched = -1; localIndex eiMatched = -1; globalIndex giMatched = -1; - integer const resElemFound = searchLocalElements( mesh, location, m_searchDepth, er, esrMatched, eiMatched, giMatched ); + integer const resElemFound = searchLocalElements( mesh, + location, + m_searchDepth, + er, + esrMatched, + eiMatched, + giMatched, + geomTol ); // if the element was found if( resElemFound ) diff --git a/src/coreComponents/mesh/WellElementSubRegion.hpp b/src/coreComponents/mesh/WellElementSubRegion.hpp index c70821d645f..26089a42838 100644 --- a/src/coreComponents/mesh/WellElementSubRegion.hpp +++ b/src/coreComponents/mesh/WellElementSubRegion.hpp @@ -245,20 +245,24 @@ class WellElementSubRegion : public ElementSubRegionBase * enum SegmentStatus. They are used to partition well elements. * @param[in] nodeOffsetGlobal the offset of the first global well node ( = offset of last global mesh node + 1 ) * @param[in] elemOffsetGlobal the offset of the first global well element ( = offset of last global mesh elem + 1 ) + * @param[in] geomTol the tolerance for geometric calculations */ void generate( MeshLevel & mesh, LineBlockABC const & lineBlock, arrayView1d< integer > & elemStatus, globalIndex nodeOffsetGlobal, - globalIndex elemOffsetGlobal ); + globalIndex elemOffsetGlobal, + real64 geomTol ); /** * @brief For each perforation, find the reservoir element that contains the perforation. * @param[in] mesh the mesh object (single level only) * @param[in] lineBlock the LineBlockABC containing the global well topology + * @param[in] geomTol the tolerance for geometric calculations */ void connectPerforationsToMeshElements( MeshLevel & mesh, - LineBlockABC const & lineBlock ); + LineBlockABC const & lineBlock, + real64 geomTol ); /** * @brief Reconstruct the (local) map nextWellElemId using nextWellElemIdGlobal after the ghost exchange. @@ -392,12 +396,14 @@ class WellElementSubRegion : public ElementSubRegionBase with the newly assigned well elements in this function. * @param[out] wellElemStatus list of current well element status. Status values are defined in * enum SegmentStatus. They are used to partition well elements. + * @param[in] geomTol the tolerance for geometric calculations */ void assignUnownedElementsInReservoir( MeshLevel & mesh, LineBlockABC const & lineBlock, SortedArray< globalIndex > const & unownedElems, SortedArray< globalIndex > & localElems, - arrayView1d< integer > & elemStatusGlobal ) const; + arrayView1d< integer > & elemStatusGlobal, + real64 geomTol ) const; /** * @brief Check that all the well elements have been assigned to a single rank. diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index fca008633b1..61de4181531 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -433,6 +433,145 @@ bool isPointInsidePolyhedron( arrayView2d< real64 const, nodes::REFERENCE_POSITI return true; } +/** + * @brief Check if a point is inside a polygon (2D version) + * @tparam POLYGON_TYPE type of @p polygon + * @tparam POINT_TYPE type of @p point + * @param[in] polygon array of polygon nodes coordinates + * @param[in] n number of polygon nodes + * @param[in] point coordinates of the query point + * @param[in] tol tolerance for coordinate comparisons + * @return whether the point is inside + */ +template< typename POLYGON_TYPE, typename POINT_TYPE > +bool isPointInPolygon2d( POLYGON_TYPE const & polygon, + integer n, + POINT_TYPE const & point, + real64 const tol = 1e-10 ) +{ + integer count = 0; + + for( integer i = 0; i < n; ++i ) + { + auto const & p1 = polygon[i]; + auto const & p2 = polygon[(i + 1) % n]; + + real64 y1 = p1[1], y2 = p2[1]; + real64 x1 = p1[0], x2 = p2[0]; + real64 py = point[1], px = point[0]; + + // quick reject in y with tolerance + if( py + tol < std::min( y1, y2 ) || py - tol > std::max( y1, y2 ) ) + continue; + + // check if point is (approximately) on the segment + // parametric t for projection on segment in y (if segment vertical-ish use x) + if( std::abs( (x2 - x1) * (py - y1) - (px - x1) * (y2 - y1) ) < tol * + ( std::hypot( x2 - x1, y2 - y1 ) + 1.0 ) ) + { + // ensure px is between x1,x2 and py between y1,y2 (with tol) + if( px + tol >= std::min( x1, x2 ) && px - tol <= std::max( x1, x2 ) && + py + tol >= std::min( y1, y2 ) && py - tol <= std::max( y1, y2 ) ) + return true; // on boundary -> consider inside + } + + // ignore nearly-horizontal edges for intersection counting + if( std::abs( y2 - y1 ) < tol ) + continue; + + // compute x coordinate of intersection of horizontal line py with segment p1-p2 + real64 xIntersect = x1 + (py - y1) * (x2 - x1) / (y2 - y1); + + // count crossing where intersection is strictly to the right of point (robust with tol) + if( px < xIntersect - tol ) + ++count; + } + + return (count % 2) == 1; +} + +/** + * @brief Check if a point is inside a polygon (3D version) + * @tparam POLYGON_TYPE type of @p polygon + * @tparam POINT_TYPE type of @p point + * @param[in] polygon array of polygon nodes coordinates + * @param[in] n number of polygon nodes + * @param[in] point coordinates of the query point + * @param[in] tol tolerance for coordinate comparisons + * @return whether the point is inside + */ +template< typename POLYGON_TYPE, typename POINT_TYPE > +bool isPointInPolygon3d( POLYGON_TYPE const & polygon, + integer const n, + POINT_TYPE const & point, + real64 const tol = 1e-10 ) +{ + // Check if the point lies in the plane of the polygon + auto const & p0 = polygon[0]; + POINT_TYPE normal = {0, 0, 0}; + for( integer i = 1; i < n - 1; i++ ) + { + auto const & p1 = polygon[i]; + auto const & p2 = polygon[i + 1]; + normal[0] += (p1[1] - p0[1]) * (p2[2] - p0[2]) - (p1[2] - p0[2]) * (p2[1] - p0[1]); + normal[1] += (p1[2] - p0[2]) * (p2[0] - p0[0]) - (p1[0] - p0[0]) * (p2[2] - p0[2]); + normal[2] += (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p1[1] - p0[1]) * (p2[0] - p0[0]); + } + + real64 const dist = normal[0] * point[0] + normal[1] * point[1] + normal[2] * point[2] -(normal[0] * p0[0] + normal[1] * p0[1] + normal[2] * p0[2]); + + if( std::abs( dist ) > tol ) + { + return false; + } + + // Determine the dominant component of the normal vector + int dominantIndex = 0; + if( std::abs( normal[1] ) > std::abs( normal[0] )) + { + dominantIndex = 1; + } + if( std::abs( normal[2] ) > std::abs( normal[dominantIndex] )) + { + dominantIndex = 2; + } + + // Project the polygon and the point onto a 2D plane + POLYGON_TYPE projectedPolygon( n ); + POINT_TYPE projectedPoint; + if( dominantIndex == 0 ) // X is dominant, project onto YZ plane + { + for( int i = 0; i < n; i++ ) + { + projectedPolygon[i][0] = polygon[i][1]; + projectedPolygon[i][1] = polygon[i][2]; + } + projectedPoint[0] = point[1]; + projectedPoint[1] = point[2]; + } + else if( dominantIndex == 1 ) // Y is dominant, project onto XZ plane + { + for( int i = 0; i < n; i++ ) + { + projectedPolygon[i][0] = polygon[i][0]; + projectedPolygon[i][1] = polygon[i][2]; + } + projectedPoint[0] = point[0]; + projectedPoint[1] = point[2]; + } + else // Z is dominant, project onto XY plane + { + for( int i = 0; i < n; i++ ) + { + projectedPolygon[i][0] = polygon[i][0]; + projectedPolygon[i][1] = polygon[i][1]; + } + projectedPoint[0] = point[0]; + projectedPoint[1] = point[1]; + } + + return isPointInPolygon2d( projectedPolygon, n, projectedPoint ); +} /** * @brief Method to perform lexicographic comparison of two nodes based on coordinates. @@ -641,7 +780,7 @@ int findTriangleRefElement( arraySlice1d< localIndex const > const & nodeElement } /** - * @brief Computes the winding number of a point with respecto to a mesh element. + * @brief Computes the winding number of a point with respect to a mesh element. * @tparam POINT_TYPE type of @p point * @param[in] element the element to be checked * @param[in] nodeCoordinates a global array of nodal coordinates @@ -839,10 +978,10 @@ bool isPointInsideConvexPolyhedronRobust( localIndex element, * @param[in] pointCoordinates the vertices coordinates. * @param[out] boxDims The dimensions of the bounding box. */ -template< typename VEC_TYPE > +template< typename NODE_MAP_TYPE, typename VEC_TYPE > GEOS_HOST_DEVICE void getBoundingBox( localIndex const elemIndex, - arrayView2d< localIndex const, cells::NODE_MAP_USD > const & pointIndices, + NODE_MAP_TYPE const & pointIndices, arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & pointCoordinates, VEC_TYPE && boxDims ) { @@ -855,7 +994,7 @@ void getBoundingBox( localIndex const elemIndex, LvArray::tensorOps::fill< 3 >( boxDims, LvArray::NumericLimits< real64 >::lowest ); // loop over all the vertices of the element to get the min and max coords - for( localIndex a = 0; a < pointIndices.size( 1 ); ++a ) + for( localIndex a = 0; a < pointIndices[elemIndex].size(); ++a ) { localIndex const id = pointIndices( elemIndex, a ); for( localIndex d = 0; d < 3; ++d )