ComputationalGeometry.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_MESH_UTILITIES_COMPUTATIONALGEOMETRY_HPP_
#define GEOS_MESH_UTILITIES_COMPUTATIONALGEOMETRY_HPP_
 
#include "common/DataTypes.hpp"
#include "common/DataLayouts.hpp"
#include "finiteElement/elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp"
#include "finiteElement/elementFormulations/H1_Pyramid_Lagrange1_Gauss5.hpp"
#include "finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp"
#include "finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp"
#include "LvArray/src/output.hpp"
#include "LvArray/src/tensorOps.hpp"
 
namespace geos
{
namespace computationalGeometry
{
 
constexpr real64 machinePrecision = LvArray::NumericLimits< real64 >::epsilon;
 
template< typename LINEDIR_TYPE,
          typename POINT_TYPE,
          typename NORMAL_TYPE,
          typename ORIGIN_TYPE,
          typename INTPOINT_TYPE >
void LinePlaneIntersection( LINEDIR_TYPE const & lineDir,
                            POINT_TYPE const & linePoint,
                            NORMAL_TYPE const & planeNormal,
                            ORIGIN_TYPE const & planeOrigin,
                            INTPOINT_TYPE & intersectionPoint )
{
  /* Find intersection line plane
   * line equation: p - (d*lineDir + linePoing) = 0;
   * plane equation: ( p - planeOrigin) * planeNormal = 0;
   * d = (planeOrigin - linePoint) * planeNormal / (lineDir * planeNormal )
   * pInt = d*lineDir+linePoint;
   */
  real64 dummy[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( planeOrigin );
  LvArray::tensorOps::subtract< 3 >( dummy, linePoint );
  real64 const d = LvArray::tensorOps::AiBi< 3 >( dummy, planeNormal ) /
                   LvArray::tensorOps::AiBi< 3 >( lineDir, planeNormal );
 
  LvArray::tensorOps::copy< 3 >( intersectionPoint, linePoint );
  LvArray::tensorOps::scaledAdd< 3 >( intersectionPoint, lineDir, d );
}
 
template< typename NORMAL_TYPE >
array1d< int >  orderPointsCCW( arrayView2d< real64 > const & points,
                                NORMAL_TYPE const & normal )
{
  localIndex const numPoints = points.size( 0 );
 
  array2d< real64 > orderedPoints( numPoints, 3 );
 
  array1d< int > indices( numPoints );
  array1d< real64 > angle( numPoints );
 
  // compute centroid of the set of points
  real64 centroid[3];
  LvArray::tensorOps::fill< 3 >( centroid, 0 );
  for( localIndex a = 0; a < numPoints; ++a )
  {
    LvArray::tensorOps::add< 3 >( centroid, points[ a ] );
    indices[ a ] = a;
  }
 
  LvArray::tensorOps::scale< 3 >( centroid, 1.0 / numPoints );
 
  real64 v0[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( centroid );
  LvArray::tensorOps::subtract< 3 >( v0, points[ 0 ] );
  LvArray::tensorOps::normalize< 3 >( v0 );
 
  // compute angles
  angle[ 0 ] = 0;
  for( localIndex a = 1; a < numPoints; ++a )
  {
    real64 v[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( centroid );
    LvArray::tensorOps::subtract< 3 >( v, points[ a ] );
    real64 const dot = LvArray::tensorOps::AiBi< 3 >( v, v0 );
 
    real64 crossProduct[ 3 ];
    LvArray::tensorOps::crossProduct( crossProduct, v, v0 );
    real64 const det = LvArray::tensorOps::AiBi< 3 >( normal, crossProduct );
 
    angle[ a ] = std::atan2( det, dot );
  }
 
  // sort the indices
  std::sort( indices.begin(), indices.end(), [&]( int i, int j ) { return angle[ i ] < angle[ j ]; } );
 
  // copy the points in the reorderedPoints array.
  for( localIndex a=0; a < numPoints; a++ )
  {
    // fill in with ordered
    LvArray::tensorOps::copy< 3 >( orderedPoints[ a ], points[ indices[ a ] ] );
  }
 
  for( localIndex a = 0; a < numPoints; a++ )
  {
    LvArray::tensorOps::copy< 3 >( points[a], orderedPoints[a] );
  }
 
  return indices;
}
 
template< typename NORMAL_TYPE >
real64 ComputeSurfaceArea( arrayView2d< real64 const > const & points,
                           NORMAL_TYPE const && normal )
{
  real64 surfaceArea = 0.0;
 
  array2d< real64 > orderedPoints( points.size( 0 ), 3 );
 
  for( localIndex a = 0; a < points.size( 0 ); a++ )
  {
    LvArray::tensorOps::copy< 3 >( orderedPoints[a], points[a] );
  }
 
  orderPointsCCW( orderedPoints, normal );
 
  for( localIndex a = 0; a < points.size( 0 ) - 2; ++a )
  {
    real64 v1[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( orderedPoints[ a + 1 ] );
    real64 v2[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( orderedPoints[ a + 2 ] );
 
    LvArray::tensorOps::subtract< 3 >( v1, orderedPoints[ 0 ] );
    LvArray::tensorOps::subtract< 3 >( v2, orderedPoints[ 0 ] );
 
    real64 triangleNormal[ 3 ];
    LvArray::tensorOps::crossProduct( triangleNormal, v1, v2 );
    surfaceArea += LvArray::tensorOps::l2Norm< 3 >( triangleNormal );
  }
 
  return surfaceArea * 0.5;
}
 
template< localIndex DIMENSION, typename POINT_COORDS_TYPE >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64 computeDiameter( POINT_COORDS_TYPE points,
                        localIndex const & numPoints )
{
  real64 diameter = 0;
  for( localIndex numPoint = 0; numPoint < numPoints; ++numPoint )
  {
    for( localIndex numOthPoint = 0; numOthPoint < numPoint; ++numOthPoint )
    {
      real64 candidateDiameter = 0.0;
      for( localIndex i = 0; i < DIMENSION; ++i )
      {
        real64 coordDiff = points[numPoint][i] - points[numOthPoint][i];
        candidateDiameter += coordDiff * coordDiff;
      }
      if( diameter < candidateDiameter )
      {
        diameter = candidateDiameter;
      }
    }
  }
  return LvArray::math::sqrt< real64 >( diameter );
}
 
template< typename CENTER_TYPE, typename NORMAL_TYPE >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64 centroid_3DPolygon( arraySlice1d< localIndex const > const pointsIndices,
                           arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & points,
                           CENTER_TYPE && center,
                           NORMAL_TYPE && normal,
                           real64 const areaTolerance = 0.0 )
{
  real64 area = 0.0;
  LvArray::tensorOps::fill< 3 >( center, 0 );
  LvArray::tensorOps::fill< 3 >( normal, 0 );
 
  localIndex const numberOfPoints = pointsIndices.size();
 
  GEOS_ERROR_IF_LT( numberOfPoints, 2 );
 
  real64 current[ 3 ], next[ 3 ], origin[ 3 ], crossProduct[ 3 ];
 
  LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ numberOfPoints - 1 ] ] );
  LvArray::tensorOps::copy< 3 >( origin, points[ pointsIndices[ 0 ]] );
 
  for( localIndex a=0; a<numberOfPoints; )
  {
    LvArray::tensorOps::copy< 3 >( current, points[ pointsIndices[ a++ ]] );
    LvArray::tensorOps::scaledAdd< 3 >( current, origin, -1. );
    LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ a % numberOfPoints ] ] );
    LvArray::tensorOps::scaledAdd< 3 >( next, origin, -1. );
 
    LvArray::tensorOps::crossProduct( crossProduct, current, next );
 
    LvArray::tensorOps::add< 3 >( normal, crossProduct );
    LvArray::tensorOps::add< 3 >( center, next );
  }
 
