From e3a68fc62a80a697d0fcf6fcbab84583842ad476 Mon Sep 17 00:00:00 2001 From: Stefano Frambati Date: Mon, 5 Aug 2024 08:25:42 +0200 Subject: [PATCH] fix: robust point location for degenerate configurations (#3202) * implemented method based on winding number * using element global index for element choice in degenerate cases * added global id criterion for points coincident with vertices, as well * using base discretization for precomputing src/rcv * enabled for high orders. Also removed extra high-order nodes where not needed --------- Co-authored-by: Aurelien Co-authored-by: acitrain <60715545+acitrain@users.noreply.github.com> --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 + .../mesh/utilities/ComputationalGeometry.hpp | 397 +++++++++++++++++- .../AcousticFirstOrderWaveEquationSEM.cpp | 45 +- .../AcousticFirstOrderWaveEquationSEM.hpp | 3 +- .../AcousticVTIWaveEquationSEM.cpp | 44 +- .../AcousticVTIWaveEquationSEM.hpp | 3 +- .../AcousticVTIWaveEquationSEMKernel.hpp | 73 ++-- .../isotropic/AcousticWaveEquationSEM.cpp | 43 +- .../isotropic/AcousticWaveEquationSEM.hpp | 3 +- .../ElasticFirstOrderWaveEquationSEM.cpp | 44 +- .../ElasticFirstOrderWaveEquationSEM.hpp | 3 +- .../isotropic/ElasticWaveEquationSEM.cpp | 38 +- .../isotropic/ElasticWaveEquationSEM.hpp | 4 +- .../PrecomputeSourcesAndReceiversKernel.hpp | 195 +++++---- .../wavePropagation/shared/WaveSolverBase.hpp | 3 +- .../shared/WaveSolverUtils.hpp | 8 +- 17 files changed, 692 insertions(+), 220 deletions(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index ed7b9c4e39e..440853cc61d 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3215-6396-9b43cc8 + baseline: integratedTests/baseline_integratedTests-pr3202-6636-2821c93 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index fff3d8f8311..dd686fc7877 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 #3202 (2024-08-03) +====================== +Acoustic VTI tests needed rebaselining after update in source and receiver location algorithm. + PR #3215 (2024-07-23) ====================== Changed the default value for massCreation and name of the wrapper. diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 0227bd7d28b..63c33ffd332 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -381,7 +381,6 @@ int sign( T const val ) return (T( 0 ) < val) - (val < T( 0 )); } - /** * @brief Check if a point is inside a convex polyhedron (3D polygon) * @tparam POINT_TYPE type of @p point @@ -435,6 +434,402 @@ bool isPointInsidePolyhedron( arrayView2d< real64 const, nodes::REFERENCE_POSITI } +/** + * @brief Method to perform lexicographic comparison of two nodes based on coordinates. + * @tparam COORD_TYPE type of coordinate + * @tparam POINT_TYPE type of point + * @param[in] ax x-coordinate of the first vertex + * @param[in] ay y-coordinate of the first vertex + * @param[in] az z-coordinate of the first vertex + * @param[in] bx x-coordinate of the second vertex + * @param[in] by y-coordinate of the second vertex + * @param[in] bz z-coordinate of the second vertex + * @return 0 if the vertices coincide, +1 if a > b and -1 if b > a + */ +template< typename COORD_TYPE, typename POINT_TYPE > +GEOS_HOST_DEVICE +int lexicographicalCompareVertex( POINT_TYPE const ax, POINT_TYPE const ay, POINT_TYPE const az, + COORD_TYPE const bx, COORD_TYPE const by, COORD_TYPE const bz ) +{ + if( ax < bx ) + return -1; + else if( ax > bx ) + return 1; + if( ay < by ) + return -1; + else if( ay > by ) + return 1; + if( az < bz ) + return -1; + else if( az > bz ) + return 1; + return 0; +} + +/** + * @brief Method to perform lexicographic comparison of a node and an edge based on coordinates. + * @tparam COORD_TYPE type of coordinate + * @tparam POINT_TYPE type of point + * @param[in] ax x-coordinate of the first vertex + * @param[in] ay y-coordinate of the first vertex + * @param[in] az z-coordinate of the first vertex + * @param[in] e1x x-coordinate of the first edge vertex + * @param[in] e1y y-coordinate of the first edge vertex + * @param[in] e1z z-coordinate of the first edge vertex + * @param[in] e2x x-coordinate of the second edge vertex + * @param[in] e2y y-coordinate of the second edge vertex + * @param[in] e2z z-coordinate of the second edge vertex + * @return 0 if the vertices lies on the edge, +1 if e > a and -1 if a > e + */ +template< typename COORD_TYPE, typename POINT_TYPE > +GEOS_HOST_DEVICE +int lexicographicalCompareEdge( POINT_TYPE const ax, POINT_TYPE const ay, POINT_TYPE const az, + COORD_TYPE const e1x, COORD_TYPE const e1y, COORD_TYPE const e1z, + COORD_TYPE const e2x, COORD_TYPE const e2y, COORD_TYPE const e2z ) +{ + return lexicographicalCompareVertex( ( e1y - ay ) * ( e2x - ax ), + ( e1z - az ) * ( e2x - ax ), + ( e1z - az ) * ( e2y - ay ), + ( e1x - ax ) * ( e2y - ay ), + ( e1x - ax ) * ( e2z - az ), + ( e1y - ay ) * ( e2z - az ) ); +} + +/** + * @brief Method to perform lexicographic comparison of a node and a triangle based on coordinates. + * @tparam COORD_TYPE type of coordinate + * @tparam POINT_TYPE type of point + * @param[in] ax x-coordinate of the first vertex + * @param[in] ay y-coordinate of the first vertex + * @param[in] az z-coordinate of the first vertex + * @param[in] t1x x-coordinate of the first triangle vertex + * @param[in] t1y y-coordinate of the first triangle vertex + * @param[in] t1z z-coordinate of the first triangle vertex + * @param[in] t2x x-coordinate of the second triangle vertex + * @param[in] t2y y-coordinate of the second triangle vertex + * @param[in] t2z z-coordinate of the second triangle vertex + * @param[in] t3x x-coordinate of the third triangle vertex + * @param[in] t3y y-coordinate of the third triangle vertex + * @param[in] t3z z-coordinate of the third triangle vertex + * @return 0 if the vertices lies on the triangle, +1 if t > a and -1 if t > e + */ +template< typename COORD_TYPE, typename POINT_TYPE > +GEOS_HOST_DEVICE +int lexicographicalCompareTriangle( POINT_TYPE const ax, POINT_TYPE const ay, POINT_TYPE const az, + COORD_TYPE const t1x, COORD_TYPE const t1y, COORD_TYPE const t1z, + COORD_TYPE const t2x, COORD_TYPE const t2y, COORD_TYPE const t2z, + COORD_TYPE const t3x, COORD_TYPE const t3y, COORD_TYPE const t3z ) +{ + COORD_TYPE v1x = t1x - ax; + COORD_TYPE v1y = t1y - ay; + COORD_TYPE v1z = t1z - az; + COORD_TYPE v2x = t2x - ax; + COORD_TYPE v2y = t2y - ay; + COORD_TYPE v2z = t2z - az; + COORD_TYPE v3x = t3x - ax; + COORD_TYPE v3y = t3y - ay; + COORD_TYPE v3z = t3z - az; + COORD_TYPE sign = ( v1x * v2y - v1y * v2x ) * v3z + + ( v2x * v3y - v2y * v3x ) * v1z + + ( v3x * v1y - v3y * v1x ) * v2z; + if( sign > 0 ) + return 1; + else if( sign < 0 ) + return -1; + return 0; +} +/** + * @brief Method to find the reference element touching a vertex. The element with the lowest global ID is chosen + * from the list. + * @param[in] nodeElements the list of elements adjacent to the vertex + * @param[in] elementGlobalIndex the global IDs for elements + * @return element touching the vertex with the least global index + */ +template< typename ... LIST_TYPE > +GEOS_HOST_DEVICE +int findVertexRefElement( arraySlice1d< localIndex const > const & nodeElements, + arrayView1d< globalIndex const > const & elementGlobalIndex ) +{ + localIndex minElement = -1; + globalIndex minElementGID = std::numeric_limits< globalIndex >::max(); + for( int i = 0; i < nodeElements.size(); i++ ) + { + localIndex e = nodeElements( i ); + if( elementGlobalIndex[ e ] < minElementGID ) + { + minElementGID = elementGlobalIndex[ e ]; + minElement = e; + } + } + return minElement; +} + +/** + * @brief Method to find the reference element for an edge. The element with the lowest global ID is chosen + * from the list. + * @param[in] nodeElements1 the list of elements adjacent to the first node + * @param[in] nodeElements2 the list of elements adjacent to the second node + * @param[in] elementGlobalIndex the global IDs for elements + * @return the element shared by the two nodes, with the minimal global index + */ +template< typename ... LIST_TYPE > +GEOS_HOST_DEVICE +int findEdgeRefElement( arraySlice1d< localIndex const > const & nodeElements1, + arraySlice1d< localIndex const > const & nodeElements2, + arrayView1d< globalIndex const > const & elementGlobalIndex ) +{ + localIndex minElement = -1; + globalIndex minElementGID = std::numeric_limits< globalIndex >::max(); + for( int i = 0; i < nodeElements1.size(); i++ ) + { + localIndex e1 = nodeElements1( i ); + for( int j = 0; j < nodeElements2.size(); j++ ) + { + localIndex e2 = nodeElements2( j ); + if( e1 == e2 ) + { + if( elementGlobalIndex[ e1 ] < minElementGID ) + { + minElementGID = elementGlobalIndex[ e1 ]; + minElement = e1; + } + } + } + } + return minElement; +} + +/** + * @brief Method to find the reference element for a triangle. The element with the lowest global ID is chosen + * from the list. + * @param[in] nodeElements1 the list of elements adjacent to the first node + * @param[in] nodeElements2 the list of elements adjacent to the second node + * @param[in] nodeElements3 the list of elements adjacent to the third node + * @param[in] elementGlobalIndex the global IDs for elements + * @return the element shared by the three nodes, with the minimal global index + */ +template< typename ... LIST_TYPE > +GEOS_HOST_DEVICE +int findTriangleRefElement( arraySlice1d< localIndex const > const & nodeElements1, + arraySlice1d< localIndex const > const & nodeElements2, + arraySlice1d< localIndex const > const & nodeElements3, + arrayView1d< globalIndex const > const & elementGlobalIndex ) +{ + localIndex minElement = -1; + globalIndex minElementGID = std::numeric_limits< globalIndex >::max(); + for( int i = 0; i < nodeElements1.size(); i++ ) + { + localIndex e1 = nodeElements1( i ); + for( int j = 0; j < nodeElements2.size(); j++ ) + { + localIndex e2 = nodeElements2( j ); + for( int k = 0; k < nodeElements3.size(); k++ ) + { + localIndex e3 = nodeElements3( k ); + if( e1 == e2 && e2 == e3 ) + { + if( elementGlobalIndex[ e1 ] < minElementGID ) + { + minElementGID = elementGlobalIndex[ e1 ]; + minElement = e1; + } + } + } + } + } + return minElement; +} + +/** + * @brief Computes the winding number of a point with respecto 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 + * @param[in] elementsToFaces map from elements to faces + * @param[in] facesToNodes map from faces to nodes + * @param[in] nodesToElements map from nodes to elements + * @param[in] nodeLocalToGlobal global indices of nodes + * @param[in] elementLocalToGlobal global indices of elements + * @param[in] elemCenter coordinates of the element centroid + * @param[in] point coordinates of the query point + * @return the signed winding number, which is positive if and only if the point is inside the mesh element. + */ +template< typename COORD_TYPE, typename POINT_TYPE > +GEOS_HOST_DEVICE +bool computeWindingNumber( localIndex element, + arrayView2d< COORD_TYPE const, nodes::REFERENCE_POSITION_USD > const & nodeCoordinates, + arrayView2d< localIndex const > const & elementsToFaces, + ArrayOfArraysView< localIndex const > const & facesToNodes, + ArrayOfArraysView< localIndex const > const & nodesToElements, + arrayView1d< globalIndex const > const & nodeLocalToGlobal, + arrayView1d< globalIndex const > const & elementLocalToGlobal, + POINT_TYPE const & elemCenter, + POINT_TYPE const & point ) +{ + arraySlice1d< localIndex const > const & faceIndices = elementsToFaces[ element ]; + localIndex const numFaces = faceIndices.size(); + int omega = 0; + for( localIndex kf = 0; kf < numFaces; ++kf ) + { + // triangulate the face. The triangulation must be done in a consistent way across ranks. + // This can be achieved by always picking the vertex with the lowest global index as root. + localIndex const faceIndex = faceIndices[kf]; + globalIndex minGlobalId = std::numeric_limits< globalIndex >::max(); + localIndex minVertex = -1; + localIndex numFaceVertices = facesToNodes[faceIndex].size(); + for( localIndex v = 0; v < numFaceVertices; v++ ) + { + localIndex vIndex = facesToNodes( faceIndex, v ); + globalIndex globalId = nodeLocalToGlobal[ vIndex ]; + if( globalId < minGlobalId ) + { + minGlobalId = globalId; + minVertex = vIndex; + } + } + // triangulate the face using the minimum-id vertex as root + localIndex vi[ 3 ] = { minVertex, -1, -1 }; + for( localIndex v = 0; v < numFaceVertices; v++ ) + { + vi[ 1 ] = facesToNodes( faceIndex, v ); + vi[ 2 ] = facesToNodes( faceIndex, (v + 1) % numFaceVertices ); + if( vi[ 1 ] != minVertex && vi[ 2 ] != minVertex ) + { + // To make the algorithm independent of rank, always take the two additional vertices in increasing global ID + if( nodeLocalToGlobal[ vi[ 1 ] ] > nodeLocalToGlobal[ vi[ 2 ] ] ) + { + localIndex temp = vi[ 1 ]; + vi[ 1 ] = vi[ 2 ]; + vi[ 2 ] = temp; + } + COORD_TYPE v1x = nodeCoordinates( vi[ 0 ], 0 ); + COORD_TYPE v1y = nodeCoordinates( vi[ 0 ], 1 ); + COORD_TYPE v1z = nodeCoordinates( vi[ 0 ], 2 ); + COORD_TYPE v2x = nodeCoordinates( vi[ 1 ], 0 ); + COORD_TYPE v2y = nodeCoordinates( vi[ 1 ], 1 ); + COORD_TYPE v2z = nodeCoordinates( vi[ 1 ], 2 ); + COORD_TYPE v3x = nodeCoordinates( vi[ 2 ], 0 ); + COORD_TYPE v3y = nodeCoordinates( vi[ 2 ], 1 ); + COORD_TYPE v3z = nodeCoordinates( vi[ 2 ], 2 ); + // check the orientation of this triangle + R1Tensor vv1 = { v2x - v1x, v2y - v1y, v2z - v1z }; + R1Tensor vv2 = { v3x - v1x, v3y - v1y, v3z - v1z }; + R1Tensor dist = { elemCenter[ 0 ] - ( v1x + v2x + v3x )/3.0, + elemCenter[ 1 ] - ( v1y + v2y + v3y )/3.0, + elemCenter[ 2 ] - ( v1z + v2z + v3z )/3.0 }; + R1Tensor norm = { }; + LvArray::tensorOps::crossProduct( norm, vv1, vv2 ); + // check if face is oriented coherently, and change sign otherwise + int sign = LvArray::tensorOps::AiBi< 3 >( norm, dist ) > 0 ? -1 : +1; + // Compute the winding number contributed by this triangle + int cmp1 = lexicographicalCompareVertex( point[ 0 ], point[ 1 ], point[ 2 ], v1x, v1y, v1z ); + if( cmp1 == 0 ) + { + return findVertexRefElement( nodesToElements[ vi[ 0 ] ], elementLocalToGlobal ) == element; + } + int cmp2 = lexicographicalCompareVertex( point[ 0 ], point[ 1 ], point[ 2 ], v2x, v2y, v2z ); + if( cmp2 == 0 ) + { + return findVertexRefElement( nodesToElements[ vi[ 1 ] ], elementLocalToGlobal ) == element; + } + int cmp3 = lexicographicalCompareVertex( point[ 0 ], point[ 1 ], point[ 2 ], v3x, v3y, v3z ); + if( cmp3 == 0 ) + { + return findVertexRefElement( nodesToElements[ vi[ 2 ] ], elementLocalToGlobal ) == element; + } + int facecmp = 0; + int edgecmp = 0; + if( cmp1 != cmp2 ) + { + edgecmp = lexicographicalCompareEdge( point[ 0 ], point[ 1 ], point[ 2 ], + v1x, v1y, v1z, + v2x, v2y, v2z ); + if( edgecmp == 0 ) + { + return findEdgeRefElement( nodesToElements[ vi[ 0 ] ], nodesToElements[ vi[ 1 ] ], elementLocalToGlobal ) == element; + } + facecmp += sign * edgecmp; + } + if( cmp2 != cmp3 ) + { + edgecmp = lexicographicalCompareEdge( point[ 0 ], point[ 1 ], point[ 2 ], + v2x, v2y, v2z, + v3x, v3y, v3z ); + if( edgecmp == 0 ) + { + return findEdgeRefElement( nodesToElements[ vi[ 1 ] ], nodesToElements[ vi[ 2 ] ], elementLocalToGlobal ) == element; + } + facecmp += sign * edgecmp; + } + if( cmp3 != cmp1 ) + { + edgecmp = lexicographicalCompareEdge( point[ 0 ], point[ 1 ], point[ 2 ], + v3x, v3y, v3z, + v1x, v1y, v1z ); + if( edgecmp == 0 ) + { + return findEdgeRefElement( nodesToElements[ vi[ 0 ] ], nodesToElements[ vi[ 2 ] ], elementLocalToGlobal ) == element; + } + facecmp += sign * edgecmp; + } + // if all edges are on the same side, this triangle does not contribute to the winding number + if( facecmp == 0 ) + continue; + facecmp = lexicographicalCompareTriangle( point[ 0 ], point[ 1 ], point[ 2 ], + v1x, v1y, v1z, + v2x, v2y, v2z, + v3x, v3y, v3z ); + + if( facecmp == 0 ) + { + return findTriangleRefElement( nodesToElements[ vi[ 0 ] ], nodesToElements[ vi[ 1 ] ], nodesToElements[ vi[ 2 ] ], elementLocalToGlobal ) == element; + } + omega += sign * facecmp; + } + } + } + + return omega; +} + +/** + * @brief Check if a point is inside a convex polyhedron (3D polygon), using a robust method + * to avoid ambiguity when the point lies on an interface. + * This method is based on the following method: + * - the winding number omega of the point with respect to the cell is used to determine if the point is inside (omega=1) or not (omega=0) + * - corner cases (point lying on a face, edge or vertex of the cell) are detected using a robust method based on lexicographical + * comparisons + * - these comparisons are made consistent across MPI ranks by consistently arranging items based on global indices (GIDs). In particular: + * - Faces are triangulated using the vertex with the smallest GID as root; + * - Edges and faces are described by vertices in increasing GID order + * - Finally, if the point lies on a (vertex, edge, face), it is assigned to the first shared element with the least global index + * @tparam POINT_TYPE type of @p point + * @param[in] element the element to be checked + * @param[in] nodeCoordinates a global array of nodal coordinates + * @param[in] elementsToFaces map from elements to faces + * @param[in] facesToNodes map from faces to nodes + * @param[in] nodesToElements map from nodes to elements + * @param[in] nodeLocalToGlobal global indices of nodes + * @param[in] elementLocalToGlobal global indices of elements + * @param[in] elemCenter coordinates of the element centroid + * @param[in] point coordinates of the query point + * @return whether the point is inside + */ +template< typename COORD_TYPE, typename POINT_TYPE > +GEOS_HOST_DEVICE +bool isPointInsideConvexPolyhedronRobust( localIndex element, + arrayView2d< COORD_TYPE const, nodes::REFERENCE_POSITION_USD > const & nodeCoordinates, + arrayView2d< localIndex const > const & elementsToFaces, + ArrayOfArraysView< localIndex const > const & facesToNodes, + ArrayOfArraysView< localIndex const > const & nodesToElements, + arrayView1d< globalIndex const > const & nodeLocalToGlobal, + arrayView1d< globalIndex const > const & elementLocalToGlobal, + POINT_TYPE const & elemCenter, + POINT_TYPE const & point ) +{ + return computeWindingNumber( element, nodeCoordinates, elementsToFaces, facesToNodes, nodesToElements, nodeLocalToGlobal, elementLocalToGlobal, elemCenter, point ) > 0; +} + /** * @brief Compute the dimensions of the bounding box containing the element * defined here by the coordinates of its vertices. diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.cpp index 32755474615..e9db0ca9617 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.cpp @@ -153,15 +153,14 @@ void AcousticFirstOrderWaveEquationSEM::postInputInitialization() m_receiverRegion.resize( numReceiversGlobal ); } -void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - NodeManager const & nodeManager = mesh.getNodeManager(); - FaceManager const & faceManager = mesh.getFaceManager(); + GEOS_MARK_FUNCTION; - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const - X = nodeManager.getField< fields::referencePosition32 >().toViewConst(); - arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - arrayView2d< real64 const > const faceCenter = faceManager.faceCenter(); + arrayView1d< globalIndex const > const nodeLocalToGlobal = baseMesh.getNodeManager().localToGlobalMap().toViewConst(); + ArrayOfArraysView< localIndex const > const nodesToElements = baseMesh.getNodeManager().elementList().toViewConst(); + ArrayOfArraysView< localIndex const > const facesToNodes = baseMesh.getFaceManager().nodeList().toViewConst(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodeCoords = baseMesh.getNodeManager().referencePosition(); arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst(); arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView(); @@ -195,8 +194,11 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev } } - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const regionIndex, - CellElementSubRegion & elementSubRegion ) + mesh.getElemManager().forElementSubRegionsComplete< CellElementSubRegion >( regionNames, [&]( localIndex const, + localIndex const regionIndex, + localIndex const esr, + ElementRegionBase &, + CellElementSubRegion & elementSubRegion ) { GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Hexahedron, getDataContext() << ": Invalid type of element, the acoustic solver is designed for hexahedral meshes only (C3D8) ", @@ -204,8 +206,10 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes = baseMesh.getElemManager().getRegion( regionIndex ).getSubRegion< CellElementSubRegion >( esr ).nodeList(); arrayView2d< real64 const > const elemCenter = elementSubRegion.getElementCenter(); arrayView1d< integer const > const elemGhostRank = elementSubRegion.ghostRank(); + arrayView1d< globalIndex const > const elemLocalToGlobal = elementSubRegion.localToGlobalMap().toViewConst(); finiteElement::FiniteElementBase const & fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); @@ -213,21 +217,21 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev { using FE_TYPE = TYPEOFREF( finiteElement ); - localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); - PreComputeSourcesAndReceivers:: Compute1DSourceAndReceiverConstantsWithElementsAndRegionStorage < EXEC_POLICY, FE_TYPE > ( elementSubRegion.size(), regionIndex, - numFacesPerElem, - X, + facesToNodes, + nodeCoords, + nodeLocalToGlobal, + elemLocalToGlobal, + nodesToElements, + baseElemsToNodes, elemGhostRank, elemsToNodes, elemsToFaces, elemCenter, - faceNormal, - faceCenter, sourceCoordinates, sourceIsAccessible, sourceElem, @@ -258,11 +262,12 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro applyFreeSurfaceBC( 0.0, domain ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - precomputeSourceAndReceiverTerm( mesh, regionNames ); + MeshLevel & baseMesh = domain.getMeshBodies().getGroup< MeshBody >( meshBodyName ).getBaseDiscretization(); + precomputeSourceAndReceiverTerm( baseMesh, mesh, regionNames ); NodeManager & nodeManager = mesh.getNodeManager(); FaceManager & faceManager = mesh.getFaceManager(); @@ -435,7 +440,7 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t { NodeManager & nodeManager = mesh.getNodeManager(); - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.getField< fields::referencePosition32 >().toViewConst(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst(); arrayView1d< real32 const > const mass = nodeManager.getField< acousticfields::AcousticMassVector >(); arrayView1d< real32 const > const damping = nodeManager.getField< acousticfields::DampingVector >(); @@ -468,7 +473,7 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t VelocityComputation< FE_TYPE > kernel( finiteElement ); kernel.template launch< EXEC_POLICY, ATOMIC_POLICY > ( elementSubRegion.size(), - X, + nodeCoords, elemsToNodes, p_np1, density, @@ -482,7 +487,7 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t ( elementSubRegion.size(), regionIndex, nodeManager.size(), - X, + nodeCoords, elemsToNodes, velocity_x, velocity_y, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.hpp index 423dc8ae5ba..52e72af8d91 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/firstOrderEqn/isotropic/AcousticFirstOrderWaveEquationSEM.hpp @@ -119,9 +119,10 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase /** * @brief Locate sources and receivers position in the mesh elements, evaluate the basis functions at each point and save them to the * corresponding elements nodes. + * @param baseMesh the level-0 mesh * @param mesh mesh of the computational domain */ - virtual void precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + virtual void precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** * @brief Apply free surface condition to the face define in the geometry box from the xml diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp index 23cd9ba2db4..1bfe1f970d8 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.cpp @@ -112,18 +112,15 @@ void AcousticVTIWaveEquationSEM::postInputInitialization() m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 ); } -void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, +void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - NodeManager const & nodeManager = mesh.getNodeManager(); - FaceManager const & faceManager = mesh.getFaceManager(); - - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords32 = - nodeManager.getField< fields::referencePosition32 >().toViewConst(); - - arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - arrayView2d< real64 const > const faceCenter = faceManager.faceCenter(); + GEOS_MARK_FUNCTION; + arrayView1d< globalIndex const > const nodeLocalToGlobal = baseMesh.getNodeManager().localToGlobalMap().toViewConst(); + ArrayOfArraysView< localIndex const > const nodesToElements = baseMesh.getNodeManager().elementList().toViewConst(); + ArrayOfArraysView< localIndex const > const facesToNodes = baseMesh.getFaceManager().nodeList().toViewConst(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodeCoords = baseMesh.getNodeManager().referencePosition().toViewConst(); arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst(); arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView(); @@ -153,8 +150,11 @@ void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & me } } - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - CellElementSubRegion & elementSubRegion ) + mesh.getElemManager().forElementSubRegionsComplete< CellElementSubRegion >( regionNames, [&]( localIndex const, + localIndex const er, + localIndex const esr, + ElementRegionBase &, + CellElementSubRegion & elementSubRegion ) { GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Hexahedron, "Invalid type of element, the acoustic solver is designed for hexahedral meshes only (C3D8), using the SEM formulation", @@ -162,8 +162,10 @@ void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & me arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes = baseMesh.getElemManager().getRegion( er ).getSubRegion< CellElementSubRegion >( esr ).nodeList(); arrayView2d< real64 const > const elemCenter = elementSubRegion.getElementCenter(); arrayView1d< integer const > const elemGhostRank = elementSubRegion.ghostRank(); + arrayView1d< globalIndex const > const elemLocalToGlobal = elementSubRegion.localToGlobalMap().toViewConst(); finiteElement::FiniteElementBase const & fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); @@ -171,20 +173,20 @@ void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & me { using FE_TYPE = TYPEOFREF( finiteElement ); - localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); - acousticVTIWaveEquationSEMKernels:: PrecomputeSourceAndReceiverKernel:: launch< EXEC_POLICY, FE_TYPE > ( elementSubRegion.size(), - numFacesPerElem, - nodeCoords32, + facesToNodes, + nodeCoords, + nodeLocalToGlobal, + elemLocalToGlobal, + nodesToElements, + baseElemsToNodes, elemGhostRank, elemsToNodes, elemsToFaces, elemCenter, - faceNormal, - faceCenter, sourceCoordinates, sourceIsAccessible, sourceNodeIds, @@ -233,19 +235,19 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups() applyFreeSurfaceBC( 0.0, domain ); precomputeSurfaceFieldIndicator( domain ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - precomputeSourceAndReceiverTerm( mesh, regionNames ); + MeshLevel & baseMesh = domain.getMeshBodies().getGroup< MeshBody >( meshBodyName ).getBaseDiscretization(); + precomputeSourceAndReceiverTerm( baseMesh, mesh, regionNames ); NodeManager & nodeManager = mesh.getNodeManager(); FaceManager & faceManager = mesh.getFaceManager(); /// get the array of indicators: 1 if the face is on the boundary; 0 otherwise arrayView1d< integer > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator(); - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = - nodeManager.getField< fields::referencePosition32 >().toViewConst(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst(); /// get face to node map ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst(); diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp index 889cb424261..d672c398387 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEM.hpp @@ -127,9 +127,10 @@ class AcousticVTIWaveEquationSEM : public WaveSolverBase /** * @brief Locate sources and receivers position in the mesh elements, evaluate the basis functions at each point and save them to the * corresponding elements nodes. + * @param baseMesh the level-0 mesh * @param mesh mesh of the computational domain */ - virtual void precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + virtual void precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** * @brief Compute the lateral and bottom surface Field indicators of the boxed domain diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp index 1acd21ac28b..e4de8b7d4ed 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIWaveEquationSEMKernel.hpp @@ -40,38 +40,43 @@ struct PrecomputeSourceAndReceiverKernel * @tparam EXEC_POLICY execution policy * @tparam FE_TYPE finite element type * @param[in] size the number of cells in the subRegion - * @param[in] numFacesPerElem number of faces per element - * @param[in] nodeCoords coordinates of the nodes - * @param[in] elemGhostRank the ghost ranks + * @param[in] baseFacesToNodes face to node map + * @param[in] baseNodeCoords coordinates of the nodes + * @param[in] baseNodeLocalToGlobal local to global index map for nodes + * @param[in] elementLocalToGlobal local to global index map for elements + * @param[in] baseNodesToElements node to element map for the base mesh + * @param[in] baseElemsToNodes element to node map for the base mesh + * @param[in] elemGhostRank rank of the ghost element * @param[in] elemsToNodes map from element to nodes * @param[in] elemsToFaces map from element to faces - * @param[in] facesToNodes map from faces to nodes * @param[in] elemCenter coordinates of the element centers * @param[in] sourceCoordinates coordinates of the source terms * @param[out] sourceIsAccessible flag indicating whether the source is accessible or not * @param[out] sourceNodeIds indices of the nodes of the element where the source is located - * @param[out] sourceNodeConstants constant part of the source terms + * @param[out] sourceConstants constant part of the source terms * @param[in] receiverCoordinates coordinates of the receiver terms * @param[out] receiverIsLocal flag indicating whether the receiver is local or not * @param[out] receiverNodeIds indices of the nodes of the element where the receiver is located - * @param[out] receiverNodeConstants constant part of the receiver term - * @param[out] sourceValue the value of the source - * @param[in] dt the time step size - * @param[in] timeSourceFrequency the time frequency of the source + * @param[out] receiverConstants constant part of the receiver term + * @param[out] sourceValue value of the temporal source (eg. Ricker) + * @param[in] dt time-step + * @param[in] timeSourceFrequency the central frequency of the source * @param[in] timeSourceDelay the time delay of the source - * @param[in] rickerOrder the order of the ricker + * @param[in] rickerOrder order of the Ricker wavelet */ template< typename EXEC_POLICY, typename FE_TYPE > static void launch( localIndex const size, - localIndex const numFacesPerElem, - arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + ArrayOfArraysView< localIndex const > const baseFacesToNodes, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const baseNodeCoords, + arrayView1d< globalIndex const > const baseNodeLocalToGlobal, + arrayView1d< globalIndex const > const elementLocalToGlobal, + ArrayOfArraysView< localIndex const > const baseNodesToElements, + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes, arrayView1d< integer const > const elemGhostRank, arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes, arrayView2d< localIndex const > const elemsToFaces, arrayView2d< real64 const > const & elemCenter, - arrayView2d< real64 const > const faceNormal, - arrayView2d< real64 const > const faceCenter, arrayView2d< real64 const > const sourceCoordinates, arrayView1d< localIndex > const sourceIsAccessible, arrayView2d< localIndex > const sourceNodeIds, @@ -106,20 +111,23 @@ struct PrecomputeSourceAndReceiverKernel sourceCoordinates[isrc][2] }; bool const sourceFound = - WaveSolverUtils::locateSourceElement( numFacesPerElem, - center, - faceNormal, - faceCenter, - elemsToFaces[k], - coords ); + computationalGeometry::isPointInsideConvexPolyhedronRobust( k, + baseNodeCoords, + elemsToFaces, + baseFacesToNodes, + baseNodesToElements, + baseNodeLocalToGlobal, + elementLocalToGlobal, + center, + coords ); if( sourceFound ) { real64 coordsOnRefElem[3]{}; WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, - elemsToNodes[k], - nodeCoords, + baseElemsToNodes[k], + baseNodeCoords, coordsOnRefElem ); sourceIsAccessible[isrc] = 1; @@ -153,19 +161,22 @@ struct PrecomputeSourceAndReceiverKernel receiverCoordinates[ircv][2] }; real64 coordsOnRefElem[3]{}; - bool const receiverFound = - WaveSolverUtils::locateSourceElement( numFacesPerElem, - center, - faceNormal, - faceCenter, - elemsToFaces[k], - coords ); + bool const receiverFound = + computationalGeometry::isPointInsideConvexPolyhedronRobust( k, + baseNodeCoords, + elemsToFaces, + baseFacesToNodes, + baseNodesToElements, + baseNodeLocalToGlobal, + elementLocalToGlobal, + center, + coords ); if( receiverFound && elemGhostRank[k] < 0 ) { WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, - elemsToNodes[k], - nodeCoords, + baseElemsToNodes[k], + baseNodeCoords, coordsOnRefElem ); receiverIsLocal[ircv] = 1; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp index 0d5c2a7a354..77e3c55098a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp @@ -118,18 +118,15 @@ void AcousticWaveEquationSEM::postInputInitialization() m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, m_receiverCoordinates.size( 0 ) + 1 ); } -void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, +void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { GEOS_MARK_FUNCTION; - NodeManager const & nodeManager = mesh.getNodeManager(); - FaceManager const & faceManager = mesh.getFaceManager(); - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const - nodeCoords32 = nodeManager.getField< fields::referencePosition32 >().toViewConst(); - - arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - arrayView2d< real64 const > const faceCenter = faceManager.faceCenter(); + arrayView1d< globalIndex const > const nodeLocalToGlobalMap = baseMesh.getNodeManager().localToGlobalMap().toViewConst(); + ArrayOfArraysView< localIndex const > const nodesToElements = baseMesh.getNodeManager().elementList().toViewConst(); + ArrayOfArraysView< localIndex const > const facesToNodes = baseMesh.getFaceManager().nodeList().toViewConst(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodeCoords = baseMesh.getNodeManager().referencePosition(); arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst(); arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView(); @@ -159,8 +156,11 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, } } - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - CellElementSubRegion & elementSubRegion ) + mesh.getElemManager().forElementSubRegionsComplete< CellElementSubRegion >( regionNames, [&]( localIndex const, + localIndex const er, + localIndex const esr, + ElementRegionBase &, + CellElementSubRegion & elementSubRegion ) { GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Hexahedron, getDataContext() << ": Invalid type of element, the acoustic solver is designed for hexahedral meshes only (C3D8), using the SEM formulation", @@ -168,8 +168,10 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes = baseMesh.getElemManager().getRegion( er ).getSubRegion< CellElementSubRegion >( esr ).nodeList(); arrayView2d< real64 const > const elemCenter = elementSubRegion.getElementCenter(); arrayView1d< integer const > const elemGhostRank = elementSubRegion.ghostRank(); + arrayView1d< globalIndex const > const elemLocalToGlobalMap = elementSubRegion.localToGlobalMap().toViewConst(); finiteElement::FiniteElementBase const & fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); @@ -177,22 +179,22 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, { using FE_TYPE = TYPEOFREF( finiteElement ); - localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); - { GEOS_MARK_SCOPE( acousticWaveEquationSEMKernels::PrecomputeSourceAndReceiverKernel ); PreComputeSourcesAndReceivers:: Compute1DSourceAndReceiverConstants < EXEC_POLICY, FE_TYPE > ( elementSubRegion.size(), - numFacesPerElem, - nodeCoords32, + facesToNodes, + nodeCoords, + nodeLocalToGlobalMap, + elemLocalToGlobalMap, + nodesToElements, + baseElemsToNodes, elemGhostRank, elemsToNodes, elemsToFaces, elemCenter, - faceNormal, - faceCenter, sourceCoordinates, sourceIsAccessible, sourceNodeIds, @@ -250,11 +252,14 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() applyFreeSurfaceBC( 0.0, domain ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - precomputeSourceAndReceiverTerm( mesh, regionNames ); + MeshLevel & baseMesh = domain.getMeshBodies().getGroup< MeshBody >( meshBodyName ).getBaseDiscretization(); + precomputeSourceAndReceiverTerm( baseMesh, mesh, regionNames ); NodeManager & nodeManager = mesh.getNodeManager(); FaceManager & faceManager = mesh.getFaceManager(); @@ -430,7 +435,7 @@ void AcousticWaveEquationSEM::initializePML() NodeManager & nodeManager = mesh.getNodeManager(); /// WARNING: the array below is one of the PML auxiliary variables arrayView1d< real32 > const indicatorPML = nodeManager.getField< acousticfields::AuxiliaryVar4PML >(); - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords32 = nodeManager.getField< fields::referencePosition32 >().toViewConst(); + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords32 = nodeManager.getField< fields::referencePosition32 >().toViewConst(); indicatorPML.zero(); real32 xInteriorMin[3]{}; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp index da36ea8dcbb..0e0729bfbd9 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp @@ -148,9 +148,10 @@ class AcousticWaveEquationSEM : public WaveSolverBase /** * @brief Locate sources and receivers position in the mesh elements, evaluate the basis functions at each point and save them to the * corresponding elements nodes. + * @param baseMesh the level-0 mesh * @param mesh mesh of the computational domain */ - virtual void precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + virtual void precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** * @brief Apply free surface condition to the face define in the geometry box from the xml diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.cpp index 07eb4802edd..98d04bc3b7a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.cpp @@ -202,14 +202,14 @@ void ElasticFirstOrderWaveEquationSEM::postInputInitialization() } -void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - NodeManager const & nodeManager = mesh.getNodeManager(); - FaceManager const & faceManager = mesh.getFaceManager(); + GEOS_MARK_FUNCTION; - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.getField< fields::referencePosition32 >().toViewConst(); - arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - arrayView2d< real64 const > const faceCenter = faceManager.faceCenter(); + arrayView1d< globalIndex const > const nodeLocalToGlobal = baseMesh.getNodeManager().localToGlobalMap().toViewConst(); + ArrayOfArraysView< localIndex const > const nodesToElements = baseMesh.getNodeManager().elementList().toViewConst(); + ArrayOfArraysView< localIndex const > const facesToNodes = baseMesh.getFaceManager().nodeList().toViewConst(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodeCoords = baseMesh.getNodeManager().referencePosition(); arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst(); arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView(); @@ -245,8 +245,11 @@ void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLeve } } - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const regionIndex, - CellElementSubRegion & elementSubRegion ) + mesh.getElemManager().forElementSubRegionsComplete< CellElementSubRegion >( regionNames, [&]( localIndex const, + localIndex const regionIndex, + localIndex const esr, + ElementRegionBase &, + CellElementSubRegion & elementSubRegion ) { GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Hexahedron, @@ -255,8 +258,10 @@ void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLeve arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes = baseMesh.getElemManager().getRegion( regionIndex ).getSubRegion< CellElementSubRegion >( esr ).nodeList(); arrayView2d< real64 const > const elemCenter = elementSubRegion.getElementCenter(); arrayView1d< integer const > const elemGhostRank = elementSubRegion.ghostRank(); + arrayView1d< globalIndex const > const elemLocalToGlobal = elementSubRegion.localToGlobalMap().toViewConst(); finiteElement::FiniteElementBase const & fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); @@ -264,21 +269,21 @@ void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLeve { using FE_TYPE = TYPEOFREF( finiteElement ); - localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); - PreComputeSourcesAndReceivers:: Compute1DSourceAndReceiverConstantsWithElementsAndRegionStorage < EXEC_POLICY, FE_TYPE > ( elementSubRegion.size(), regionIndex, - numFacesPerElem, - X, + facesToNodes, + nodeCoords, + nodeLocalToGlobal, + elemLocalToGlobal, + nodesToElements, + baseElemsToNodes, elemGhostRank, elemsToNodes, elemsToFaces, elemCenter, - faceNormal, - faceCenter, sourceCoordinates, sourceIsAccessible, sourceElem, @@ -310,11 +315,12 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou real64 const time = 0.0; applyFreeSurfaceBC( time, domain ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - precomputeSourceAndReceiverTerm( mesh, regionNames ); + MeshLevel & baseMesh = domain.getMeshBodies().getGroup< MeshBody >( meshBodyName ).getBaseDiscretization(); + precomputeSourceAndReceiverTerm( baseMesh, mesh, regionNames ); NodeManager & nodeManager = mesh.getNodeManager(); FaceManager & faceManager = mesh.getFaceManager(); @@ -500,7 +506,7 @@ real64 ElasticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & ti NodeManager & nodeManager = mesh.getNodeManager(); - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.getField< fields::referencePosition32 >().toViewConst(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst(); arrayView1d< real32 const > const mass = nodeManager.getField< elasticfields::ElasticMassVector >(); arrayView1d< real32 > const dampingx = nodeManager.getField< elasticfields::DampingVectorx >(); @@ -549,7 +555,7 @@ real64 ElasticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & ti kernel.template launch< EXEC_POLICY, ATOMIC_POLICY > ( elementSubRegion.size(), regionIndex, - X, + nodeCoords, elemsToNodes, ux_np1, uy_np1, @@ -578,7 +584,7 @@ real64 ElasticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & ti kernel2.template launch< EXEC_POLICY, ATOMIC_POLICY > ( elementSubRegion.size(), nodeManager.size(), - X, + nodeCoords, elemsToNodes, stressxx, stressyy, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.hpp index 4ab8c383b0b..575226716f9 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/firstOrderEqn/isotropic/ElasticFirstOrderWaveEquationSEM.hpp @@ -125,10 +125,11 @@ class ElasticFirstOrderWaveEquationSEM : public WaveSolverBase /** * @brief Locate sources and receivers position in the mesh elements, evaluate the basis functions at each point and save them to the * corresponding elements nodes. + * @param baseMesh the level-0 mesh * @param mesh mesh of the computational domain * @param regionNames name of the region you are currently on */ - virtual void precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + virtual void precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** * @brief Apply free surface condition to the face define in the geometry box from the xml diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp index 6a2d0d0d7d4..234915bc7f0 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp @@ -223,15 +223,14 @@ void ElasticWaveEquationSEM::postInputInitialization() } -void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - NodeManager const & nodeManager = mesh.getNodeManager(); - FaceManager const & faceManager = mesh.getFaceManager(); + GEOS_MARK_FUNCTION; - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X = - nodeManager.getField< fields::referencePosition32 >().toViewConst(); - arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - arrayView2d< real64 const > const faceCenter = faceManager.faceCenter(); + arrayView1d< globalIndex const > const nodeLocalToGlobal = baseMesh.getNodeManager().localToGlobalMap().toViewConst(); + ArrayOfArraysView< localIndex const > const nodesToElements = baseMesh.getNodeManager().elementList().toViewConst(); + ArrayOfArraysView< localIndex const > const facesToNodes = baseMesh.getFaceManager().nodeList().toViewConst(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodeCoords = baseMesh.getNodeManager().referencePosition(); arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst(); arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView(); @@ -267,8 +266,11 @@ void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, } } - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - CellElementSubRegion & elementSubRegion ) + mesh.getElemManager().forElementSubRegionsComplete< CellElementSubRegion >( regionNames, [&]( localIndex const, + localIndex const er, + localIndex const esr, + ElementRegionBase &, + CellElementSubRegion & elementSubRegion ) { GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Hexahedron, @@ -277,8 +279,10 @@ void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes = baseMesh.getElemManager().getRegion( er ).getSubRegion< CellElementSubRegion >( esr ).nodeList(); arrayView2d< real64 const > const elemCenter = elementSubRegion.getElementCenter(); arrayView1d< integer const > const elemGhostRank = elementSubRegion.ghostRank(); + arrayView1d< globalIndex const > const elemLocalToGlobal = elementSubRegion.localToGlobalMap().toViewConst(); finiteElement::FiniteElementBase const & fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); @@ -286,19 +290,20 @@ void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, { using FE_TYPE = TYPEOFREF( finiteElement ); - localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); PreComputeSourcesAndReceivers:: Compute3DSourceAndReceiverConstantsWithDAS < EXEC_POLICY, FE_TYPE > ( elementSubRegion.size(), - numFacesPerElem, - X, + facesToNodes, + nodeCoords, + nodeLocalToGlobal, + elemLocalToGlobal, + nodesToElements, + baseElemsToNodes, elemGhostRank, elemsToNodes, elemsToFaces, elemCenter, - faceNormal, - faceCenter, sourceCoordinates, sourceIsAccessible, sourceNodeIds, @@ -363,11 +368,12 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() applyFreeSurfaceBC( 0.0, domain ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - precomputeSourceAndReceiverTerm( mesh, regionNames ); + MeshLevel & baseMesh = domain.getMeshBodies().getGroup< MeshBody >( meshBodyName ).getBaseDiscretization(); + precomputeSourceAndReceiverTerm( baseMesh, mesh, regionNames ); NodeManager & nodeManager = mesh.getNodeManager(); FaceManager & faceManager = mesh.getFaceManager(); diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp index bdb15dcff18..2c9ec3c37ee 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.hpp @@ -132,6 +132,7 @@ class ElasticWaveEquationSEM : public WaveSolverBase * @param time_n time at the beginning of the step * @param dt the perscribed timestep * @param cycleNumber the current cycle number + * @param meshBodyName the name of the mesh body * @param domain the domain object * @return return the timestep that was achieved during the step. */ @@ -173,10 +174,11 @@ class ElasticWaveEquationSEM : public WaveSolverBase /** * @brief Locate sources and receivers position in the mesh elements, evaluate the basis functions at each point and save them to the * corresponding elements nodes. + * @param baseMesh the level-0 mesh * @param mesh mesh of the computational domain * @param regionNames the names of the region you loop on */ - virtual void precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + virtual void precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** * @brief Apply free surface condition to the face define in the geometry box from the xml diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp index c688e457c6d..c6fbf016edd 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp @@ -33,14 +33,16 @@ struct PreComputeSourcesAndReceivers * @tparam EXEC_POLICY execution policy * @tparam FE_TYPE finite element type * @param[in] size the number of cells in the subRegion - * @param[in] numFacesPerElem number of faces per element - * @param[in] nodeCoords coordinates of the nodes + * @param[in] baseFacesToNodes face to node map + * @param[in] baseNodeCoords coordinates of the nodes + * @param[in] baseNodeLocalToGlobal local to global index map for nodes + * @param[in] elementLocalToGlobal local to global index map for elements + * @param[in] baseNodesToElements node to element map for the base mesh + * @param[in] baseElemsToNodes element to node map for the base mesh * @param[in] elemGhostRank rank of the ghost element * @param[in] elemsToNodes map from element to nodes * @param[in] elemsToFaces map from element to faces * @param[in] elemCenter coordinates of the element centers - * @param[in] faceNormal normal of each faces - * @param[in] faceCenter coordinates of the center of a face * @param[in] sourceCoordinates coordinates of the source terms * @param[out] sourceIsAccessible flag indicating whether the source is accessible or not * @param[out] sourceNodeIds indices of the nodes of the element where the source is located @@ -58,14 +60,16 @@ struct PreComputeSourcesAndReceivers template< typename EXEC_POLICY, typename FE_TYPE > static void Compute1DSourceAndReceiverConstants( localIndex const size, - localIndex const numFacesPerElem, - arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + ArrayOfArraysView< localIndex const > const baseFacesToNodes, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const baseNodeCoords, + arrayView1d< globalIndex const > const baseNodeLocalToGlobal, + arrayView1d< globalIndex const > const elementLocalToGlobal, + ArrayOfArraysView< localIndex const > const baseNodesToElements, + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes, arrayView1d< integer const > const elemGhostRank, arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes, arrayView2d< localIndex const > const elemsToFaces, arrayView2d< real64 const > const & elemCenter, - arrayView2d< real64 const > const faceNormal, - arrayView2d< real64 const > const faceCenter, arrayView2d< real64 const > const sourceCoordinates, arrayView1d< localIndex > const sourceIsAccessible, arrayView2d< localIndex > const sourceNodeIds, @@ -98,22 +102,24 @@ struct PreComputeSourcesAndReceivers real64 const coords[3] = { sourceCoordinates[isrc][0], sourceCoordinates[isrc][1], sourceCoordinates[isrc][2] }; - bool const sourceFound = - WaveSolverUtils::locateSourceElement( numFacesPerElem, - center, - faceNormal, - faceCenter, - elemsToFaces[k], - coords ); + computationalGeometry::isPointInsideConvexPolyhedronRobust( k, + baseNodeCoords, + elemsToFaces, + baseFacesToNodes, + baseNodesToElements, + baseNodeLocalToGlobal, + elementLocalToGlobal, + center, + coords ); if( sourceFound ) { real64 coordsOnRefElem[3]{}; WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, - elemsToNodes[k], - nodeCoords, + baseElemsToNodes[k], + baseNodeCoords, coordsOnRefElem ); sourceIsAccessible[isrc] = 1; @@ -148,18 +154,21 @@ struct PreComputeSourcesAndReceivers real64 coordsOnRefElem[3]{}; bool const receiverFound = - WaveSolverUtils::locateSourceElement( numFacesPerElem, - center, - faceNormal, - faceCenter, - elemsToFaces[k], - coords ); + computationalGeometry::isPointInsideConvexPolyhedronRobust( k, + baseNodeCoords, + elemsToFaces, + baseFacesToNodes, + baseNodesToElements, + baseNodeLocalToGlobal, + elementLocalToGlobal, + center, + coords ); if( receiverFound && elemGhostRank[k] < 0 ) { WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, - elemsToNodes[k], - nodeCoords, + baseElemsToNodes[k], + baseNodeCoords, coordsOnRefElem ); receiverIsLocal[ircv] = 1; @@ -187,14 +196,16 @@ struct PreComputeSourcesAndReceivers * @tparam EXEC_POLICY execution policy * @tparam FE_TYPE finite element type * @param[in] size the number of cells in the subRegion - * @param[in] numFacesPerElem number of faces per element - * @param[in] nodeCoords coordinates of the nodes + * @param[in] baseFacesToNodes face to node map of the base mesh + * @param[in] baseNodeCoords coordinates of the nodes of the base mesh + * @param[in] baseNodeLocalToGlobal local to global index map for nodes of the base mesh + * @param[in] elementLocalToGlobal local to global index map for elements (for the base or high order mesh) + * @param[in] baseNodesToElements local node to element map for the base mesh + * @param[in] baseElemsToNodes element to node map for the base mesh * @param[in] elemGhostRank rank of the ghost element * @param[in] elemsToNodes map from element to nodes * @param[in] elemsToFaces map from element to faces * @param[in] elemCenter coordinates of the element centers - * @param[in] faceNormal normal of each faces - * @param[in] faceCenter coordinates of the center of a face * @param[in] sourceCoordinates coordinates of the source terms * @param[out] sourceIsAccessible flag indicating whether the source is accessible or not * @param[out] sourceElem element where a source is located @@ -215,14 +226,16 @@ struct PreComputeSourcesAndReceivers static void Compute1DSourceAndReceiverConstantsWithElementsAndRegionStorage( localIndex const size, localIndex const regionIndex, - localIndex const numFacesPerElem, - arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + ArrayOfArraysView< localIndex const > const baseFacesToNodes, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const baseNodeCoords, + arrayView1d< globalIndex const > const baseNodeLocalToGlobal, + arrayView1d< globalIndex const > const elementLocalToGlobal, + ArrayOfArraysView< localIndex const > const baseNodesToElements, + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes, arrayView1d< integer const > const elemGhostRank, arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes, arrayView2d< localIndex const > const elemsToFaces, arrayView2d< real64 const > const & elemCenter, - arrayView2d< real64 const > const faceNormal, - arrayView2d< real64 const > const faceCenter, arrayView2d< real64 const > const sourceCoordinates, arrayView1d< localIndex > const sourceIsAccessible, arrayView1d< localIndex > const sourceElem, @@ -261,20 +274,22 @@ struct PreComputeSourcesAndReceivers sourceCoordinates[isrc][2] }; bool const sourceFound = - WaveSolverUtils::locateSourceElement( numFacesPerElem, - center, - faceNormal, - faceCenter, - elemsToFaces[k], - coords ); - + computationalGeometry::isPointInsideConvexPolyhedronRobust( k, + baseNodeCoords, + elemsToFaces, + baseFacesToNodes, + baseNodesToElements, + baseNodeLocalToGlobal, + elementLocalToGlobal, + center, + coords ); if( sourceFound ) { real64 coordsOnRefElem[3]{}; WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, - elemsToNodes[k], - nodeCoords, + baseElemsToNodes[k], + baseNodeCoords, coordsOnRefElem ); sourceIsAccessible[isrc] = 1; @@ -311,18 +326,21 @@ struct PreComputeSourcesAndReceivers real64 coordsOnRefElem[3]{}; bool const receiverFound = - WaveSolverUtils::locateSourceElement( numFacesPerElem, - center, - faceNormal, - faceCenter, - elemsToFaces[k], - coords ); + computationalGeometry::isPointInsideConvexPolyhedronRobust( k, + baseNodeCoords, + elemsToFaces, + baseFacesToNodes, + baseNodesToElements, + baseNodeLocalToGlobal, + elementLocalToGlobal, + center, + coords ); if( receiverFound && elemGhostRank[k] < 0 ) { WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, - elemsToNodes[k], - nodeCoords, + baseElemsToNodes[k], + baseNodeCoords, coordsOnRefElem ); receiverIsLocal[ircv] = 1; receiverElem[ircv] = k; @@ -352,14 +370,16 @@ struct PreComputeSourcesAndReceivers * @tparam EXEC_POLICY execution policy * @tparam FE_TYPE finite element type * @param[in] size the number of cells in the subRegion - * @param[in] numFacesPerElem number of face on an element - * @param[in] nodeCoords coordinates of the nodes + * @param[in] baseFacesToNodes face to node map + * @param[in] baseNodeCoords coordinates of the nodes + * @param[in] baseNodeLocalToGlobal local to global index map for nodes + * @param[in] elementLocalToGlobal local to global index map for elements + * @param[in] baseNodesToElements node to element map for the base mesh + * @param[in] baseElemsToNodes element to node map for the base mesh * @param[in] elemGhostRank array containing the ghost rank * @param[in] elemsToNodes map from element to nodes * @param[in] elemsToFaces map from element to faces * @param[in] elemCenter coordinates of the element centers - * @param[in] faceNormal array containing the normal of all faces - * @param[in] faceCenter array containing the center of all faces * @param[in] sourceCoordinates coordinates of the source terms * @param[out] sourceIsAccessible flag indicating whether the source is accessible or not * @param[out] sourceNodeIds indices of the nodes of the element where the source is located @@ -385,14 +405,16 @@ struct PreComputeSourcesAndReceivers template< typename EXEC_POLICY, typename FE_TYPE > static void Compute3DSourceAndReceiverConstantsWithDAS( localIndex const size, - localIndex const numFacesPerElem, - arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + ArrayOfArraysView< localIndex const > const baseFacesToNodes, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const baseNodeCoords, + arrayView1d< globalIndex const > const baseNodeLocalToGlobal, + arrayView1d< globalIndex const > const elementLocalToGlobal, + ArrayOfArraysView< localIndex const > const baseNodesToElements, + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & baseElemsToNodes, arrayView1d< integer const > const elemGhostRank, arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes, arrayView2d< localIndex const > const elemsToFaces, arrayView2d< real64 const > const & elemCenter, - arrayView2d< real64 const > const faceNormal, - arrayView2d< real64 const > const faceCenter, arrayView2d< real64 const > const sourceCoordinates, arrayView1d< localIndex > const sourceIsAccessible, arrayView2d< localIndex > const sourceNodeIds, @@ -444,18 +466,21 @@ struct PreComputeSourcesAndReceivers { for( localIndex i=0; i<3; ++i ) { - xLocal[a][i] = nodeCoords( elemsToNodes( k, a ), i ); + xLocal[a][i] = baseNodeCoords( elemsToNodes( k, a ), i ); } } bool const sourceFound = - WaveSolverUtils::locateSourceElement( numFacesPerElem, - center, - faceNormal, - faceCenter, - elemsToFaces[k], - coords ); + computationalGeometry::isPointInsideConvexPolyhedronRobust( k, + baseNodeCoords, + elemsToFaces, + baseFacesToNodes, + baseNodesToElements, + baseNodeLocalToGlobal, + elementLocalToGlobal, + center, + coords ); if( sourceFound ) { @@ -463,8 +488,8 @@ struct PreComputeSourcesAndReceivers WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, - elemsToNodes[k], - nodeCoords, + baseElemsToNodes[k], + baseNodeCoords, coordsOnRefElem ); sourceIsAccessible[isrc] = 1; @@ -556,12 +581,15 @@ struct PreComputeSourcesAndReceivers receiverCenter[ 1 ] + receiverVector[ 1 ] * receiverLength * samplePointLocations[ iSample ], receiverCenter[ 2 ] + receiverVector[ 2 ] * receiverLength * samplePointLocations[ iSample ] }; bool const sampleFound = - WaveSolverUtils::locateSourceElement( numFacesPerElem, - center, - faceNormal, - faceCenter, - elemsToFaces[k], - coords ); + computationalGeometry::isPointInsideConvexPolyhedronRobust( k, + baseNodeCoords, + elemsToFaces, + baseFacesToNodes, + baseNodesToElements, + baseNodeLocalToGlobal, + elementLocalToGlobal, + center, + coords ); if( sampleFound && elemGhostRank[k] < 0 ) { real64 coordsOnRefElem[3]{}; @@ -571,13 +599,13 @@ struct PreComputeSourcesAndReceivers { for( localIndex i=0; i<3; ++i ) { - xLocal[a][i] = nodeCoords( elemsToNodes( k, a ), i ); + xLocal[a][i] = baseNodeCoords( elemsToNodes( k, a ), i ); } } WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, - elemsToNodes[k], - nodeCoords, + baseElemsToNodes[k], + baseNodeCoords, coordsOnRefElem ); real64 N[numNodesPerElem]; real64 gradN[numNodesPerElem][3]; @@ -603,12 +631,15 @@ struct PreComputeSourcesAndReceivers // determine if the current rank is the owner of this receiver real64 const coords[3] = { receiverCenter[ 0 ], receiverCenter[ 1 ], receiverCenter[ 2 ] }; bool const receiverFound = - WaveSolverUtils::locateSourceElement( numFacesPerElem, - center, - faceNormal, - faceCenter, - elemsToFaces[k], - coords ); + computationalGeometry::isPointInsideConvexPolyhedronRobust( k, + baseNodeCoords, + elemsToFaces, + baseFacesToNodes, + baseNodesToElements, + baseNodeLocalToGlobal, + elementLocalToGlobal, + center, + coords ); if( receiverFound && elemGhostRank[k] < 0 ) { receiverIsLocal[ ircv ] = 1; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp index 6cb5a459a0d..fe17f6d6a9d 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp @@ -207,9 +207,10 @@ class WaveSolverBase : public SolverBase /** * @brief Locate sources and receivers positions in the mesh elements, evaluate the basis functions at each point and save them to the * corresponding elements nodes. + * @param baseMesh the level-0 mesh * @param mesh mesh of the computational domain */ - virtual void precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) = 0; + virtual void precomputeSourceAndReceiverTerm( MeshLevel & baseMesh, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) = 0; /** * @brief Perform forward explicit step diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverUtils.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverUtils.hpp index b4087dcb2fb..962337e6b40 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverUtils.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverUtils.hpp @@ -364,8 +364,8 @@ struct WaveSolverUtils * @brief Convert a mesh element point coordinate into a coordinate on the reference element * @tparam FE_TYPE finite element type * @param[in] coords coordinate of the point - * @param[in] elemsToNodes map to obtaint global nodes from element index - * @param[in] nodeCoords array of mesh nodes coordinates + * @param[in] elemsToNodes element to node map for the base mesh + * @param[in] nodeCoords array of base mesh nodes coordinates * @param[out] coordsOnRefElem to contain the coordinate computed in the reference element */ template< typename FE_TYPE > @@ -373,14 +373,14 @@ struct WaveSolverUtils static void computeCoordinatesOnReferenceElement( real64 const (&coords)[3], arraySlice1d< localIndex const, cells::NODE_MAP_USD - 1 > const elemsToNodes, - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodeCoords, real64 (& coordsOnRefElem)[3] ) { // only the eight corners of the mesh cell are needed to compute the Jacobian real64 xLocal[8][3]{}; for( localIndex a = 0; a < 8; ++a ) { - LvArray::tensorOps::copy< 3 >( xLocal[a], nodeCoords[ elemsToNodes[ FE_TYPE::meshIndexToLinearIndex3D( a )] ] ); + LvArray::tensorOps::copy< 3 >( xLocal[a], nodeCoords[ elemsToNodes[ a ] ] ); } // coordsOnRefElem = invJ*(coords-coordsNode_0) real64 invJ[3][3]{};