From b7988a567d83b46c21fb061509ba1e4bfb801fa3 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Wed, 18 Sep 2024 09:17:10 -0500 Subject: [PATCH 01/10] wip --- ...Poromechanics_FaultModel_well_new_base.xml | 7 +- .../mesh/ElementSubRegionBase.hpp | 1 + .../mesh/EmbeddedSurfaceSubRegion.cpp | 8 -- .../mesh/SurfaceElementRegion.cpp | 4 + .../mesh/WellElementSubRegion.cpp | 128 ++++++++++++++++-- .../mesh/utilities/ComputationalGeometry.hpp | 1 - 6 files changed, 128 insertions(+), 21 deletions(-) diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml index f91d3f1e8f..0d7f4af036 100644 --- a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml +++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml @@ -188,14 +188,15 @@ logLevel="1" wellRegionName="wellRegion2" wellControlsName="wellControls2" - polylineNodeCoords="{ { -5400, -5400, 0 }, - { -5400, -5400, -2500 } }" + polylineNodeCoords="{ { 1200, -600, 0 }, + { 1200, -600, -2500 } }" polylineSegmentConn="{ { 0, 1 } }" radius="0.1" numElementsPerSegment="1"> + targetRegion="Fault" + distanceFromHead="1689"/> diff --git a/src/coreComponents/mesh/ElementSubRegionBase.hpp b/src/coreComponents/mesh/ElementSubRegionBase.hpp index 74bd710a31..76c2e115ed 100644 --- a/src/coreComponents/mesh/ElementSubRegionBase.hpp +++ b/src/coreComponents/mesh/ElementSubRegionBase.hpp @@ -287,6 +287,7 @@ class ElementSubRegionBase : public ObjectManagerBase } LvArray::tensorOps::scale< 3 >( elementCenters[k], 1.0 / numNodes ); + //std::cout << k << " " << elementCenters[k] << std::endl; } ); } }; diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index 83f603d999..17fb78a032 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/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." ); diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index 78233703bd..47f8993b32 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -122,8 +122,12 @@ localIndex SurfaceElementRegion::addToFractureMesh( real64 const time_np1, localIndex const kfe = subRegion.size() - 1; ruptureTime( kfe ) = time_np1; +//std::cout << kfe << " " << elemCenter[ kfe ] << " " << faceCenter[ faceIndices[ 0 ] ] << std::endl; + LvArray::tensorOps::copy< 3 >( elemCenter[ kfe ], faceCenter[ faceIndices[ 0 ] ] ); +//std::cout << kfe << " " << elemCenter[ kfe ] << " " << faceCenter[ faceIndices[ 0 ] ] << std::endl; + faceMap.resizeArray( kfe, 2 ); faceMap[kfe][0] = faceIndices[0]; faceMap[kfe][1] = faceIndices[1]; diff --git a/src/coreComponents/mesh/WellElementSubRegion.cpp b/src/coreComponents/mesh/WellElementSubRegion.cpp index dd4a4e6c8e..1555c66331 100644 --- a/src/coreComponents/mesh/WellElementSubRegion.cpp +++ b/src/coreComponents/mesh/WellElementSubRegion.cpp @@ -143,7 +143,6 @@ bool isPointInsideElement( SUBREGION_TYPE const & 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,10 +153,13 @@ 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, @@ -165,6 +167,108 @@ bool isPointInsideElement( CellElementSubRegion const & subRegion, location ); } +bool isPointInPolygon2D(real64** polygon, integer n, real64* point) { + integer count = 0; + + for (integer i = 0; i < n; i++) { + real64* p1 = polygon[i]; + real64* p2 = polygon[(i + 1) % n]; + + if ((point[1] > std::min(p1[1], p2[1])) && (point[1] <= std::max(p1[1], p2[1])) && (point[0] <= std::max(p1[0], p2[0]))) { + real64 xIntersect = (point[1] - p1[1]) * (p2[0] - p1[0]) / (p2[1] - p1[1]) + p1[0]; + if (std::abs(p1[0] - p2[0]) < 1e-10 || point[0] <= xIntersect) { + count++; + } + } + } + + return count % 2 == 1; +} + +bool isPointInPolygon3D(real64** polygon, integer n, real64* point) { + // Check if the point lies in the plane of the polygon + real64* p0 = polygon[0]; + real64 normal[3] = {0, 0, 0}; + for (integer i = 1; i < n - 1; i++) { + real64* p1 = polygon[i]; + real64* 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 d = -(normal[0] * p0[0] + normal[1] * p0[1] + normal[2] * p0[2]); + real64 dist = normal[0] * point[0] + normal[1] * point[1] + normal[2] * point[2] + d; + + if (std::abs(dist) > 1e-6) { + return false; + } + + // Project the polygon and the point onto a 2D plane + real64** projectedPolygon = new real64*[n]; + for (integer i = 0; i < n; i++) { + projectedPolygon[i] = new real64[2]; + projectedPolygon[i][0] = polygon[i][1]; + projectedPolygon[i][1] = polygon[i][2]; + } + real64 projectedPoint[2] = {point[1], point[2]}; + + bool result = isPointInPolygon2D(projectedPolygon, n, projectedPoint); + + for (integer i = 0; i < n; i++) { + delete[] projectedPolygon[i]; + } + delete[] projectedPolygon; + + return result; +} + +template< typename POINT_TYPE > +bool isPointInsidePolygon( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodeCoordinates, + localIndex const & eiLocal, + integer const & nV, + SurfaceElementSubRegion::NodeMapType const & facesToNodes, + POINT_TYPE const & point_ ) +{ + real64** polygon = new real64*[nV]; + real64* point = new real64[3]; + std::cout <<"checking point = "; + for( integer j = 0; j < 3; ++j ) + { + point[j] = point_[j]; + std::cout << point[j] << " "; + } + std::cout << " polygon = "; + for( integer i = 0; i < nV; ++i ) + { + polygon[i] = new real64[3]; + std::cout << "( "; + for( integer j = 0; j < 3; ++j ) + { + polygon[i][j] = nodeCoordinates[facesToNodes( eiLocal, i )][j]; + std::cout << polygon[i][j] << " "; + } + std::cout << ")"; + } + std::cout << std::endl; + return isPointInPolygon3D(polygon, nV, point); +} + +bool isPointInsideElement( SurfaceElementSubRegion const & subRegion, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & referencePosition, + localIndex const & eiLocal, + ArrayOfArraysView< localIndex const > const & facesToNodes, + real64 const (&location)[3] ) +{ + std::cout << eiLocal << " " << subRegion.numNodesPerElement( eiLocal ) << std::endl; + SurfaceElementSubRegion::NodeMapType const & nodeList = subRegion.nodeList(); + return isPointInsidePolygon( referencePosition, + eiLocal, + 4,//subRegion.numNodesPerElement( eiLocal ), + nodeList, + location ); +} + /** * @brief Search the reservoir elements that can be accessed from the set "nodes". Stop if a reservoir element containing the perforation is found. @@ -333,8 +437,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 @@ -342,13 +444,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 ) ) { eiMatched = eiLocal; matched = true; @@ -427,7 +525,16 @@ void initializeLocalSearch( MeshLevel const & mesh, localIndex const & targetSubRegionIndex, localIndex & eiInit ) { + { + mesh.getElemManager().getRegion( targetRegionIndex ).forElementSubRegionsIndex< CellElementSubRegion, FaceElementSubRegion >( [&] ( localIndex const esr, auto const & subRegion ) + { + if(targetSubRegionIndex == esr) + subRegion.calculateElementCenters(mesh.getNodeManager().referencePosition()); + } ); + } + 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() ); @@ -438,12 +545,15 @@ void initializeLocalSearch( MeshLevel const & mesh, 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 ); + std::cout << resElemCenter[targetRegionIndex][targetSubRegionIndex][ei] << " " << dist << std::endl; return dist; } ); // save the index of the reservoir element // note that this reservoir element does not necessarily contains "location" eiInit = ret.second; + + std::cout << "init " << targetRegionIndex << " " << targetSubRegionIndex << " " << eiInit << " " << resElemCenter[targetRegionIndex][targetSubRegionIndex][eiInit] << std::endl; } /** diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 23e3fe6588..fccdf22d97 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -433,7 +433,6 @@ bool isPointInsidePolyhedron( arrayView2d< real64 const, nodes::REFERENCE_POSITI return true; } - /** * @brief Method to perform lexicographic comparison of two nodes based on coordinates. * @tparam COORD_TYPE type of coordinate From 5e61ed8bcd27af7e8154205fbf3b88d638dcb841 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 19 Sep 2024 11:14:01 -0500 Subject: [PATCH 02/10] this should work --- ...Poromechanics_FaultModel_well_new_base.xml | 6 +- .../mesh/WellElementSubRegion.cpp | 136 ++++++------------ .../mesh/utilities/ComputationalGeometry.hpp | 101 +++++++++++++ 3 files changed, 146 insertions(+), 97 deletions(-) diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml index 0d7f4af036..f6c8fae153 100644 --- a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml +++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml @@ -188,15 +188,15 @@ logLevel="1" wellRegionName="wellRegion2" wellControlsName="wellControls2" - polylineNodeCoords="{ { 1200, -600, 0 }, - { 1200, -600, -2500 } }" + polylineNodeCoords="{ { 1200, -700, 0 }, + { 1200, -700, -2500 } }" polylineSegmentConn="{ { 0, 1 } }" radius="0.1" numElementsPerSegment="1"> + distanceFromHead="1700"/> diff --git a/src/coreComponents/mesh/WellElementSubRegion.cpp b/src/coreComponents/mesh/WellElementSubRegion.cpp index 1555c66331..c58399cd1e 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 { @@ -158,8 +159,8 @@ bool isPointInsideElement( CellElementSubRegion const & subRegion, 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] }; + elemCenters[eiLocal][1], + elemCenters[eiLocal][2] }; return computationalGeometry::isPointInsidePolyhedron( referencePosition, elemsToFaces[eiLocal], facesToNodes, @@ -167,92 +168,28 @@ bool isPointInsideElement( CellElementSubRegion const & subRegion, location ); } -bool isPointInPolygon2D(real64** polygon, integer n, real64* point) { - integer count = 0; - - for (integer i = 0; i < n; i++) { - real64* p1 = polygon[i]; - real64* p2 = polygon[(i + 1) % n]; - - if ((point[1] > std::min(p1[1], p2[1])) && (point[1] <= std::max(p1[1], p2[1])) && (point[0] <= std::max(p1[0], p2[0]))) { - real64 xIntersect = (point[1] - p1[1]) * (p2[0] - p1[0]) / (p2[1] - p1[1]) + p1[0]; - if (std::abs(p1[0] - p2[0]) < 1e-10 || point[0] <= xIntersect) { - count++; - } - } - } - - return count % 2 == 1; -} - -bool isPointInPolygon3D(real64** polygon, integer n, real64* point) { - // Check if the point lies in the plane of the polygon - real64* p0 = polygon[0]; - real64 normal[3] = {0, 0, 0}; - for (integer i = 1; i < n - 1; i++) { - real64* p1 = polygon[i]; - real64* 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 d = -(normal[0] * p0[0] + normal[1] * p0[1] + normal[2] * p0[2]); - real64 dist = normal[0] * point[0] + normal[1] * point[1] + normal[2] * point[2] + d; - - if (std::abs(dist) > 1e-6) { - return false; - } - - // Project the polygon and the point onto a 2D plane - real64** projectedPolygon = new real64*[n]; - for (integer i = 0; i < n; i++) { - projectedPolygon[i] = new real64[2]; - projectedPolygon[i][0] = polygon[i][1]; - projectedPolygon[i][1] = polygon[i][2]; - } - real64 projectedPoint[2] = {point[1], point[2]}; - - bool result = isPointInPolygon2D(projectedPolygon, n, projectedPoint); - - for (integer i = 0; i < n; i++) { - delete[] projectedPolygon[i]; - } - delete[] projectedPolygon; - - return result; -} +// 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 > -bool isPointInsidePolygon( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodeCoordinates, - localIndex const & eiLocal, - integer const & nV, - SurfaceElementSubRegion::NodeMapType const & facesToNodes, - POINT_TYPE const & point_ ) +struct PointsEqual { - real64** polygon = new real64*[nV]; - real64* point = new real64[3]; - std::cout <<"checking point = "; - for( integer j = 0; j < 3; ++j ) - { - point[j] = point_[j]; - std::cout << point[j] << " "; - } - std::cout << " polygon = "; - for( integer i = 0; i < nV; ++i ) + bool operator()( POINT_TYPE const & p1, POINT_TYPE const & p2 ) const { - polygon[i] = new real64[3]; - std::cout << "( "; - for( integer j = 0; j < 3; ++j ) - { - polygon[i][j] = nodeCoordinates[facesToNodes( eiLocal, i )][j]; - std::cout << polygon[i][j] << " "; - } - std::cout << ")"; + 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); } - std::cout << std::endl; - return isPointInPolygon3D(polygon, nV, point); -} +}; bool isPointInsideElement( SurfaceElementSubRegion const & subRegion, arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & referencePosition, @@ -260,13 +197,27 @@ bool isPointInsideElement( SurfaceElementSubRegion const & subRegion, ArrayOfArraysView< localIndex const > const & facesToNodes, real64 const (&location)[3] ) { - std::cout << eiLocal << " " << subRegion.numNodesPerElement( eiLocal ) << std::endl; + typedef std::array< real64, 3 > Point3d; + + // collect element nodes + integer const nV = subRegion.numNodesPerElement( eiLocal ); SurfaceElementSubRegion::NodeMapType const & nodeList = subRegion.nodeList(); - return isPointInsidePolygon( referencePosition, - eiLocal, - 4,//subRegion.numNodesPerElement( eiLocal ), - nodeList, - location ); + 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 ); } /** @@ -528,8 +479,8 @@ void initializeLocalSearch( MeshLevel const & mesh, { mesh.getElemManager().getRegion( targetRegionIndex ).forElementSubRegionsIndex< CellElementSubRegion, FaceElementSubRegion >( [&] ( localIndex const esr, auto const & subRegion ) { - if(targetSubRegionIndex == esr) - subRegion.calculateElementCenters(mesh.getNodeManager().referencePosition()); + if( targetSubRegionIndex == esr ) + subRegion.calculateElementCenters( mesh.getNodeManager().referencePosition()); } ); } @@ -545,15 +496,12 @@ void initializeLocalSearch( MeshLevel const & mesh, 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 ); - std::cout << resElemCenter[targetRegionIndex][targetSubRegionIndex][ei] << " " << dist << std::endl; return dist; } ); // save the index of the reservoir element // note that this reservoir element does not necessarily contains "location" eiInit = ret.second; - - std::cout << "init " << targetRegionIndex << " " << targetSubRegionIndex << " " << eiInit << " " << resElemCenter[targetRegionIndex][targetSubRegionIndex][eiInit] << std::endl; } /** diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index fccdf22d97..25a7b2bd9c 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -433,6 +433,107 @@ bool isPointInsidePolyhedron( arrayView2d< real64 const, nodes::REFERENCE_POSITI return true; } +template< typename POLYGON_TYPE, typename POINT_TYPE > +bool isPointInPolygon2d( POLYGON_TYPE const & polygon, integer n, POINT_TYPE const & point ) +{ + integer count = 0; + + for( integer i = 0; i < n; i++ ) + { + auto const & p1 = polygon[i]; + auto const & p2 = polygon[(i + 1) % n]; + + if((point[1] > std::min( p1[1], p2[1] )) && + (point[1] <= std::max( p1[1], p2[1] )) && + (point[0] <= std::max( p1[0], p2[0] ))) + { + real64 const xIntersect = (point[1] - p1[1]) * (p2[0] - p1[0]) / (p2[1] - p1[1]) + p1[0]; + if( std::abs( p1[0] - p2[0] ) < 1e-10 || point[0] <= xIntersect ) + { + count++; + } + } + } + + return count % 2 == 1; +} + +template< typename POLYGON_TYPE, typename POINT_TYPE > +bool isPointInPolygon3d( POLYGON_TYPE const & polygon, integer const n, POINT_TYPE const & point ) +{ + // 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 d = -(normal[0] * p0[0] + normal[1] * p0[1] + normal[2] * p0[2]); + real64 dist = normal[0] * point[0] + normal[1] * point[1] + normal[2] * point[2] + d; + + if( std::abs( dist ) > 1e-6 ) + { + 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; + for( integer i = 0; i < n; i++ ) + { + projectedPolygon[i][0] = polygon[i][1]; + projectedPolygon[i][1] = polygon[i][2]; + } + 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. * @tparam COORD_TYPE type of coordinate From 7858fd786b5bdfe3936bf8a13ddd5ac4b4eccc77 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 19 Sep 2024 11:18:55 -0500 Subject: [PATCH 03/10] Update ElementSubRegionBase.hpp --- src/coreComponents/mesh/ElementSubRegionBase.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/coreComponents/mesh/ElementSubRegionBase.hpp b/src/coreComponents/mesh/ElementSubRegionBase.hpp index 76c2e115ed..74bd710a31 100644 --- a/src/coreComponents/mesh/ElementSubRegionBase.hpp +++ b/src/coreComponents/mesh/ElementSubRegionBase.hpp @@ -287,7 +287,6 @@ class ElementSubRegionBase : public ObjectManagerBase } LvArray::tensorOps::scale< 3 >( elementCenters[k], 1.0 / numNodes ); - //std::cout << k << " " << elementCenters[k] << std::endl; } ); } }; From 14422f226a8d058301a18ae882ca826ee7ad1c5c Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 19 Sep 2024 11:20:10 -0500 Subject: [PATCH 04/10] Update SurfaceElementRegion.cpp --- src/coreComponents/mesh/SurfaceElementRegion.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index 47f8993b32..78233703bd 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -122,12 +122,8 @@ localIndex SurfaceElementRegion::addToFractureMesh( real64 const time_np1, localIndex const kfe = subRegion.size() - 1; ruptureTime( kfe ) = time_np1; -//std::cout << kfe << " " << elemCenter[ kfe ] << " " << faceCenter[ faceIndices[ 0 ] ] << std::endl; - LvArray::tensorOps::copy< 3 >( elemCenter[ kfe ], faceCenter[ faceIndices[ 0 ] ] ); -//std::cout << kfe << " " << elemCenter[ kfe ] << " " << faceCenter[ faceIndices[ 0 ] ] << std::endl; - faceMap.resizeArray( kfe, 2 ); faceMap[kfe][0] = faceIndices[0]; faceMap[kfe][1] = faceIndices[1]; From 1210866ee18c7e63cf9473f9cd822f5bf5a8f534 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 19 Sep 2024 12:25:09 -0500 Subject: [PATCH 05/10] fix hack --- src/coreComponents/mesh/DomainPartition.cpp | 5 +++-- src/coreComponents/mesh/FaceManager.cpp | 15 ++++++++------- src/coreComponents/mesh/FaceManager.hpp | 4 +++- src/coreComponents/mesh/WellElementSubRegion.cpp | 8 -------- 4 files changed, 14 insertions(+), 18 deletions(-) diff --git a/src/coreComponents/mesh/DomainPartition.cpp b/src/coreComponents/mesh/DomainPartition.cpp index a586b6d9c7..68c137c8fd 100644 --- a/src/coreComponents/mesh/DomainPartition.cpp +++ b/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 )) ) { diff --git a/src/coreComponents/mesh/FaceManager.cpp b/src/coreComponents/mesh/FaceManager.cpp index 9a3923ed35..7e4c25bdda 100644 --- a/src/coreComponents/mesh/FaceManager.cpp +++ b/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() ); diff --git a/src/coreComponents/mesh/FaceManager.hpp b/src/coreComponents/mesh/FaceManager.hpp index 7d7693ab4d..310fdb853c 100644 --- a/src/coreComponents/mesh/FaceManager.hpp +++ b/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. diff --git a/src/coreComponents/mesh/WellElementSubRegion.cpp b/src/coreComponents/mesh/WellElementSubRegion.cpp index c58399cd1e..62f60dfc27 100644 --- a/src/coreComponents/mesh/WellElementSubRegion.cpp +++ b/src/coreComponents/mesh/WellElementSubRegion.cpp @@ -476,14 +476,6 @@ void initializeLocalSearch( MeshLevel const & mesh, localIndex const & targetSubRegionIndex, localIndex & eiInit ) { - { - mesh.getElemManager().getRegion( targetRegionIndex ).forElementSubRegionsIndex< CellElementSubRegion, FaceElementSubRegion >( [&] ( localIndex const esr, auto const & subRegion ) - { - if( targetSubRegionIndex == esr ) - subRegion.calculateElementCenters( mesh.getNodeManager().referencePosition()); - } ); - } - ElementSubRegionBase const & subRegion = mesh.getElemManager().getRegion( targetRegionIndex ).getSubRegion( targetSubRegionIndex ); ElementRegionManager::ElementViewAccessor< arrayView2d< real64 const > > From 4d8a986f094338ba5986dabcdbd86e5b43ca2691 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 19 Sep 2024 14:33:37 -0500 Subject: [PATCH 06/10] wip --- src/coreComponents/mesh/PerforationData.cpp | 37 ++++++++++--------- .../mesh/utilities/ComputationalGeometry.hpp | 9 +++-- 2 files changed, 26 insertions(+), 20 deletions(-) diff --git a/src/coreComponents/mesh/PerforationData.cpp b/src/coreComponents/mesh/PerforationData.cpp index e91970382c..e53e98d31b 100644 --- a/src/coreComponents/mesh/PerforationData.cpp +++ b/src/coreComponents/mesh/PerforationData.cpp @@ -231,23 +231,26 @@ 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 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/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 25a7b2bd9c..cd0d6b117e 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -939,13 +939,15 @@ 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 ) { + std::cout << "getBoundingBox " << elemIndex << std::endl; + // This holds the min coordinates of the set in each direction R1Tensor minCoords = { LvArray::NumericLimits< real64 >::max, LvArray::NumericLimits< real64 >::max, @@ -955,7 +957,8 @@ 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 ) + std::cout << pointIndices[elemIndex].size() << std::endl; + for( localIndex a = 0; a < pointIndices[elemIndex].size(); ++a ) { localIndex const id = pointIndices( elemIndex, a ); for( localIndex d = 0; d < 3; ++d ) From 418f18f409e0f72ad26eefcee1f93d138bcde7ed Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 19 Sep 2024 16:57:08 -0500 Subject: [PATCH 07/10] seems to work --- src/coreComponents/mesh/PerforationData.cpp | 26 ++++++++++++++----- .../mesh/utilities/ComputationalGeometry.hpp | 3 --- 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/src/coreComponents/mesh/PerforationData.cpp b/src/coreComponents/mesh/PerforationData.cpp index e53e98d31b..3879dc4b37 100644 --- a/src/coreComponents/mesh/PerforationData.cpp +++ b/src/coreComponents/mesh/PerforationData.cpp @@ -239,17 +239,31 @@ void PerforationData::getReservoirElementDimensions( MeshLevel const & mesh, // compute the bounding box of the element real64 boxDims[ 3 ]; computationalGeometry::getBoundingBox( ei, - subRegion.nodeList(), - nodeManager.referencePosition(), - boxDims ); + subRegion.nodeList(), + nodeManager.referencePosition(), + boxDims ); // dx and dz from bounding box dx = boxDims[ 0 ]; dy = boxDims[ 1 ]; + dz = boxDims[ 2 ]; - // dz is computed as vol / (dx * dy) - dz = subRegion.getElementVolume()[ei]; - dz /= dx * dy; + 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; + } } ); } diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index cd0d6b117e..31506293bd 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -946,8 +946,6 @@ void getBoundingBox( localIndex const elemIndex, arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & pointCoordinates, VEC_TYPE && boxDims ) { - std::cout << "getBoundingBox " << elemIndex << std::endl; - // This holds the min coordinates of the set in each direction R1Tensor minCoords = { LvArray::NumericLimits< real64 >::max, LvArray::NumericLimits< real64 >::max, @@ -957,7 +955,6 @@ 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 - std::cout << pointIndices[elemIndex].size() << std::endl; for( localIndex a = 0; a < pointIndices[elemIndex].size(); ++a ) { localIndex const id = pointIndices( elemIndex, a ); From cd306951b8eceb51561c02df83eedf6f41a2f50c Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 20 Sep 2024 09:00:18 -0500 Subject: [PATCH 08/10] test case update --- .../singlePhasePoromechanics_FaultModel_well_new_base.xml | 2 +- .../singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml index f6c8fae153..c6bbab3527 100644 --- a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml +++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml @@ -29,7 +29,7 @@ + permeabilityComponents="{ 1.0e-14, 1.0e-14, 1.0e-14 }" /> From 0c6bcef5241f91d8cac4ea11e7c4a9290c038600 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 20 Sep 2024 09:14:28 -0500 Subject: [PATCH 09/10] Update WellElementSubRegion.cpp --- src/coreComponents/mesh/WellElementSubRegion.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/mesh/WellElementSubRegion.cpp b/src/coreComponents/mesh/WellElementSubRegion.cpp index 62f60dfc27..ebaa020f6e 100644 --- a/src/coreComponents/mesh/WellElementSubRegion.cpp +++ b/src/coreComponents/mesh/WellElementSubRegion.cpp @@ -194,7 +194,7 @@ struct PointsEqual bool isPointInsideElement( SurfaceElementSubRegion const & subRegion, arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & referencePosition, localIndex const & eiLocal, - ArrayOfArraysView< localIndex const > const & facesToNodes, + ArrayOfArraysView< localIndex const > const & GEOS_UNUSED_PARAM( facesToNodes ), real64 const (&location)[3] ) { typedef std::array< real64, 3 > Point3d; From d2fc2a63e4dc36e191b8251feec455aa3571ecc6 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 20 Sep 2024 09:17:09 -0500 Subject: [PATCH 10/10] Update ComputationalGeometry.hpp --- .../mesh/utilities/ComputationalGeometry.hpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 31506293bd..3ee5d8d48a 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -433,6 +433,15 @@ 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 ploygon nodes coordinates + * @param[in] n number of polygon nodes + * @param[in] point coordinates of the query point + * @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 ) { @@ -458,6 +467,15 @@ bool isPointInPolygon2d( POLYGON_TYPE const & polygon, integer n, POINT_TYPE con 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 ploygon nodes coordinates + * @param[in] n number of polygon nodes + * @param[in] point coordinates of the query point + * @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 ) {