  area = LvArray::tensorOps::l2Norm< 3 >( normal );
  LvArray::tensorOps::scale< 3 >( center, 1.0 / numberOfPoints );
  LvArray::tensorOps::scaledAdd< 3 >( center, origin, 1. );
 
  if( area > areaTolerance )
  {
    LvArray::tensorOps::normalize< 3 >( normal );
    area *= 0.5;
  }
  else if( area < -areaTolerance )
  {
    for( localIndex a=0; a<numberOfPoints; ++a )
    {
      GEOS_LOG_RANK( "Points: " << points[ pointsIndices[ a ] ] << " " << pointsIndices[ a ] );
    }
#if defined(GEOS_DEVICE_COMPILE)
    GEOS_ERROR( "Negative area found" );
#else
    GEOS_ERROR( GEOS_FMT( "Negative area found : {}", area ) );
#endif
  }
  else
  {
    for( localIndex a=0; a<numberOfPoints; ++a )
    {
      GEOS_LOG_RANK( "Points: " << points[ pointsIndices[ a ] ] << " " << pointsIndices[ a ] );
    }
#if defined(GEOS_DEVICE_COMPILE)
    GEOS_ERROR( "Null area found" );
#else
    GEOS_ERROR( GEOS_FMT( "Null area found : {}", area ) );
#endif
 
    return 0.0;
  }
 
