H1_Wedge_Lagrange1_Gauss6.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_H1WEDGELAGRANGE1GAUSS6_HPP_
#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1WEDGELAGRANGE1GAUSS6_HPP_
 
#include "FiniteElementBase.hpp"
#include "LagrangeBasis1.hpp"
 
 
namespace geos
{
namespace finiteElement
{
 
class H1_Wedge_Lagrange1_Gauss6 final : public FiniteElementBase
{
public:
 
  using BASIS = LagrangeBasis1;
 
  constexpr static localIndex numNodes = 6;
 
  constexpr static localIndex numFaces = 5;
 
  constexpr static localIndex maxSupportPoints = numNodes;
 
  constexpr static localIndex numQuadraturePoints = 6;
 
  constexpr static int numSamplingPoints = numSamplingPointsPerDirection * numSamplingPointsPerDirection * numSamplingPointsPerDirection;
 
  GEOS_HOST_DEVICE
  virtual ~H1_Wedge_Lagrange1_Gauss6() 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 * 2 * step;
    if( (r+s) <= 1 )
    {
      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 - r;
      samplingPointCoord[1] = 1 - s;
      samplingPointCoord[2] = t;
    }
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcN( real64 const (&coords)[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 xi = pointCoord[2];
 
    N[0] = 4.0 * (1.0 - r - s) * s * LagrangeBasis1::valueBubble( xi );
    N[1] = 4.0 * (1.0 - r - s) * r * LagrangeBasis1::valueBubble( xi );
    N[2] = (1.0 - r - s) * r * s * LagrangeBasis1::value( 0, xi );
    N[3] = (1.0 - r - s) * r * s * LagrangeBasis1::value( 1, xi );
    N[4] = 4.0 * r * s * LagrangeBasis1::valueBubble( xi );
 
  }
 
  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] )
  {
    jacobianTransformation( q, X, J );
    return LvArray::tensorOps::invert< 3 >( J );
  }
 
 
private:
  constexpr static real64 parentVolume = 1.0;
 
  constexpr static real64 weight = parentVolume / numQuadraturePoints;
 
  constexpr static real64 quadratureCrossSectionCoord = 1.0 / 6.0;
 
  constexpr static real64 quadratureLongitudinalCoord = 1.0 / 1.732050807568877293528;
 
