ComputationalGeometry.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2018-2020 Total, S.A
 * Copyright (c) 2019-     GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */

#ifndef GEOSX_MESHUTILITIES_COMPUTATIONALGEOMETRY_HPP_
#define GEOSX_MESHUTILITIES_COMPUTATIONALGEOMETRY_HPP_

#include "common/DataTypes.hpp"
#include "common/DataLayouts.hpp"
#include "LvArray/src/output.hpp"
#include "LvArray/src/tensorOps.hpp"

namespace geosx
{
namespace computationalGeometry
{

constexpr real64 machinePrecision = std::numeric_limits< 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 >
void orderPointsCCW( arrayView2d< real64 > const & points,
                     NORMAL_TYPE const & normal )
{
  localIndex const numPoints = points.size( 0 );

  array2d< real64 > orderedPoints( numPoints, 3 );

  std::vector< int > indices( numPoints );
  std::vector< 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] );
  }
}

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< typename CENTER_TYPE, typename NORMAL_TYPE >
GEOSX_HOST_DEVICE
GEOSX_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 );

  GEOSX_ERROR_IF_LT( pointsIndices.size(), 2 );
  for( localIndex a=0; a<(pointsIndices.size()-2); ++a )
  {
    real64 v1[ 3 ], v2[ 3 ], vc[ 3 ];

    LvArray::tensorOps::copy< 3 >( v1, points[ pointsIndices[ a + 1 ] ] );
    LvArray::tensorOps::copy< 3 >( v2, points[ pointsIndices[ a + 2 ] ] );

    LvArray::tensorOps::copy< 3 >( vc, points[ pointsIndices[ 0 ] ] );
    LvArray::tensorOps::add< 3 >( vc, v1 );
    LvArray::tensorOps::add< 3 >( vc, v2 );

    LvArray::tensorOps::subtract< 3 >( v1, points[ pointsIndices[ 0 ] ] );
    LvArray::tensorOps::subtract< 3 >( v2, points[ pointsIndices[ 0 ] ] );

    real64 triangleNormal[ 3 ];
    LvArray::tensorOps::crossProduct( triangleNormal, v1, v2 );
    real64 const triangleArea = LvArray::tensorOps::l2Norm< 3 >( triangleNormal );

    LvArray::tensorOps::add< 3 >( normal, triangleNormal );

    area += triangleArea;
    LvArray::tensorOps::scaledAdd< 3 >( center, vc, triangleArea );
  }
  if( area > areaTolerance )
  {
    LvArray::tensorOps::scale< 3 >( center, 1.0 / ( area * 3.0 ) );
    LvArray::tensorOps::normalize< 3 >( normal );
    area *= 0.5;
  }
  else if( area < -areaTolerance )
  {
    for( localIndex a=0; a<pointsIndices.size(); ++a )
    {
      GEOSX_LOG_RANK( "Points: " << points[ pointsIndices[ a ] ] << " " << pointsIndices[ a ] );
    }
    GEOSX_ERROR( "Negative area found : " << area );
  }
  else
  {
    return 0.0;
  }

  return area;
}

template< typename NORMAL_TYPE >
GEOSX_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 >
GEOSX_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 ];

  GEOSX_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 >
GEOSX_HOST_DEVICE
GEOSX_FORCE_INLINE
int sign( T const val )
{
  return (T( 0 ) < val) - (val < T( 0 ));
}

template< typename POINT_TYPE >
GEOSX_HOST_DEVICE
bool IsPointInsidePolyhedron( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodeCoordinates,
                              array1d< array1d< localIndex > > const & faceNodeIndicies,
                              POINT_TYPE const & point,
                              real64 const areaTolerance = 0.0 )
{
  localIndex const numFaces = faceNodeIndicies.size( 0 );
  R1Tensor faceCenter, faceNormal;
  int prev_sign = 0;

  for( localIndex kf = 0; kf < numFaces; ++kf )
  {
    Centroid_3DPolygon( faceNodeIndicies[kf], nodeCoordinates, faceCenter, faceNormal, areaTolerance );

    LvArray::tensorOps::subtract< 3 >( faceCenter, point );
    int const s = sign( LvArray::tensorOps::AiBi< 3 >( faceNormal, faceCenter ) );

    // all dot products should be non-negative (for outward normals) or non-positive (for inward normals)
    if( prev_sign * s < 0 )
    {
      return false;
    }
    prev_sign = s;
  }

  return true;
}

template< typename VEC_TYPE >
GEOSX_HOST_DEVICE
void GetBoundingBox( localIndex const elemIndex,
                     arrayView2d< localIndex const, cells::NODE_MAP_USD > 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 >::min );

  // loop over all the vertices of the element to get the min and max coords
  for( localIndex a = 0; a < pointIndices.size( 1 ); ++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 );
}