  return area;
}
 
template< typename NORMAL_TYPE >
GEOS_HOST_DEVICE
void FixNormalOrientation_3D( NORMAL_TYPE && normal )
{
  real64 const orientationTolerance = 10 * machinePrecision;
 
  // Orient local normal in global sense.
  // First check: align with z direction
  if( normal[ 2 ] <= -orientationTolerance )
  {
    LvArray::tensorOps::scale< 3 >( normal, -1.0 );
  }
  else if( std::fabs( normal[ 2 ] ) < orientationTolerance )
  {
    // If needed, second check: align with y direction
    if( normal[ 1 ] <= -orientationTolerance )
    {
      LvArray::tensorOps::scale< 3 >( normal, -1.0 );
    }
    else if( fabs( normal[ 1 ] ) < orientationTolerance )
    {
      // If needed, third check: align with x direction
      if( normal[ 0 ] <= -orientationTolerance )
      {
        LvArray::tensorOps::scale< 3 >( normal, -1.0 );
      }
    }
  }
}
 
template< typename NORMAL_TYPE, typename MATRIX_TYPE >
GEOS_HOST_DEVICE
void RotationMatrix_3D( NORMAL_TYPE const & normal,
                        MATRIX_TYPE && rotationMatrix )
{
  real64 m1[ 3 ] = { normal[ 2 ], 0.0, -normal[ 0 ] };
  real64 m2[ 3 ] = { 0.0, normal[ 2 ], -normal[ 1 ] };
  real64 const norm_m1 = LvArray::tensorOps::l2Norm< 3 >( m1 );
  real64 const norm_m2 = LvArray::tensorOps::l2Norm< 3 >( m2 );
 
  // If present, looks for a vector with 0 norm
  // Fix the uncertain case of norm_m1 very close to norm_m2
  if( norm_m1+1.e+2*machinePrecision > norm_m2 )
  {
    LvArray::tensorOps::crossProduct( m2, normal, m1 );
    LvArray::tensorOps::normalize< 3 >( m2 );
    LvArray::tensorOps::normalize< 3 >( m1 );
  }
  else
  {
    LvArray::tensorOps::crossProduct( m1, normal, m2 );
    LvArray::tensorOps::scale< 3 >( m1, -1 );
    LvArray::tensorOps::normalize< 3 >( m1 );
    LvArray::tensorOps::normalize< 3 >( m2 );
  }
 
  // Save everything in the standard form (3x3 rotation matrix)
  rotationMatrix[ 0 ][ 0 ] = normal[ 0 ];
  rotationMatrix[ 1 ][ 0 ] = normal[ 1 ];
  rotationMatrix[ 2 ][ 0 ] = normal[ 2 ];
  rotationMatrix[ 0 ][ 1 ] = m1[ 0 ];
  rotationMatrix[ 1 ][ 1 ] = m1[ 1 ];
  rotationMatrix[ 2 ][ 1 ] = m1[ 2 ];
  rotationMatrix[ 0 ][ 2 ] = m2[ 0 ];
  rotationMatrix[ 1 ][ 2 ] = m2[ 1 ];
  rotationMatrix[ 2 ][ 2 ] = m2[ 2 ];
 
  GEOS_ERROR_IF( fabs( LvArray::tensorOps::determinant< 3 >( rotationMatrix ) - 1.0 ) > 1.e+1 * machinePrecision,
                 "Rotation matrix with determinant different from +1.0" );
}
 
template< typename T >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
int sign( T const val )
{
  return (T( 0 ) < val) - (val < T( 0 ));
}
 
template< typename POINT_TYPE >
GEOS_HOST_DEVICE
bool isPointInsidePolyhedron( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodeCoordinates,
                              arraySlice1d< localIndex const > const & faceIndices,
                              ArrayOfArraysView< localIndex const > const & facesToNodes,
                              POINT_TYPE const & elemCenter,
                              POINT_TYPE const & point,
                              real64 const areaTolerance = 0.0 )
{
  localIndex const numFaces = faceIndices.size();
  R1Tensor faceCenter, faceNormal, cellToFaceVec;
 
  for( localIndex kf = 0; kf < numFaces; ++kf )
  {
    // compute the face normal at this face
    localIndex const faceIndex = faceIndices[kf];
    centroid_3DPolygon( facesToNodes[faceIndex], nodeCoordinates, faceCenter, faceNormal, areaTolerance );
 
    // make sure that the normal is outward pointing
    LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter );
    LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter );
    if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 )
    {
      LvArray::tensorOps::scale< 3 >( faceNormal, -1 );
    }
 
    // compute the vector face center to query point
    LvArray::tensorOps::subtract< 3 >( faceCenter, point );
    int const s = sign( LvArray::tensorOps::AiBi< 3 >( faceNormal, faceCenter ) );
 
    // all dot products should be non-negative (we enforce outward normals)
    if( s < 0 )
    {
      return false;
    }
  }
  return true;
}
 
