Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix: robust point location for degenerate configurations #3202

Merged
merged 42 commits into from
Aug 5, 2024
Merged
Show file tree
Hide file tree
Changes from 25 commits
Commits
Show all changes
42 commits
Select commit Hold shift + click to select a range
73f8348
work in progress
sframba Jun 12, 2024
637cf71
Merge branch 'develop' of https://github.com/GEOS-DEV/GEOS into develop
sframba Jun 12, 2024
24b87bc
implemented method based on winding number
sframba Jul 2, 2024
ded94c2
Merge branch 'develop' into feature/robustInsideOutside
sframba Jul 2, 2024
a0aea88
using element global index for element choice in degenerate cases
sframba Jul 3, 2024
897217f
added global id criterion for points coincident with vertices, as well
sframba Jul 3, 2024
be87fa1
wrong condition for vertices
sframba Jul 3, 2024
1f6d82e
Merge branch 'develop' into feature/robustInsideOutside
sframba Jul 3, 2024
0421543
uncrustify
sframba Jul 3, 2024
4589bdc
Merge branch 'feature/robustInsideOutside' of https://github.com/GEOS…
sframba Jul 3, 2024
de21789
doxygen
sframba Jul 3, 2024
4075951
typo
sframba Jul 3, 2024
c617d2f
removed unused variable
sframba Jul 3, 2024
de5c110
another unused variable
sframba Jul 4, 2024
4f25702
Merge branch 'develop' into feature/robustInsideOutside
sframba Jul 10, 2024
023811b
Fix some typos in comments
acitrain Jul 10, 2024
e2fd255
Merge branch 'feature/robustInsideOutside' of https://github.com/GEOS…
acitrain Jul 10, 2024
ba49f23
Merge branch 'develop' into feature/robustInsideOutside
sframba Jul 12, 2024
89af49e
Merge branch 'develop' into feature/robustInsideOutside
sframba Jul 17, 2024
c845d9b
Merge branch 'develop' into feature/robustInsideOutside
sframba Jul 22, 2024
3a613b9
Merge branch 'develop' into feature/robustInsideOutside
acitrain Jul 24, 2024
4915455
Merge branch 'develop' into feature/robustInsideOutside
acitrain Jul 25, 2024
bca995d
review comments + high-order case
sframba Jul 29, 2024
4dffd52
Merge branch 'develop' into feature/robustInsideOutside
sframba Jul 29, 2024
c0375b6
restored omega initial value (changed by mistake)
sframba Jul 29, 2024
8dccedb
Changed order of methods
sframba Jul 29, 2024
51787c8
using base discretization for precomputing src/rcv
sframba Jul 30, 2024
eb2011e
Revert "using base discretization for precomputing src/rcv"
sframba Jul 31, 2024
6a97776
fixed issue for higher orders
sframba Jul 31, 2024
2252f8d
Merge branch 'develop' into feature/robustInsideOutside
sframba Jul 31, 2024
aa82ad3
Revert "fixed issue for higher orders"
sframba Jul 31, 2024
8a2187a
Merge branch 'feature/robustInsideOutside' of https://github.com/GEOS…
sframba Jul 31, 2024
4648f8f
enabled for high orders. Also removed extra high-order nodes where no…
sframba Aug 1, 2024
f473ff6
restored default node coordinate behaviour outside of precompute src …
sframba Aug 3, 2024
c1b553b
uncrustify + spurious changes
sframba Aug 3, 2024
a5a7b46
leftover spurious change
sframba Aug 3, 2024
00cf346
Merge branch 'develop' into feature/robustInsideOutside
sframba Aug 3, 2024
e51b95e
unused vars
sframba Aug 3, 2024
baa23ce
Merge branch 'feature/robustInsideOutside' of https://github.com/GEOS…
sframba Aug 3, 2024
2821c93
unused vars
sframba Aug 3, 2024
a458470
rebaseline
sframba Aug 3, 2024
71f5930
Merge branch 'develop' into feature/robustInsideOutside
rrsettgast Aug 4, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
399 changes: 398 additions & 1 deletion src/coreComponents/mesh/utilities/ComputationalGeometry.hpp

Large diffs are not rendered by default.

Original file line number Diff line number Diff line change
Expand Up @@ -160,8 +160,9 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev

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 = nodeManager.localToGlobalMap().toViewConst();
ArrayOfArraysView< localIndex const > const nodesToElements = nodeManager.elementList().toViewConst();
ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();

arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst();
arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView();
Expand Down Expand Up @@ -206,28 +207,28 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev
arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.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() );
finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
using FE_TYPE = TYPEOFREF( finiteElement );

localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement();

PreComputeSourcesAndReceivers::
Compute1DSourceAndReceiverConstantsWithElementsAndRegionStorage
< EXEC_POLICY, FE_TYPE >
( elementSubRegion.size(),
regionIndex,
numFacesPerElem,
facesToNodes,
X,
nodeLocalToGlobal,
elemLocalToGlobal,
nodesToElements,
elemGhostRank,
elemsToNodes,
elemsToFaces,
elemCenter,
faceNormal,
faceCenter,
sourceCoordinates,
sourceIsAccessible,
sourceElem,
Expand Down Expand Up @@ -262,7 +263,7 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro
MeshLevel & mesh,
arrayView1d< string const > const & regionNames )
{
precomputeSourceAndReceiverTerm( mesh, regionNames );
precomputeSourceAndReceiverTerm( mesh.getBaseDiscretization(), regionNames );

NodeManager & nodeManager = mesh.getNodeManager();
FaceManager & faceManager = mesh.getFaceManager();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -120,9 +120,9 @@ void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & me

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 = nodeManager.localToGlobalMap().toViewConst();
ArrayOfArraysView< localIndex const > const nodesToElements = nodeManager.elementList().toViewConst();
ArrayOfArraysView< localIndex const > const facesToNodesMap = faceManager.nodeList().toViewConst();


arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst();
Expand Down Expand Up @@ -164,27 +164,27 @@ void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & me
arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.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() );
finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
using FE_TYPE = TYPEOFREF( finiteElement );

localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement();

acousticVTIWaveEquationSEMKernels::
PrecomputeSourceAndReceiverKernel::
launch< EXEC_POLICY, FE_TYPE >
( elementSubRegion.size(),
numFacesPerElem,
facesToNodesMap,
nodeCoords32,
nodeLocalToGlobalMap,
elemLocalToGlobalMap,
nodesToElements,
elemGhostRank,
elemsToNodes,
elemsToFaces,
elemCenter,
faceNormal,
faceCenter,
sourceCoordinates,
sourceIsAccessible,
sourceNodeIds,
Expand Down Expand Up @@ -237,7 +237,7 @@ void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
MeshLevel & mesh,
arrayView1d< string const > const & regionNames )
{
precomputeSourceAndReceiverTerm( mesh, regionNames );
precomputeSourceAndReceiverTerm( mesh.getBaseDiscretization(), regionNames );

NodeManager & nodeManager = mesh.getNodeManager();
FaceManager & faceManager = mesh.getFaceManager();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -40,12 +40,13 @@ 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] facesToNodes face to node map
* @param[in] nodeCoords coordinates of the nodes
* @param[in] nodeLocalToGlobal local to global map for nodes
* @param[in] elementLocalToGlobal local to global map for elements
* @param[in] elemGhostRank the ghost ranks
* @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
Expand All @@ -64,14 +65,15 @@ struct PrecomputeSourceAndReceiverKernel
template< typename EXEC_POLICY, typename FE_TYPE >
static void
launch( localIndex const size,
localIndex const numFacesPerElem,
ArrayOfArraysView< localIndex const > const facesToNodes,
arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
arrayView1d< globalIndex const > const nodeLocalToGlobal,
arrayView1d< globalIndex const > const elementLocalToGlobal,
ArrayOfArraysView< localIndex const > const nodesToElements,
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,
Expand Down Expand Up @@ -106,12 +108,15 @@ struct PrecomputeSourceAndReceiverKernel
sourceCoordinates[isrc][2] };

bool const sourceFound =
WaveSolverUtils::locateSourceElement( numFacesPerElem,
center,
faceNormal,
faceCenter,
elemsToFaces[k],
coords );
computationalGeometry::isPointInsideConvexPolyhedronRobust( k,
nodeCoords,
elemsToFaces,
facesToNodes,
nodesToElements,
nodeLocalToGlobal,
elementLocalToGlobal,
center,
coords );
if( sourceFound )
{
real64 coordsOnRefElem[3]{};
Expand Down Expand Up @@ -153,14 +158,17 @@ 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,
nodeCoords,
elemsToFaces,
facesToNodes,
nodesToElements,
nodeLocalToGlobal,
elementLocalToGlobal,
center,
coords );
if( receiverFound && elemGhostRank[k] < 0 )
{
WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -127,9 +127,10 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,

arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const
nodeCoords32 = nodeManager.getField< fields::referencePosition32 >().toViewConst();
arrayView1d< globalIndex const > const nodeLocalToGlobalMap = nodeManager.localToGlobalMap().toViewConst();
ArrayOfArraysView< localIndex const > const nodesToElements = nodeManager.elementList().toViewConst();

arrayView2d< real64 const > const faceNormal = faceManager.faceNormal();
arrayView2d< real64 const > const faceCenter = faceManager.faceCenter();
ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();

arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst();
arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView();
Expand Down Expand Up @@ -170,29 +171,29 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,
arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.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() );
finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
using FE_TYPE = TYPEOFREF( finiteElement );

localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement();

{
GEOS_MARK_SCOPE( acousticWaveEquationSEMKernels::PrecomputeSourceAndReceiverKernel );
PreComputeSourcesAndReceivers::
Compute1DSourceAndReceiverConstants
< EXEC_POLICY, FE_TYPE >
( elementSubRegion.size(),
numFacesPerElem,
facesToNodes,
nodeCoords32,
nodeLocalToGlobalMap,
elemLocalToGlobalMap,
nodesToElements,
elemGhostRank,
elemsToNodes,
elemsToFaces,
elemCenter,
faceNormal,
faceCenter,
sourceCoordinates,
sourceIsAccessible,
sourceNodeIds,
Expand Down Expand Up @@ -254,7 +255,7 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
MeshLevel & mesh,
arrayView1d< string const > const & regionNames )
{
precomputeSourceAndReceiverTerm( mesh, regionNames );
precomputeSourceAndReceiverTerm( mesh.getBaseDiscretization(), regionNames );

NodeManager & nodeManager = mesh.getNodeManager();
FaceManager & faceManager = mesh.getFaceManager();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -208,8 +208,9 @@ void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLeve
FaceManager const & faceManager = mesh.getFaceManager();

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 = nodeManager.localToGlobalMap().toViewConst();
ArrayOfArraysView< localIndex const > const nodesToElements = nodeManager.elementList().toViewConst();
ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();

arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst();
arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView();
Expand Down Expand Up @@ -257,28 +258,28 @@ void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLeve
arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.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() );
finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
using FE_TYPE = TYPEOFREF( finiteElement );

localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement();

PreComputeSourcesAndReceivers::
Compute1DSourceAndReceiverConstantsWithElementsAndRegionStorage
< EXEC_POLICY, FE_TYPE >
( elementSubRegion.size(),
regionIndex,
numFacesPerElem,
facesToNodes,
X,
nodeLocalToGlobal,
elemLocalToGlobal,
nodesToElements,
elemGhostRank,
elemsToNodes,
elemsToFaces,
elemCenter,
faceNormal,
faceCenter,
sourceCoordinates,
sourceIsAccessible,
sourceElem,
Expand Down Expand Up @@ -314,7 +315,7 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou
MeshLevel & mesh,
arrayView1d< string const > const & regionNames )
{
precomputeSourceAndReceiverTerm( mesh, regionNames );
precomputeSourceAndReceiverTerm( mesh.getBaseDiscretization(), regionNames );

NodeManager & nodeManager = mesh.getNodeManager();
FaceManager & faceManager = mesh.getFaceManager();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -230,8 +230,9 @@ void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,

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 = nodeManager.localToGlobalMap().toViewConst();
ArrayOfArraysView< localIndex const > const nodesToElements = nodeManager.elementList().toViewConst();
ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();

arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst();
arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView();
Expand Down Expand Up @@ -279,26 +280,27 @@ void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,
arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.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() );
finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
using FE_TYPE = TYPEOFREF( finiteElement );

localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement();
PreComputeSourcesAndReceivers::
Compute3DSourceAndReceiverConstantsWithDAS
< EXEC_POLICY, FE_TYPE >
( elementSubRegion.size(),
numFacesPerElem,
facesToNodes,
X,
nodeLocalToGlobal,
elemLocalToGlobal,
nodesToElements,
elemGhostRank,
elemsToNodes,
elemsToFaces,
elemCenter,
faceNormal,
faceCenter,
sourceCoordinates,
sourceIsAccessible,
sourceNodeIds,
Expand Down Expand Up @@ -367,7 +369,7 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
MeshLevel & mesh,
arrayView1d< string const > const & regionNames )
{
precomputeSourceAndReceiverTerm( mesh, regionNames );
precomputeSourceAndReceiverTerm( mesh.getBaseDiscretization(), regionNames );

NodeManager & nodeManager = mesh.getNodeManager();
FaceManager & faceManager = mesh.getFaceManager();
Expand Down
Loading
Loading