GEOSX_HOST_DEVICE
inline
real64 HexVolume( real64 const X[][3] )
{
  real64 X7_X1[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[7] );
  LvArray::tensorOps::subtract< 3 >( X7_X1, X[1] );

  real64 X6_X0[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[6] );
  LvArray::tensorOps::subtract< 3 >( X6_X0, X[0] );

  real64 X7_X2[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[7] );
  LvArray::tensorOps::subtract< 3 >( X7_X2, X[2] );

  real64 X3_X0[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[3] );
  LvArray::tensorOps::subtract< 3 >( X3_X0, X[0] );

  real64 X5_X0[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[5] );
  LvArray::tensorOps::subtract< 3 >( X5_X0, X[0] );

  real64 X7_X4[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[7] );
  LvArray::tensorOps::subtract< 3 >( X7_X4, X[4] );

  real64 X7_X1plusX6_X0[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X7_X1 );
  LvArray::tensorOps::add< 3 >( X7_X1plusX6_X0, X6_X0 );

  real64 X7_X2plusX5_X0[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X7_X2 );
  LvArray::tensorOps::add< 3 >( X7_X2plusX5_X0, X5_X0 );

  real64 X7_X4plusX3_X0[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X7_X4 );
  LvArray::tensorOps::add< 3 >( X7_X4plusX3_X0, X3_X0 );

  real64 X7_X2crossX3_X0[3];
  LvArray::tensorOps::crossProduct( X7_X2crossX3_X0, X7_X2, X3_X0 );

  real64 X7_X2plusX5_X0crossX7_X4[3];
  LvArray::tensorOps::crossProduct( X7_X2plusX5_X0crossX7_X4, X7_X2plusX5_X0, X7_X4 );

  real64 X5_X0crossX7_X4plusX3_X0[3];
  LvArray::tensorOps::crossProduct( X5_X0crossX7_X4plusX3_X0, X5_X0, X7_X4plusX3_X0 );

  return 1.0/12.0 * ( LvArray::tensorOps::AiBi< 3 >( X7_X1plusX6_X0, X7_X2crossX3_X0 ) +
                      LvArray::tensorOps::AiBi< 3 >( X6_X0, X7_X2plusX5_X0crossX7_X4 ) +
                      LvArray::tensorOps::AiBi< 3 >( X7_X1, X5_X0crossX7_X4plusX3_X0 ) );
}

GEOSX_HOST_DEVICE
inline
real64 TetVolume( real64 const X[][3] )
{
  real64 X1_X0[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[1] );
  LvArray::tensorOps::subtract< 3 >( X1_X0, X[0] );

  real64 X2_X0[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[2] );
  LvArray::tensorOps::subtract< 3 >( X2_X0, X[0] );

  real64 X3_X0[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3( X[3] );
  LvArray::tensorOps::subtract< 3 >( X3_X0, X[0] );

  real64 X2_X0crossX3_X0[ 3 ];
  LvArray::tensorOps::crossProduct( X2_X0crossX3_X0, X2_X0, X3_X0 );

  return std::fabs( LvArray::tensorOps::AiBi< 3 >( X1_X0, X2_X0crossX3_X0 ) / 6.0 );
}

GEOSX_HOST_DEVICE
inline
real64 WedgeVolume( real64 const X[][3] )
{
  real64 const tet1[4][3] = { LVARRAY_TENSOROPS_INIT_LOCAL_3( X[0] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[1] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[2] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[4] ) };

  real64 const tet2[4][3] = { LVARRAY_TENSOROPS_INIT_LOCAL_3( X[0] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[2] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[4] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[5] ) };

  real64 const tet3[4][3] = { LVARRAY_TENSOROPS_INIT_LOCAL_3( X[0] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[3] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[4] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[5] ) };

  return TetVolume( tet1 ) + TetVolume( tet2 ) + TetVolume( tet3 );
}

GEOSX_HOST_DEVICE
inline
real64 PyramidVolume( real64 const X[][3] )
{
  real64 const tet1[4][3] = { LVARRAY_TENSOROPS_INIT_LOCAL_3( X[0] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[1] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[2] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[4] ) };

  real64 const tet2[4][3] = { LVARRAY_TENSOROPS_INIT_LOCAL_3( X[0] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[2] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[3] ),
                              LVARRAY_TENSOROPS_INIT_LOCAL_3( X[4] ) };

  return TetVolume( tet1 ) + TetVolume( tet2 );
}

} // namespace computationalGeometry
} // namespace geosx

#endif /* GEOSX_MESHUTILITIES_COMPUTATIONALGEOMETRY_HPP_ */