template< typename POLYGON_TYPE, typename POINT_TYPE >
bool isPointInPolygon2d( POLYGON_TYPE const & polygon,
                         integer n,
                         POINT_TYPE const & point,
                         real64 const tol = 1e-10 )
{
  integer count = 0;
 
  for( integer i = 0; i < n; ++i )
  {
    auto const & p1 = polygon[i];
    auto const & p2 = polygon[(i + 1) % n];
 
    real64 y1 = p1[1], y2 = p2[1];
    real64 x1 = p1[0], x2 = p2[0];
    real64 py = point[1], px = point[0];
 
    // quick reject in y with tolerance
    if( py + tol < std::min( y1, y2 ) || py - tol > std::max( y1, y2 ) )
      continue;
 
    // check if point is (approximately) on the segment
    // parametric t for projection on segment in y (if segment vertical-ish use x)
    if( std::abs( (x2 - x1) * (py - y1) - (px - x1) * (y2 - y1) ) < tol *
        ( std::hypot( x2 - x1, y2 - y1 ) + 1.0 ) )
    {
      // ensure px is between x1,x2 and py between y1,y2 (with tol)
      if( px + tol >= std::min( x1, x2 ) && px - tol <= std::max( x1, x2 ) &&
          py + tol >= std::min( y1, y2 ) && py - tol <= std::max( y1, y2 ) )
        return true; // on boundary -> consider inside
    }
 
    // ignore nearly-horizontal edges for intersection counting
    if( std::abs( y2 - y1 ) < tol )
      continue;
 
    // compute x coordinate of intersection of horizontal line py with segment p1-p2
    real64 xIntersect = x1 + (py - y1) * (x2 - x1) / (y2 - y1);
 
    // count crossing where intersection is strictly to the right of point (robust with tol)
    if( px < xIntersect - tol )
      ++count;
  }
 
  return (count % 2) == 1;
}
 
