Pk_Pyramid_BCD.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 Total, S.A
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2018-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_PKPYRAMIDBCD_HPP_
#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_PKPYRAMIDBCD_HPP_
 
#include "FiniteElementBase.hpp"
#include "denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp"
#include <utility>
 
 
 
namespace geos
{
namespace finiteElement
{
 
 
template< int ORDER >
class Pk_Pyramid_BCD final : public FiniteElementBase
{
public:
 
  static constexpr int order = ORDER;
 
  constexpr static localIndex numNodes = ( ORDER + 1 ) * ( ORDER + 2 ) * ( 2 * ORDER + 3 ) / 6;
 
  constexpr static localIndex numFaces = 5;
 
  constexpr static localIndex maxSupportPoints = numNodes;
 
  constexpr static localIndex numQuadraturePoints = numNodes;
 
  constexpr static localIndex numModes = numNodes;
 
  USING_FINITEELEMENTBASE
  virtual ~Pk_Pyramid_BCD() = default;
 
  GEOS_HOST_DEVICE
  virtual localIndex getNumQuadraturePoints() const override
  {
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static localIndex getNumQuadraturePoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  virtual localIndex getNumSupportPoints() const override
  {
    return numNodes;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  virtual localIndex getMaxSupportPoints() const override
  {
    return maxSupportPoints;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static localIndex getNumSupportPoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numNodes;
  }
 
 
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static localIndex getNumModes( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numNodes;
  }
 
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static localIndex linearIndex3DVal( const localIndex qa, localIndex const qb, localIndex const qc )
  {
    localIndex index = 0;
 
    for( int l = 0; l < qc; ++l )
    {
      int n = order + 1 - l;
      index += n * n;
    }
    int n_k = order + 1 - qc;
    index += qa + qb * n_k;
 
    return index;
  }
 
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr void generatePointsCoordinates( real64 (& coords)[numNodes][3] )
  {
    // Generate the coordinates of the support points based on the order
    if constexpr (ORDER == 1)
    {
      coords[0][0] = -1.0; coords[0][1] = 1.0; coords[0][2] = 0.0;
      coords[1][0] =  -1.0; coords[1][1] = -1.0; coords[1][2] = 0.0;
      coords[2][0] =  1.0; coords[2][1] =  -1.0; coords[2][2] = 0.0;
      coords[3][0] = 1.0; coords[3][1] =  1.0; coords[3][2] = 0.0;
      coords[4][0] =  0.0; coords[4][1] =  0.0; coords[4][2] = 1.0;
    }
    else if constexpr (ORDER == 2)
    {
      //TODO
    }
 
  }
 
 
  template< typename FUNC >
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr void generateIndexes( FUNC && func )
  {
 
    for( localIndex i = 0; i <= ORDER; ++i )
    {
      for( localIndex j = 0; j <= ORDER; ++j )
      {
        localIndex maxIj = LvArray::math::max( i, j );
        for( localIndex k = 0; k <= ORDER - maxIj; ++k )
        {
          func( i, j, k );
        }
      }
 
    }
 
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static real64 jacobianDeterminant( real64 const (&X)[5][3] )
 
  {
    real64 m[3][3] = {};
    for( int i = 0; i < 4; i++ )
    {
      for( int j = 0; j < 3; j++ )
      {
        m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
      }
    }
    return LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) );
  }
 
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static real64 faceJacobianDeterminant( localIndex face,
                                         real64 const (&X)[5][3] )
  {
    int i1 = ( face + 1 ) % 4;
    int i2 = ( face + 2 ) % 4;
    int i3 = ( face + 3 ) % 4;
 
    real64 ab[3] = { X[ i2 ][ 0 ] - X[ i1 ][ 0 ],
                     X[ i2 ][ 1 ] - X[ i1 ][ 1 ],
                     X[ i2 ][ 2 ] - X[ i1 ][ 2 ] };
    real64 ac[3] = { X[ i3 ][ 0 ] - X[ i1 ][ 0 ],
                     X[ i3 ][ 1 ] - X[ i1 ][ 1 ],
                     X[ i3 ][ 2 ] - X[ i1 ][ 2 ] };
    real64 term1 = ab[1] * ac[2] - ab[2] * ac[1];
    real64 term2 = ab[2] * ac[0] - ab[0] * ac[2];
    real64 term3 = ab[0] * ac[1] - ab[1] * ac[0];
    return LvArray::math::sqrt( ( term1 * term1 + term2 * term2 + term3 * term3 ) );
 
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static void calcN( localIndex const q,
                     StackVariables const & stack,
                     real64 ( & N )[numNodes] )
  {
    GEOS_UNUSED_VAR( stack );
    return calcN( q, N );
  }
 
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr real64 EvaluateJacobiPolynomial( localIndex const n, real64 const alpha, real64 const beta, real64 const x )
  {
 
    if( n == 0 )
    {
      real64 val = 1.0;
      return val;
    }
 
    if( n == 1 )
    {
      real64 val = 0.5 * (alpha - beta + (alpha + beta + 2) * x);
      return val;
    }
 
    // Initial conditions for the recurrence relation
    real64 P_prev2 = 1.0;  // P_0
    real64 P_prev1 = 0.5 * (alpha - beta + (alpha + beta + 2) * x);  // P_1
    real64 P_current = 0.0;
 
    // Recurrence relation to compute the Jacobi polynomial
    for( int k = 2; k < n+1; ++k )
    {
      real64 a_k = 2 * k * (k + alpha + beta) * (2 * k + alpha + beta - 2);
      real64 b_k = (2 * k + alpha + beta - 1) * (alpha * alpha - beta * beta);
      real64 c_k = (2 * k + alpha + beta - 1) * (2 * k + alpha + beta) * (2 * k + alpha + beta - 2);
      real64 d_k = 2 * (k + alpha - 1) * (k + beta - 1) * (2 * k + alpha + beta);
 
      P_current = ((b_k + c_k * x) * P_prev1 - d_k * P_prev2) / a_k;
 
      // Mise à jour pour l'itération suivante
      P_prev2 = P_prev1;
      P_prev1 = P_current;
    }
 
    return P_current;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr real64 EvaluateJacobiPolynomialDerivative( localIndex const n, real64 const alpha, real64 const beta, real64 const x )
  {
    // Particular case for n = 0
    if( n == 0 )
    {
      return 0.0;
    }
 
    // dP_n^{(alpha,beta)}(x)/dx = (n+alpha+beta+1)/2 * P_{n-1}^{(alpha+1,beta+1)}(x)
    real64 coeff = 0.5 * (n + alpha + beta + 1.0);
    real64 jacobi_nm1 = EvaluateJacobiPolynomial( n-1, alpha+1.0, beta+1.0, x );
 
    return coeff * jacobi_nm1;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr real64 calcModal( localIndex const i, localIndex const j, localIndex const k, real64 const (&X)[3] )
  {
 
    real64 const x = X[0];
    real64 const y = X[1];
    real64 const z = X[2];
 
    real64 const epsilon = 1e-10; // Small value to avoid compairing floating point numbers directly
    if( LvArray::math::abs( z-1.0 ) < epsilon )
    {
      if( i == 0 && j == 0 )
      {
        return (k+2)*(k+1)/2.0;
      }
      else
      {
        return 0.0;
      }
    }
 
    real64 xi = x / (1.0 - z);
    real64 eta = y / (1.0 - z);
    real64 chi = 2.0 * z - 1.0;
    localIndex max_ij = LvArray::math::max( i, j );
    real64 P_i = EvaluateJacobiPolynomial( i, 0.0, 0.0, xi );
    real64 P_j = EvaluateJacobiPolynomial( j, 0.0, 0.0, eta );
    real64 P_k = EvaluateJacobiPolynomial( k, 2.0 * max_ij+2.0, 0.0, chi );
 
    return P_i * P_j * std::pow( 1.0 - z, max_ij ) * P_k;
 
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr void calcGradModal( int i, int j, int k, real64 const (&X)[3],
                                       real64 (& gradPsiX)[3] )
  {
    real64 x = X[0];
    real64 y = X[1];
    real64 z = X[2];
 
    real64 const epsilon = 1e-10; // Small value to avoid compairing floating point numbers directly
    if( LvArray::math::abs( z-1.0 ) < epsilon )
    {
      gradPsiX[0] = 0.0;
      gradPsiX[1] = 0.0;
      gradPsiX[2] = 0.0;
      return;
    }
 
 
    real64 xi = x / (1.0 - z);
    real64 eta = y / (1.0 - z);
    real64 chi = 2.0 * z - 1.0;
    localIndex m = LvArray::math::max( i, j );
    real64 P_i   = EvaluateJacobiPolynomial( i, 0.0, 0.0, xi );
    real64 P_j   = EvaluateJacobiPolynomial( j, 0.0, 0.0, eta );
    real64 P_k   = EvaluateJacobiPolynomial( k, 2.0 * m + 2.0, 0.0, chi );
    real64 dP_i  = EvaluateJacobiPolynomialDerivative( i, 0.0, 0.0, xi );
    real64 dP_j  = EvaluateJacobiPolynomialDerivative( j, 0.0, 0.0, eta );
    real64 dP_k  = EvaluateJacobiPolynomialDerivative( k, 2.0 * m + 2.0, 0.0, chi );
 
    real64 f = pow( 1.0 - z, m );
    // Gradient en x
    gradPsiX[0] = (1.0 / (1.0 - z)) * dP_i * P_j * f * P_k;
 
    // Gradient en y
    gradPsiX[1] = (1.0 / (1.0 - z)) * dP_j * P_i * f * P_k;
 
    // Gradient en z
    real64 dxi_dz  = x / std::pow( 1.0 - z, 2 );
    real64 deta_dz = y / std::pow( 1.0 - z, 2 );
    real64 dchi_dz = 2.0;
    real64 df_dz   = -m * std::pow( 1.0 - z, m - 1 );
    gradPsiX[2] =
      dP_i * dxi_dz * P_j * f * P_k +
      P_i * dP_j * deta_dz * f * P_k +
      P_i * P_j * df_dz * P_k +
      P_i * P_j * f * dP_k * dchi_dz;
 
 
  }
 
  GEOS_FORCE_INLINE
  static array2d< real64 > computeVanderMondeMatrixInverse( )
  {
    array2d< real64 > VDM;
    VDM.resize( numNodes, numNodes );
    real64 PsiX[numNodes] = {};
    real64 coords[numNodes][3] = {};
 
    generatePointsCoordinates( coords );
 
    for( int j = 0; j < numNodes; ++j )
    {
      localIndex count = 0;
      generateIndexes( [&]( localIndex const p, localIndex const r, localIndex const s )
      {
        PsiX[count]=calcModal( p, r, s, coords[j] );
        VDM[count][j] = PsiX[count];
        ++count;
      } );
      //   }
    }
    array2d< real64 > VDM_inv;
    VDM_inv.resize( numNodes, numNodes );
    // Inversion of VanDerMonde matrix
    BlasLapackLA::matrixInverse( VDM, VDM_inv );
    return VDM_inv;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr void calcN( localIndex q, arrayView2d< real64 const > VDM_inv, real64 (& N)[numNodes] )
  {
 
    real64 PsiX[numNodes] = {0.0};
    real64 coords[numNodes][3] = {{0.0}};
 
    //Generate the coordinates of the support points
    generatePointsCoordinates( coords );
 
 
    localIndex count = 0;
    //Compute the modal functions at the quadrature point
    generateIndexes( [&]( localIndex const p, localIndex const r, localIndex const s )
    {
      PsiX[count] = calcModal( p, r, s, coords[q] );
      ++count;
    } );
 
    //Compute the nodal functions at the quadrature point
    for( int i = 0; i < numNodes; ++i )
    {
      N[i] = 0.0;
      for( int j = 0; j < numNodes; ++j )
      {
        N[i] += VDM_inv[i][j] * PsiX[j];
      }
    }
 
  }
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr void calcN( real64 const (&X)[3], arrayView2d< real64 const > VDM_inv, real64 (& N)[numNodes] )
  {
 
    real64 PsiX[numNodes] = {};
 
    localIndex count = 0;
    //Compute the modal functions at the given point
    generateIndexes( [&]( localIndex const p, localIndex const r, localIndex const s )
    {
      PsiX[count] = calcModal( p, r, s, X );
      ++count;
    } );
 
    //Compute the nodal functions at the given points
    for( int i = 0; i < numNodes; ++i )
    {
      N[i] = 0.0;
      for( int j = 0; j < numNodes; ++j )
      {
        N[i] += VDM_inv[i][j] * PsiX[j];
      }
    }
 
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr void calcGradN( localIndex q, arrayView2d< real64 const > VDM_inv, real64 (& gradN)[numNodes][3] )
  {
 
    real64 gradModal[numNodes][3] = {{}};
    real64 coords[numNodes][3] = {};
 
    //Generate the coordinates of the support points
    generatePointsCoordinates( coords );
 
    localIndex count = 0;
    //Compute the modal functions derivatives at the quadrature point
    generateIndexes( [&]( localIndex const p, localIndex const r, localIndex const s )
    {
      calcGradModal( p, r, s, coords[q], gradModal[count] );
      ++count;
    } );
 
 
    //Compute the nodal functions derivatives at the quadrature point
    for( int i = 0; i < numNodes; ++i )
    {
      gradN[i][0] = 0.0;
      gradN[i][1] = 0.0;
      gradN[i][2] = 0.0;
      for( int j = 0; j < numNodes; ++j )
      {
        gradN[i][0] += VDM_inv[i][j] * gradModal[j][0];
        gradN[i][1] += VDM_inv[i][j] * gradModal[j][1];
        gradN[i][2] += VDM_inv[i][j] * gradModal[j][2];
      }
    }
 
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr void calcGradN( real64 const (&X)[3], arrayView2d< real64 const > VDM_inv, real64 (& gradN)[numNodes][3] )
  {
 
    real64 gradModal[numNodes][3] = {{}};
 
    localIndex count = 0;
    //Compute the modal functions derivatives at the given point
    generateIndexes( [&]( localIndex const p, localIndex const r, localIndex const s )
    {
      calcGradModal( p, r, s, X, gradModal[count] );
      ++count;
    } );
 
    //Compute the nodal functions derivatives at the given points
    for( int i = 0; i < numNodes; ++i )
    {
      gradN[i][0] = 0.0;
      gradN[i][1] = 0.0;
      gradN[i][2] = 0.0;
      for( int j = 0; j < numNodes; ++j )
      {
        gradN[i][0] += VDM_inv[i][j] * gradModal[j][0];
        gradN[i][1] += VDM_inv[i][j] * gradModal[j][1];
        gradN[i][2] += VDM_inv[i][j] * gradModal[j][2];
      }
    }
 
  }
 
  template< typename FUNC >
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static constexpr void computeMassTerm( arrayView2d< real64 const > VDM_inv, FUNC && func )
  {
 
    // Gauss-Legendre points and weights for the quadrature
    constexpr real64 GLeQuadraturePoints[3] = { -0.7745966692, 0.0, 0.7745966692 };
    constexpr real64 GLeQuadratureWeights[3] = { 0.5555555556, 0.8888888889, 0.5555555556 };
    // Gauss-Jacobi points and weights for the quadrature
    constexpr real64 GJQuadraturePoints[3] = { 0.07299, 0.34700, 0.70500 };
    constexpr real64 GJQuadratureWeights[3] = { 0.15714, 0.14625, 0.02995 };
 
 
 
    real64 N[numNodes] = {0.0};
    for( localIndex a = 0; a < numNodes; ++a )
    {
      for( localIndex b = 0; b < numNodes; ++b )
      {
        real64 val = 0.0;
        for( int i = 0; i < 3; ++i )
        {
          for( int j = 0; j < 3; ++j )
          {
            for( int k = 0; k < 3; ++k )
            {
 
              real64 xi = GLeQuadraturePoints[i];
              real64 eta = GLeQuadraturePoints[j];
              real64 chi = GJQuadraturePoints[k];
 
              real64 weight = GLeQuadratureWeights[i] * GLeQuadratureWeights[j] * GJQuadratureWeights[k];
              real64 x_i = (1-chi)*xi;
              real64 y_j = (1-chi)*eta;
              real64 z_k = chi;
              real64 X[3] = {x_i, y_j, z_k};
 
              calcN( X, VDM_inv, N );
 
              val += weight * N[a] * N[b];
 
            }
          }
        }
        func( a, b, val );
      }
    }
  }
 
  template< typename FUNC >
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static void computeStiffnessTerm( arrayView2d< real64 const > VDM_inv, FUNC && func )
  {
    // Gauss-Legendre points and weights for the quadrature
    constexpr real64 GLeQuadraturePoints[3] = { -0.7745966692, 0.0, 0.7745966692 };
    constexpr real64 GLeQuadratureWeights[3] = { 0.5555555556, 0.8888888889, 0.5555555556 };
    // Gauss-Jacobi points and weights for the quadrature
    constexpr real64 GJQuadraturePoints[3] = { 0.07299, 0.34700, 0.70500 };
    constexpr real64 GJQuadratureWeights[3] = { 0.15714, 0.14625, 0.02995 };
 
 
 
    real64 gradN[numNodes][3] = {{0.0}};
 
    for( localIndex a = 0; a < numNodes; ++a )
    {
      for( localIndex b = 0; b < numNodes; ++b )
      {
        real64 val = 0.0;
        for( int i = 0; i < 3; ++i )
        {
          for( int j = 0; j < 3; ++j )
          {
            for( int k = 0; k < 3; ++k )
            {
              real64 xi = GLeQuadraturePoints[i];
              real64 eta = GLeQuadraturePoints[j];
              real64 chi = GJQuadraturePoints[k];
 
              real64 weight = GLeQuadratureWeights[i] * GLeQuadratureWeights[j] * GJQuadratureWeights[k];
              real64 x_i = (1-chi)*xi;
              real64 y_j = (1-chi)*eta;
              real64 z_k = chi;
              real64 X[3] = {x_i, y_j, z_k};
 
              calcGradN( X, VDM_inv, gradN );
              real64 dot = gradN[a][0] * gradN[b][0]
                           + gradN[a][1] * gradN[b][1]
                           + gradN[a][2] * gradN[b][2];
              val += weight * dot;
            }
          }
        }
        func( a, b, val );
      }
    }
  }
 
};
 
using P1_Pyramid_BCD = Pk_Pyramid_BCD< 1 >;
 
 
}
}
 
#endif // GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_PKPYRAMIDBCD_HPP_