  template< typename T >
  GEOS_HOST_DEVICE
  inline
  constexpr static T linearMap( T const indexT, T const indexL )
  {
    return 2 * indexT + indexL;
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 parentCoords0( localIndex const a )
  {
    return 0.5 * ( a & 2 );
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 parentCoords1( localIndex const a )
  {
    return 0.25 * ( a & 4 );
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 parentCoords2( localIndex const a )
  {
    return -1.0 + 2 * ( a & 1 );
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords0( localIndex const q )
  {
    return quadratureCrossSectionCoord + 0.5  * parentCoords0( q );
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords1( localIndex const q )
  {
    return quadratureCrossSectionCoord + 0.5  * parentCoords1( q );
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords2( localIndex const q )
  {
    return parentCoords2( q ) * quadratureLongitudinalCoord;
  }
 
  GEOS_HOST_DEVICE
  static void jacobianTransformation( int const q,
                                      real64 const (&X)[numNodes][3],
                                      real64 ( &J )[3][3] );
 
  GEOS_HOST_DEVICE
  static void
    applyJacobianTransformationToShapeFunctionsDerivatives( int const q,
                                                            real64 const ( &invJ )[3][3],
                                                            real64 ( &gradN )[numNodes][3] );
 
};
 
 
GEOS_HOST_DEVICE
inline
void
H1_Wedge_Lagrange1_Gauss6::
  jacobianTransformation( int const q,
                          real64 const (&X)[numNodes][3],
                          real64 ( & J )[3][3] )
{
  real64 const r  = quadratureParentCoords0( q );
  real64 const s  = quadratureParentCoords1( q );
  real64 const xi = quadratureParentCoords2( q );
 
  real64 const psiTRI[3] = { 1.0 - r - s, r, s };
  real64 const psiLIN[2] = { 0.5 - 0.5*xi, 0.5 + 0.5*xi };
  constexpr real64 dpsiTRI[2][3] = { { -1.0, 1.0, 0.0 }, { -1.0, 0.0, 1.0 } };
  constexpr real64 dpsiLIN[2] = { -0.5, 0.5 };
 
  for( localIndex a=0; a<3; ++a )
  {
    for( localIndex b=0; b<2; ++b )
    {
      real64 const dNdXi[3] = { dpsiTRI[0][a] * psiLIN[b],
                                dpsiTRI[1][a] * psiLIN[b],
                                psiTRI[a] * dpsiLIN[b] };
      localIndex const nodeIndex = linearMap( a, b );
      for( int i = 0; i < 3; ++i )
      {
        for( int j = 0; j < 3; ++j )
        {
          J[i][j] = J[i][j] + dNdXi[ j ] * X[nodeIndex][i];
        }
      }
    }
  }
}
 
//*************************************************************************************************
 
GEOS_HOST_DEVICE
inline
void
H1_Wedge_Lagrange1_Gauss6::
  applyJacobianTransformationToShapeFunctionsDerivatives( int const q,
                                                          real64 const ( &invJ )[3][3],
                                                          real64 (& gradN)[numNodes][3] )
{
  real64 const r  = quadratureParentCoords0( q );
  real64 const s  = quadratureParentCoords1( q );
  real64 const xi = quadratureParentCoords2( q );
 
  real64 const psiTRI[3] = { 1.0 - r - s, r, s };
  real64 const psiLIN[2] = { 0.5 - 0.5*xi, 0.5 + 0.5*xi };
  constexpr real64 dpsiTRI[2][3] = { { -1.0, 1.0, 0.0 }, { -1.0, 0.0, 1.0 } };
  constexpr real64 dpsiLIN[2] = { -0.5, 0.5 };
 
  for( localIndex a=0; a<3; ++a )
  {
    for( localIndex b=0; b<2; ++b )
    {
      real64 const dNdXi[3] = { dpsiTRI[0][a] * psiLIN[b],
                                dpsiTRI[1][a] * psiLIN[b],
                                psiTRI[a] * dpsiLIN[b] };
      localIndex const nodeIndex = linearMap( a, b );
      for( int i = 0; i < 3; ++i )
      {
        gradN[nodeIndex][i] = 0.0;
        for( int j = 0; j < 3; ++j )
        {
          gradN[nodeIndex][i] = gradN[nodeIndex][i] + dNdXi[ j ] * invJ[j][i];
        }
      }
    }
  }
}
 
//*************************************************************************************************
 
GEOS_HOST_DEVICE
inline
void
H1_Wedge_Lagrange1_Gauss6::calcN( real64 const (&coords)[3],
                                  real64 (& N)[numNodes] )
{
  real64 const r  = coords[0];
  real64 const s  = coords[1];
  real64 const xi = coords[2];
 
  N[0] = 0.5*( 1.0 - r - s ) * ( 1.0 - xi );
  N[1] = 0.5*( 1.0 - r - s ) * ( 1.0 + xi );
  N[2] = 0.5* r * ( 1.0 - xi );
  N[3] = 0.5* r * ( 1.0 + xi );
  N[4] = 0.5* s * ( 1.0 - xi );
  N[5] = 0.5* s * ( 1.0 + xi );
 
}
 
 
GEOS_HOST_DEVICE
inline
void
H1_Wedge_Lagrange1_Gauss6::
  calcN( localIndex const q,
         real64 (& N)[numNodes] )
{
  real64 const pointCoord[3] = {quadratureParentCoords0( q ),
                                quadratureParentCoords1( q ),
                                quadratureParentCoords2( q )};
 
  calcN( pointCoord, N );
}
 
GEOS_HOST_DEVICE
inline
void H1_Wedge_Lagrange1_Gauss6::
  calcN( localIndex const q,
         StackVariables const & GEOS_UNUSED_PARAM( stack ),
         real64 ( & N )[numNodes] )
{
  return calcN( q, N );
}
 
//*************************************************************************************************
 
GEOS_HOST_DEVICE
inline
real64
H1_Wedge_Lagrange1_Gauss6::
  calcGradN( localIndex const q,
             real64 const (&X)[numNodes][3],
             real64 (& gradN)[numNodes][3] )
{
  real64 J[3][3] = {{0}};
 
  jacobianTransformation( q, X, J );
 
  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
 
  applyJacobianTransformationToShapeFunctionsDerivatives( q, J, gradN );
 
  return detJ * weight;
}
 
GEOS_HOST_DEVICE
inline
real64 H1_Wedge_Lagrange1_Gauss6::
  calcGradN( localIndex const q,
             real64 const (&X)[numNodes][3],
             StackVariables const & GEOS_UNUSED_PARAM( stack ),
             real64 ( & gradN )[numNodes][3] )
{
  return calcGradN( q, X, gradN );
}
 
GEOS_HOST_DEVICE
inline
real64
H1_Wedge_Lagrange1_Gauss6::calcGradFaceBubbleN( localIndex const q,
                                                real64 const (&X)[numNodes][3],
                                                real64 (& gradN)[numFaces][3] )
{
 
  real64 J[3][3] = {{0}};
 
  jacobianTransformation( q, X, J );
 
  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 xi = quadratureParentCoords2( q );
 
  dNdXi[0][0] = -4.0 * s * LagrangeBasis1::valueBubble( xi );                  // dN0/dr
  dNdXi[0][1] = 4.0 * (1.0 - r - 2.0 * s) * LagrangeBasis1::valueBubble( xi ); // dN0/ds
  dNdXi[0][2] = 4.0 *(1 - r - s) * s * LagrangeBasis1::gradientBubble( xi );   // dN0/dxi
 
  dNdXi[1][0] =  4.0 * (1.0 - 2.0 * r - s) * LagrangeBasis1::valueBubble( xi ); // dN1/dr
  dNdXi[1][1] =  -4.0 * r * LagrangeBasis1::valueBubble( xi );                  // dN1/ds
  dNdXi[1][2] =  4.0 * (1 - r - s) * r * LagrangeBasis1::gradientBubble( xi );  // dN1/dxi
 
  dNdXi[2][0] = (1.0 - 2.0 * r - s) * s * LagrangeBasis1::value( 0, xi );  // dN2/dr
  dNdXi[2][1] = (1.0 - r - 2.0 * s) * r * LagrangeBasis1::value( 0, xi );  // dN2/ds
  dNdXi[2][2] = (1 - r - s) * r * s * LagrangeBasis1::gradient( 0, xi );   // dN2/dxi
 
  dNdXi[3][0] = (1.0 - 2.0 * r - s) * s * LagrangeBasis1::value( 1, xi );  // dN3/dr
  dNdXi[3][1] = (1.0 - r - 2.0 * s) * r * LagrangeBasis1::value( 1, xi );  // dN3/ds
  dNdXi[3][2] = (1 - r - s) * r * s * LagrangeBasis1::gradient( 1, xi );   // dN3/dxi
 
  dNdXi[4][0] = 4.0 * s * LagrangeBasis1::valueBubble( xi );                    // dN4/dr
  dNdXi[4][1] = 4.0 * r * LagrangeBasis1::valueBubble( xi );                    // dN4/ds
  dNdXi[4][2] = 4.0 * r * s * LagrangeBasis1::gradientBubble( xi );             // dN4/dxi
 
  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;
}
 
//*************************************************************************************************
 
GEOS_HOST_DEVICE
inline
real64
H1_Wedge_Lagrange1_Gauss6::
  transformedQuadratureWeight( localIndex const q,
                               real64 const (&X)[numNodes][3] )
{
  real64 J[3][3] = {{0}};
 
  jacobianTransformation( q, X, J );
 
  return LvArray::tensorOps::determinant< 3 >( J ) * weight;
}
 
 
}
}
 
#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1WEDGELAGRANGE1GAUSS6_HPP_