template< typename POLYGON_TYPE, typename POINT_TYPE >
bool isPointInPolygon3d( POLYGON_TYPE const & polygon,
                         integer const n,
                         POINT_TYPE const & point,
                         real64 const tol = 1e-10 )
{
  // Check if the point lies in the plane of the polygon
  auto const & p0 = polygon[0];
  POINT_TYPE normal = {0, 0, 0};
  for( integer i = 1; i < n - 1; i++ )
  {
    auto const & p1 = polygon[i];
    auto const & p2 = polygon[i + 1];
    normal[0] += (p1[1] - p0[1]) * (p2[2] - p0[2]) - (p1[2] - p0[2]) * (p2[1] - p0[1]);
    normal[1] += (p1[2] - p0[2]) * (p2[0] - p0[0]) - (p1[0] - p0[0]) * (p2[2] - p0[2]);
    normal[2] += (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p1[1] - p0[1]) * (p2[0] - p0[0]);
  }
 
  real64 const dist = normal[0] * point[0] + normal[1] * point[1] + normal[2] * point[2] -(normal[0] * p0[0] + normal[1] * p0[1] + normal[2] * p0[2]);
 
  if( std::abs( dist ) > tol )
  {
    return false;
  }
 
  // Determine the dominant component of the normal vector
  int dominantIndex = 0;
  if( std::abs( normal[1] ) > std::abs( normal[0] ))
  {
    dominantIndex = 1;
  }
  if( std::abs( normal[2] ) > std::abs( normal[dominantIndex] ))
  {
    dominantIndex = 2;
  }
 
  // Project the polygon and the point onto a 2D plane
  POLYGON_TYPE projectedPolygon( n );
  POINT_TYPE projectedPoint;
  if( dominantIndex == 0 )  // X is dominant, project onto YZ plane
  {
    for( int i = 0; i < n; i++ )
    {
      projectedPolygon[i][0] = polygon[i][1];
      projectedPolygon[i][1] = polygon[i][2];
    }
    projectedPoint[0] = point[1];
    projectedPoint[1] = point[2];
  }
  else if( dominantIndex == 1 )  // Y is dominant, project onto XZ plane
  {
    for( int i = 0; i < n; i++ )
    {
      projectedPolygon[i][0] = polygon[i][0];
      projectedPolygon[i][1] = polygon[i][2];
    }
    projectedPoint[0] = point[0];
    projectedPoint[1] = point[2];
  }
  else  // Z is dominant, project onto XY plane
  {
    for( int i = 0; i < n; i++ )
    {
      projectedPolygon[i][0] = polygon[i][0];
      projectedPolygon[i][1] = polygon[i][1];
    }
    projectedPoint[0] = point[0];
    projectedPoint[1] = point[1];
  }
 
  return isPointInPolygon2d( projectedPolygon, n, projectedPoint );
}
 
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;
}
 
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 ) );
}
 
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;
}
template< typename ... LIST_TYPE >
GEOS_HOST_DEVICE
int findVertexRefElement( arraySlice1d< localIndex const > const & nodeElements,
                          arrayView1d< globalIndex const > const & elementGlobalIndex )
{
  localIndex minElement = -1;
  globalIndex minElementGID = LvArray::NumericLimits< 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;
}
 
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 = LvArray::NumericLimits< 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;
}
 
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 = LvArray::NumericLimits< 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;
}
 
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 = LvArray::NumericLimits< 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;
}
 
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;
}
 
template< typename NODE_MAP_TYPE, typename VEC_TYPE >
GEOS_HOST_DEVICE
void getBoundingBox( localIndex const elemIndex,
                     NODE_MAP_TYPE const & pointIndices,
                     arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & pointCoordinates,
                     VEC_TYPE && boxDims )
{
  // This holds the min coordinates of the set in each direction
  R1Tensor minCoords = { LvArray::NumericLimits< real64 >::max,
                         LvArray::NumericLimits< real64 >::max,
                         LvArray::NumericLimits< real64 >::max };
 
  // boxDims is used to hold the max coordinates.
  LvArray::tensorOps::fill< 3 >( boxDims, LvArray::NumericLimits< real64 >::lowest );
 
  // loop over all the vertices of the element to get the min and max coords
  for( localIndex a = 0; a < pointIndices[elemIndex].size(); ++a )
  {
    localIndex const id = pointIndices( elemIndex, a );
    for( localIndex d = 0; d < 3; ++d )
    {
      minCoords[ d ] = fmin( minCoords[ d ], pointCoordinates( id, d ) );
      boxDims[ d ] = fmax( boxDims[ d ], pointCoordinates( id, d ) );
    }
  }
 
  LvArray::tensorOps::subtract< 3 >( boxDims, minCoords );
}
 
