H1_Tetrahedron_Lagrange1_Gauss.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_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS_HPP_
#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS_HPP_
 
#include "FiniteElementBase.hpp"
#include "LagrangeBasis1.hpp"
 
 
namespace geos
{
namespace finiteElement
{
 
template< typename NUM_Q_POINTS >
class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase
{
public:
 
  static_assert( ( NUM_Q_POINTS::value == 1 ||
                   NUM_Q_POINTS::value == 5 ||
                   NUM_Q_POINTS::value == 14 ),
                 "NUM_Q_POINTS::value must be 1, 5, or 14!" );
 
  using BASIS = LagrangeBasis1;
 
  constexpr static localIndex numNodes = 4;
 
  constexpr static localIndex numFaces = 4;
 
  constexpr static localIndex maxSupportPoints = numNodes;
 
  constexpr static localIndex numQuadraturePoints = NUM_Q_POINTS::value;
 
  constexpr static int numSamplingPoints = numSamplingPointsPerDirection * numSamplingPointsPerDirection * numSamplingPointsPerDirection;
 
  GEOS_HOST_DEVICE
  virtual ~H1_Tetrahedron_Lagrange1_Gauss() override
  {}
 
  GEOS_HOST_DEVICE
  virtual localIndex getNumQuadraturePoints() const override
  {
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  static localIndex
  getNumQuadraturePoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  virtual localIndex getNumSupportPoints() const override
  {
    return numNodes;
  }
 
  GEOS_HOST_DEVICE
  virtual localIndex getMaxSupportPoints() const override
  {
    return maxSupportPoints;
  }
 
  GEOS_HOST_DEVICE
  static localIndex getNumSupportPoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numNodes;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static void getSamplingPointCoordInParentSpace( int const & linearIndex,
                                                  real64 (& samplingPointCoord)[3] )
  {
    int const i0 = linearIndex % numSamplingPointsPerDirection;
    int const i1 = ( (linearIndex - i0)/numSamplingPointsPerDirection ) % numSamplingPointsPerDirection;
    int const i2 = ( (linearIndex - i0)/numSamplingPointsPerDirection - i1 ) / numSamplingPointsPerDirection;
 
    real64 const step = 1. / ( numSamplingPointsPerDirection - 1 );
 
    real64 const r = i0 * step;
    real64 const s = i1 * step;
    real64 const t = i2 * step;
    if( (r+s) <= t )
    {
      samplingPointCoord[0] = r;
      samplingPointCoord[1] = s;
      samplingPointCoord[2] = t;
    }
    else
    {
      // if outside of the triangle need to reproject it. Points will be doubled though.
      samplingPointCoord[0] = 1 - t - r;
      samplingPointCoord[1] = 1 - t - s;
      samplingPointCoord[2] = t;
    }
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static void calcN( real64 const (&pointCoord)[3],
                     real64 ( &N )[numNodes] );
 
  GEOS_HOST_DEVICE
  static void calcN( localIndex const q,
                     real64 ( &N )[numNodes] );
 
  GEOS_HOST_DEVICE
  inline
  static void calcN( localIndex const q,
                     StackVariables const & stack,
                     real64 ( &N )[numNodes] );
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static void calcFaceBubbleN( real64 const (&pointCoord)[3],
                               real64 (& N)[numFaces] )
  {
 
    real64 const r = pointCoord[0];
    real64 const s = pointCoord[1];
    real64 const t = pointCoord[2];
 
    N[0] = (1 - r - s - t) * r * t;
    N[1] = (1 - r - s - t) * r * s;
    N[2] = (1 - r - s - t) * t * s;
    N[3] = r * s * t;
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcFaceBubbleN( localIndex const q,
                               real64 (& N)[numFaces] )
  {
 
    real64 const pointCoord[3] = {quadratureParentCoords0( q ),
                                  quadratureParentCoords1( q ),
                                  quadratureParentCoords2( q )};
 
    calcFaceBubbleN( pointCoord, N );
  }
 
  GEOS_HOST_DEVICE
  static real64 calcGradN( localIndex const q,
                           real64 const (&X)[numNodes][3],
                           real64 ( &gradN )[numNodes][3] );
 
  GEOS_HOST_DEVICE
  inline
  static real64 calcGradN( localIndex const q,
                           real64 const (&X)[numNodes][3],
                           StackVariables const & stack,
                           real64 ( &gradN )[numNodes][3] );
 
  GEOS_HOST_DEVICE
  static real64 calcGradFaceBubbleN( localIndex const q,
                                     real64 const (&X)[numNodes][3],
                                     real64 ( &gradN )[numFaces][3] );
 
  GEOS_HOST_DEVICE
  static real64 transformedQuadratureWeight( localIndex const q,
                                             real64 const (&X)[numNodes][3] );
 
  GEOS_HOST_DEVICE
  static real64 transformedQuadratureWeight( localIndex const q,
                                             real64 const (&X)[numNodes][3],
                                             StackVariables const & stack )
  { GEOS_UNUSED_VAR( stack ); return transformedQuadratureWeight( q, X ); }
 
  GEOS_HOST_DEVICE
  static real64 invJacobianTransformation( int const q,
                                           real64 const (&X)[numNodes][3],
                                           real64 ( & J )[3][3] )
  {
    GEOS_UNUSED_VAR( q, X );
//    jacobianTransformation( q, X, J );
    return LvArray::tensorOps::invert< 3 >( J );
  }
private:
  constexpr static real64 parentVolume = 1.0 / 6.0;
 
  constexpr static real64 weight = parentVolume;
  //constexpr static real64 weight = parentVolume / numQuadraturePoints;
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureWeight( localIndex const q )
  {
    if constexpr (numQuadraturePoints == 1)
    {
      real64 const w[numQuadraturePoints] = { 1.0 };
      return w[q];
    }
    else if constexpr (numQuadraturePoints == 5)
    {
      real64 const w[numQuadraturePoints] = {-4.0/5.0, 9.0/20.0, 9.0/20.0, 9.0/20.0, 9.0/20.0 };
      return w[q];
    }
    else if constexpr (numQuadraturePoints == 14)
    {
      real64 const w[numQuadraturePoints] = { 0.073493043116361949544,
                                              0.073493043116361949544,
                                              0.073493043116361949544,
                                              0.073493043116361949544,
                                              0.11268792571801585080,
                                              0.11268792571801585080,
                                              0.11268792571801585080,
                                              0.11268792571801585080,
                                              0.042546020777081466438,
                                              0.042546020777081466438,
                                              0.042546020777081466438,
                                              0.042546020777081466438,
                                              0.042546020777081466438,
                                              0.042546020777081466438 };
      return w[q];
    }
 
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords0( localIndex const q )
  {
 
    if constexpr (numQuadraturePoints == 1)
    {
      real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 5)
    {
      real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/2.0, 1.0/6.0, 1.0/6.0, 1.0/6.0 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 14)
    {
      real64 const qCoords[numQuadraturePoints] = { 0.72179424906732632079,
                                                    0.092735250310891226402,
                                                    0.092735250310891226402,
                                                    0.092735250310891226402,
                                                    0.067342242210098170608,
                                                    0.31088591926330060980,
                                                    0.31088591926330060980,
                                                    0.31088591926330060980,
                                                    0.045503704125649649492,
                                                    0.045503704125649649492,
                                                    0.045503704125649649492,
                                                    0.45449629587435035051,
                                                    0.45449629587435035051,
                                                    0.45449629587435035051 };
      return qCoords[q];
    }
 
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords1( localIndex const q )
  {
 
    if constexpr (numQuadraturePoints == 1)
    {
      real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 5)
    {
      real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/6.0, 1.0/2.0, 1.0/6.0, 1.0/6.0 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 14)
    {
      real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402,
                                                    0.72179424906732632079,
                                                    0.092735250310891226402,
                                                    0.092735250310891226402,
                                                    0.31088591926330060980,
                                                    0.067342242210098170608,
                                                    0.31088591926330060980,
                                                    0.31088591926330060980,
                                                    0.045503704125649649492,
                                                    0.45449629587435035051,
                                                    0.45449629587435035051,
                                                    0.045503704125649649492,
                                                    0.045503704125649649492,
                                                    0.45449629587435035051 };
      return qCoords[q];
    }
 
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords2( localIndex const q )
  {
 
    if constexpr (numQuadraturePoints == 1)
    {
      real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 5)
    {
      real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/6.0, 1.0/6.0, 1.0/2.0, 1.0/6.0 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 14)
    {
      real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402,
                                                    0.092735250310891226402,
                                                    0.72179424906732632079,
                                                    0.092735250310891226402,
                                                    0.31088591926330060980,
                                                    0.31088591926330060980,
                                                    0.067342242210098170608,
                                                    0.31088591926330060980,
                                                    0.45449629587435035051,
                                                    0.045503704125649649492,
                                                    0.45449629587435035051,
                                                    0.045503704125649649492,
                                                    0.45449629587435035051,
                                                    0.045503704125649649492 };
 
      return qCoords[q];
    }
 
  }
 
  GEOS_HOST_DEVICE
  static real64 determinantJacobianTransformation( real64 const (&X)[numNodes][3] );
 
};
 
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
real64
H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::
determinantJacobianTransformation( real64 const (&X)[numNodes][3] )
{
  return ( X[1][0] - X[0][0] )*( ( X[2][1] - X[0][1] )*( X[3][2] - X[0][2] ) - ( X[3][1] - X[0][1] )*( X[2][2] - X[0][2] ) )
         + ( X[1][1] - X[0][1] )*( ( X[3][0] - X[0][0] )*( X[2][2] - X[0][2] ) - ( X[2][0] - X[0][0] )*( X[3][2] - X[0][2] ) )
         + ( X[1][2] - X[0][2] )*( ( X[2][0] - X[0][0] )*( X[3][1] - X[0][1] ) - ( X[3][0] - X[0][0] )*( X[2][1] - X[0][1] ) );
}
 
//*************************************************************************************************
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
void
H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::
calcN( real64 const (&coords)[3],
       real64 (& N)[numNodes] )
{
  real64 const r = coords[0];
  real64 const s = coords[1];
  real64 const t = coords[2];
 
  // single quadrature point (centroid), i.e.  r = s = t = 1/4
  N[0] = 1 - r - s - t;
  N[1] = r;
  N[2] = s;
  N[3] = t;
}
 
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
void
H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::
calcN( localIndex const q,
       real64 (& N)[numNodes] )
{
  real64 const pointCoord[3] = {quadratureParentCoords0( q ),
                                quadratureParentCoords1( q ),
                                quadratureParentCoords2( q )};
 
  calcN( pointCoord, N );
}
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
void H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::
calcN( localIndex const q,
       StackVariables const & GEOS_UNUSED_PARAM( stack ),
       real64 ( & N )[numNodes] )
{
  return calcN( q, N );
}
 
//*************************************************************************************************
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
real64
H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::
calcGradN( localIndex const q,
           real64 const (&X)[numNodes][3],
           real64 (& gradN)[numNodes][3] )
{
 
  gradN[0][0] =  X[1][1]*( X[3][2] - X[2][2] ) - X[2][1]*( X[3][2] - X[1][2] ) + X[3][1]*( X[2][2] - X[1][2] );
  gradN[0][1] = -X[1][0]*( X[3][2] - X[2][2] ) + X[2][0]*( X[3][2] - X[1][2] ) - X[3][0]*( X[2][2] - X[1][2] );
  gradN[0][2] =  X[1][0]*( X[3][1] - X[2][1] ) - X[2][0]*( X[3][1] - X[1][1] ) + X[3][0]*( X[2][1] - X[1][1] );
 
  gradN[1][0] = -X[0][1]*( X[3][2] - X[2][2] ) + X[2][1]*( X[3][2] - X[0][2] ) - X[3][1]*( X[2][2] - X[0][2] );
  gradN[1][1] =  X[0][0]*( X[3][2] - X[2][2] ) - X[2][0]*( X[3][2] - X[0][2] ) + X[3][0]*( X[2][2] - X[0][2] );
  gradN[1][2] = -X[0][0]*( X[3][1] - X[2][1] ) + X[2][0]*( X[3][1] - X[0][1] ) - X[3][0]*( X[2][1] - X[0][1] );
 
  gradN[2][0] =  X[0][1]*( X[3][2] - X[1][2] ) - X[1][1]*( X[3][2] - X[0][2] ) + X[3][1]*( X[1][2] - X[0][2] );
  gradN[2][1] = -X[0][0]*( X[3][2] - X[1][2] ) + X[1][0]*( X[3][2] - X[0][2] ) - X[3][0]*( X[1][2] - X[0][2] );
  gradN[2][2] =  X[0][0]*( X[3][1] - X[1][1] ) - X[1][0]*( X[3][1] - X[0][1] ) + X[3][0]*( X[1][1] - X[0][1] );
 
  gradN[3][0] = -X[0][1]*( X[2][2] - X[1][2] ) + X[1][1]*( X[2][2] - X[0][2] ) - X[2][1]*( X[1][2] - X[0][2] );
  gradN[3][1] =  X[0][0]*( X[2][2] - X[1][2] ) - X[1][0]*( X[2][2] - X[0][2] ) + X[2][0]*( X[1][2] - X[0][2] );
  gradN[3][2] = -X[0][0]*( X[2][1] - X[1][1] ) + X[1][0]*( X[2][1] - X[0][1] ) - X[2][0]*( X[1][1] - X[0][1] );
 
  real64 detJ = determinantJacobianTransformation( X );
  real64 factor = 1.0 / ( detJ );
 
  for( int i = 0; i < numNodes; ++i )
  {
    for( int j = 0; j < 3; ++j )
    {
      gradN[i][j] *= factor;
    }
  }
 
  return detJ * weight * quadratureWeight( q );
}
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
real64 H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::
calcGradN( localIndex const q,
           real64 const (&X)[numNodes][3],
           StackVariables const & GEOS_UNUSED_PARAM( stack ),
           real64 ( & gradN )[numNodes][3] )
{
  return calcGradN( q, X, gradN );
}
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
real64
H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::calcGradFaceBubbleN( localIndex const q,
                                                                     real64 const (&X)[numNodes][3],
                                                                     real64 (& gradN)[numFaces][3] )
{
 
  //real64 detJ = determinantJacobianTransformation( X );
 
  real64 J[3][3] = {{0}};
 
  J[0][0] = X[1][0] - X[0][0];
  J[0][1] = X[2][0] - X[0][0];
  J[0][2] = X[3][0] - X[0][0];
 
  J[1][0] = X[1][1] - X[0][1];
  J[1][1] = X[2][1] - X[0][1];
  J[1][2] = X[3][1] - X[0][1];
 
  J[2][0] = X[1][2] - X[0][2];
  J[2][1] = X[2][2] - X[0][2];
  J[2][2] = X[3][2] - X[0][2];
 
  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
 
  real64 dNdXi[numFaces][3] = {{0}};
 
  real64 const r = quadratureParentCoords0( q );
  real64 const s = quadratureParentCoords1( q );
  real64 const t = quadratureParentCoords2( q );
 
  dNdXi[0][0] = ( 1 - 2 * r - s - t ) * t; // dN0/dr
  dNdXi[0][1] = -r * t;                    // dN0/ds
  dNdXi[0][2] = ( 1 - r - s - 2 * t ) * r; // dN0/dt
 
  dNdXi[1][0] = ( 1 - 2 * r - s - t ) * s; // dN1/dr
  dNdXi[1][1] = ( 1 - r - 2 * s - t ) * r; // dN1/ds
  dNdXi[1][2] = -r * s;                    // dN1/dt
 
  dNdXi[2][0] = -t * s;                    // dN2/dr
  dNdXi[2][1] = ( 1 - r - 2 * s - t ) * t; // dN2/ds
  dNdXi[2][2] = ( 1 - r - s - 2 * t ) * s; // dN2/dt
 
  dNdXi[3][0] = t * s;                     // dN3/dr
  dNdXi[3][1] = r * t;                     // dN3/ds
  dNdXi[3][2] = r * s;                     // dN3/dt
 
  for( int fi=0; fi<numFaces; ++fi )
  {
    gradN[fi][0] = dNdXi[fi][0] * J[0][0] + dNdXi[fi][1] * J[1][0] + dNdXi[fi][2] * J[2][0];
    gradN[fi][1] = dNdXi[fi][0] * J[0][1] + dNdXi[fi][1] * J[1][1] + dNdXi[fi][2] * J[2][1];
    gradN[fi][2] = dNdXi[fi][0] * J[0][2] + dNdXi[fi][1] * J[1][2] + dNdXi[fi][2] * J[2][2];
  }
 
  return detJ * weight * quadratureWeight( q );
}
 
//*************************************************************************************************
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
real64
H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::
transformedQuadratureWeight( localIndex const q,
                             real64 const (&X)[numNodes][3] )
{
 
  real64 detJ =  determinantJacobianTransformation( X );
 
  return detJ * weight * quadratureWeight( q );
}
 
 
using H1_Tetrahedron_Lagrange1_Gauss1 = H1_Tetrahedron_Lagrange1_Gauss< std::integral_constant< int, 1 > >;
using H1_Tetrahedron_Lagrange1_Gauss5 = H1_Tetrahedron_Lagrange1_Gauss< std::integral_constant< int, 5 > >;
using H1_Tetrahedron_Lagrange1_Gauss14 = H1_Tetrahedron_Lagrange1_Gauss< std::integral_constant< int, 14 > >;
 
}
}
 
#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS_HPP_