H1_Pyramid_Lagrange1_Gauss5.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_H1PYRAMIDLAGRANGE1GAUSS5_HPP_
#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1PYRAMIDLAGRANGE1GAUSS5_HPP_
 
#include "FiniteElementBase.hpp"
#include "LagrangeBasis1.hpp"
 
 
namespace geos
{
namespace finiteElement
{
 
class H1_Pyramid_Lagrange1_Gauss5_impl : public FiniteElementBase_impl< 5, 5, 5 >
{
public:
 
  using BASIS = LagrangeBasis1;
 
 
  struct StackVariables
  {};
 
  template< typename SUBREGION_TYPE >
  struct MeshData {};
 
  H1_Pyramid_Lagrange1_Gauss5_impl() = default;
  ~H1_Pyramid_Lagrange1_Gauss5_impl() = default;
  H1_Pyramid_Lagrange1_Gauss5_impl( H1_Pyramid_Lagrange1_Gauss5_impl const & ) = default;
  H1_Pyramid_Lagrange1_Gauss5_impl & operator=( H1_Pyramid_Lagrange1_Gauss5_impl const & ) = default;
  H1_Pyramid_Lagrange1_Gauss5_impl( H1_Pyramid_Lagrange1_Gauss5_impl && ) = default;
  H1_Pyramid_Lagrange1_Gauss5_impl & operator=( H1_Pyramid_Lagrange1_Gauss5_impl && ) = default;
 
  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] )
  {
    GEOS_UNUSED_VAR( pointCoord, N );
    GEOS_ERROR( "Unsupported bubble functions for pyramid elements" );
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcFaceBubbleN( localIndex const q,
                               real64 (& N)[numFaces] )
  {
    GEOS_UNUSED_VAR( q, N );
    GEOS_ERROR( "Unsupported bubble functions for pyramid elements" );
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static
  real64 calcJacobian( localIndex const q,
                       real64 const (&X)[numNodes][3],
                       real64 ( &J )[3][3] );
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static
  real64 calcJacobian( localIndex const q,
                       real64 const (&X)[numNodes][3],
                       StackVariables const &stack,
                       real64 ( &J )[3][3] );
 
  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 weight = 81.0 / 100.0;
 
  constexpr static real64 weightDelta  = 125.0 / 27.0 - weight;
 
  constexpr static real64 quadratureCrossSectionCoord = 0.584237394672177;
 
  constexpr static real64 quadratureLongitudinalCoordNeg = -2.0 / 3.0;
 
  constexpr static real64 quadratureLongitudinalCoordDelta = 16.0 / 15.0;
 
  template< typename T >
  GEOS_HOST_DEVICE
  inline
  constexpr static T linearMap( T const i, T const j )
  {
    return i + 2 * j;
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 parentCoords0( localIndex const a )
  {
    return -1.0 + 2.0 * ( a & 1 ) + 0.25 * ( a & 4 );
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 parentCoords1( localIndex const a )
  {
    return -1.0 + ( a & 2 ) + 0.25 * ( a & 4 );
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 parentCoords2( localIndex const a )
  {
    return -1.0 + 0.5 * ( a & 4 );
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords0( localIndex const q )
  {
    return parentCoords0( q ) * quadratureCrossSectionCoord;
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords1( localIndex const q )
  {
    return parentCoords1( q ) * quadratureCrossSectionCoord;
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords2( localIndex const q )
  {
    return quadratureLongitudinalCoordNeg + 0.5 * ( 1 + parentCoords2( q ) ) * quadratureLongitudinalCoordDelta;
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureWeight( localIndex const q )
  {
    return weight + 0.5 * ( 1 + parentCoords2( q ) ) * weightDelta;
  }
 
  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_Pyramid_Lagrange1_Gauss5_impl::
  jacobianTransformation( int const q,
                          real64 const (&X)[numNodes][3],
                          real64 ( & J )[3][3] )
{
  real64 const quadratureCoords[3] = { quadratureParentCoords0( q ),
                                       quadratureParentCoords1( q ),
                                       quadratureParentCoords2( q ) };
 
  real64 const psi0[2] = { 0.5 - 0.5*quadratureCoords[0],
                           0.5 + 0.5*quadratureCoords[0] };
  real64 const psi1[2] = { 0.5 - 0.5*quadratureCoords[1],
                           0.5 + 0.5*quadratureCoords[1] };
  real64 const psi2 = 0.5 - 0.5*quadratureCoords[2];
  constexpr real64 dpsi[2] = { -0.5, 0.5 };
 
  // Contributions from basis functions paired with base nodes
  for( localIndex a=0; a<2; ++a )
  {
    for( localIndex b=0; b<2; ++b )
    {
      real64 const dNdXi[3] = { dpsi[a] * psi1[b] * psi2,
                                psi0[a] * dpsi[b] * psi2,
                                psi0[a] * psi1[b] * dpsi[0] };
      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];
        }
      }
    }
  }
 
  // Contribution from the basis function paired with the apex nodes
  for( int i = 0; i < 3; ++i )
  {
    J[i][2] = J[i][2] + dpsi[1] * X[4][i];
  }
}
 
//*************************************************************************************************
 
GEOS_HOST_DEVICE
inline
void
H1_Pyramid_Lagrange1_Gauss5_impl::
  applyJacobianTransformationToShapeFunctionsDerivatives( int const q,
                                                          real64 const ( &invJ )[3][3],
                                                          real64 (& gradN)[numNodes][3] )
{
  real64 const quadratureCoords[3] = { quadratureParentCoords0( q ),
                                       quadratureParentCoords1( q ),
                                       quadratureParentCoords2( q ) };
 
  real64 const psi0[2] = { 0.5*( 1.0 - quadratureCoords[0] ),
                           0.5*( 1.0 + quadratureCoords[0] ) };
  real64 const psi1[2] = { 0.5*( 1.0 - quadratureCoords[1] ),
                           0.5*( 1.0 + quadratureCoords[1] ) };
  real64 const psi2 = 0.5*( 1.0 - quadratureCoords[2]);
  constexpr real64 dpsi[2] = { -0.5, 0.5 };
 
  // Contributions from basis functions paired with base nodes
  for( localIndex a=0; a<2; ++a )
  {
    for( localIndex b=0; b<2; ++b )
    {
      real64 const dNdXi[3] = { dpsi[a] * psi1[b] * psi2,
                                psi0[a] * dpsi[b] * psi2,
                                psi0[a] * psi1[b] * dpsi[0] };
      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];
        }
      }
    }
  }
 
  // Contribution from the basis function paired with the apex nodes
  for( int i = 0; i < 3; ++i )
  {
    gradN[4][i] = dpsi[1] * invJ[2][i];
  }
}
 
//*************************************************************************************************
 
GEOS_HOST_DEVICE
inline
void
H1_Pyramid_Lagrange1_Gauss5_impl::
  calcN( real64 const ( &pointCoord )[3],
         real64 ( & N )[numNodes] )
{
  N[0] = 0.125*( 1.0 - pointCoord[0] ) * ( 1.0 - pointCoord[1] ) * ( 1.0 - pointCoord[2] );
  N[1] = 0.125*( 1.0 + pointCoord[0] ) * ( 1.0 - pointCoord[1] ) * ( 1.0 - pointCoord[2] );
  N[2] = 0.125*( 1.0 - pointCoord[0] ) * ( 1.0 + pointCoord[1] ) * ( 1.0 - pointCoord[2] );
  N[3] = 0.125*( 1.0 + pointCoord[0] ) * ( 1.0 + pointCoord[1] ) * ( 1.0 - pointCoord[2] );
  N[4] = 0.5*( 1.0 + pointCoord[2] );
}
 
GEOS_HOST_DEVICE
inline
void
H1_Pyramid_Lagrange1_Gauss5_impl::
  calcN( localIndex const q,
         real64 ( & N )[numNodes] )
{
  real64 const xi[3] = { quadratureParentCoords0( q ),
                         quadratureParentCoords1( q ),
                         quadratureParentCoords2( q ) };
 
  calcN( xi, N );
}
 
GEOS_HOST_DEVICE
inline
void H1_Pyramid_Lagrange1_Gauss5_impl::
  calcN( localIndex const q,
         StackVariables const & GEOS_UNUSED_PARAM( stack ),
         real64 ( & N )[numNodes] )
{
  return calcN( q, N );
}
 
//*************************************************************************************************
 
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64
H1_Pyramid_Lagrange1_Gauss5_impl::calcJacobian( localIndex const q,
                                                real64 const (&X)[numNodes][3],
                                                real64 (& J)[3][3] )
{
  for( int i = 0; i < 3; ++i )
  {
    for( int j = 0; j < 3; ++j )
    {
      J[i][j] = 0;
    }
  }
  jacobianTransformation( q, X, J );
  real64 const detJ = LvArray::tensorOps::determinant< 3 >( J );
  return detJ * quadratureWeight( q );
}
 
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64
H1_Pyramid_Lagrange1_Gauss5_impl::calcJacobian( localIndex const q,
                                                real64 const (&X)[numNodes][3],
                                                StackVariables const & GEOS_UNUSED_PARAM( stack ),
                                                real64 (& J)[3][3] )
{
  return calcJacobian( q, X, J );
}
 
GEOS_HOST_DEVICE
inline
real64 H1_Pyramid_Lagrange1_Gauss5_impl::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 * quadratureWeight( q );
}
 
GEOS_HOST_DEVICE
inline
real64 H1_Pyramid_Lagrange1_Gauss5_impl::
  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_Pyramid_Lagrange1_Gauss5_impl::calcGradFaceBubbleN( localIndex const q,
                                                       real64 const (&X)[numNodes][3],
                                                       real64 (& gradN)[numFaces][3] )
{
  GEOS_UNUSED_VAR( q, X, gradN );
  GEOS_ERROR( "Unsupported bubble functions for pyramid elements" );
  return 0.0;
}
 
//*************************************************************************************************
 
GEOS_HOST_DEVICE
inline
real64
H1_Pyramid_Lagrange1_Gauss5_impl::
  transformedQuadratureWeight( localIndex const q,
                               real64 const (&X)[numNodes][3] )
{
  real64 J[3][3] = {{0}};
 
  jacobianTransformation( q, X, J );
 
  return LvArray::tensorOps::determinant< 3 >( J ) * quadratureWeight( q );
}
 
 
class H1_Pyramid_Lagrange1_Gauss5 final : public H1_Pyramid_Lagrange1_Gauss5_impl,
  public FiniteElementBase
{
public:
 
  using BASIS = LagrangeBasis1;
 
  using ImplType = H1_Pyramid_Lagrange1_Gauss5_impl;
 
 
  H1_Pyramid_Lagrange1_Gauss5():
    FiniteElementBase( numNodes,
                       maxSupportPoints,
                       numQuadraturePoints )
  {}
 
  virtual ~H1_Pyramid_Lagrange1_Gauss5() override final = default;
 
  H1_Pyramid_Lagrange1_Gauss5_impl * getImpl()
  {
    return static_cast< H1_Pyramid_Lagrange1_Gauss5_impl * >(this);
  }
 
  H1_Pyramid_Lagrange1_Gauss5_impl const * getImpl() const
  {
    return static_cast< H1_Pyramid_Lagrange1_Gauss5_impl const * >(this);
  }
 
};
 
}
}
#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1PYRAMIDLAGRANGE1GAUSS5_HPP_