template< typename FE_TYPE >
GEOS_HOST_DEVICE inline
real64 elementVolume( real64 const (&X)[FE_TYPE::numNodes][3] )
{
  real64 result{};
  for( localIndex q=0; q<FE_TYPE::numQuadraturePoints; ++q )
  {
    result = result + FE_TYPE::transformedQuadratureWeight( q, X );
  }
  return result;
}
 
GEOS_HOST_DEVICE
inline
real64 hexahedronVolume( real64 const (&X)[8][3] )
{
  return elementVolume< finiteElement::H1_Hexahedron_Lagrange1_GaussLegendre2 >( X );
}
 
GEOS_HOST_DEVICE
inline
real64 tetrahedronVolume( real64 const (&X)[4][3] )
{
  return elementVolume< finiteElement::H1_Tetrahedron_Lagrange1_Gauss1 >( X );
}
 
GEOS_HOST_DEVICE
inline
real64 wedgeVolume( real64 const (&X)[6][3] )
{
  return elementVolume< finiteElement::H1_Wedge_Lagrange1_Gauss6 >( X );
}
 
GEOS_HOST_DEVICE
inline
real64 pyramidVolume( real64 const (&X)[5][3] )
{
  return elementVolume< finiteElement::H1_Pyramid_Lagrange1_Gauss5 >( X );
}
 
template< integer N >
GEOS_HOST_DEVICE
inline
real64 prismVolume( real64 const (&X)[2*N][3] )
{
  static_assert( N > 4,
                 "Function prismVolume can be called for a prism with N-sided polygon base where N > 5." );
 
  real64 result{};
 
  // Compute the barycenters of the prism bases
  real64 XGBot[3]{};
  real64 XGTop[3]{};
  for( integer a = 0; a < N; ++a )
  {
    LvArray::tensorOps::add< 3 >( XGBot, X[a] );
  }
  for( integer a = N; a < 2 * N; ++a )
  {
    LvArray::tensorOps::add< 3 >( XGTop, X[a] );
  }
  LvArray::tensorOps::scale< 3 >( XGBot, 1.0 / N );
  LvArray::tensorOps::scale< 3 >( XGTop, 1.0 / N );
 
  real64 XWedge[6][3];
  for( int a = 0; a < N - 1; ++a )
  {
 
    LvArray::tensorOps::copy< 3 >( XWedge[0], X[a] );
    LvArray::tensorOps::copy< 3 >( XWedge[1], X[a+N] );
    LvArray::tensorOps::copy< 3 >( XWedge[2], X[a+1] );
    LvArray::tensorOps::copy< 3 >( XWedge[3], X[a+1+N] );
    LvArray::tensorOps::copy< 3 >( XWedge[4], XGBot );
    LvArray::tensorOps::copy< 3 >( XWedge[5], XGTop );
    result = result + computationalGeometry::elementVolume< finiteElement::H1_Wedge_Lagrange1_Gauss6 >( XWedge );
  }
  LvArray::tensorOps::copy< 3 >( XWedge[0], X[N-1] );
  LvArray::tensorOps::copy< 3 >( XWedge[1], X[2*N-1] );
  LvArray::tensorOps::copy< 3 >( XWedge[2], X[0] );
  LvArray::tensorOps::copy< 3 >( XWedge[3], X[N] );
  LvArray::tensorOps::copy< 3 >( XWedge[4], XGBot );
  LvArray::tensorOps::copy< 3 >( XWedge[5], XGTop );
  result = result + computationalGeometry::elementVolume< finiteElement::H1_Wedge_Lagrange1_Gauss6 >( XWedge );
  return result;
}
 
} /* namespace computationalGeometry */
} /* namespace geos */
 
#endif /* GEOS_MESH_UTILITIES_COMPUTATIONALGEOMETRY_HPP_ */