H1_QuadrilateralFace_Lagrange1_GaussLegendre2.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_H1QUADRILATERALFACELAGRANGE1GAUSSLEGENDRE2_HPP_
#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1QUADRILATERALFACELAGRANGE1GAUSSLEGENDRE2_HPP_
 
#include "FiniteElementBase.hpp"
#include "LagrangeBasis1.hpp"
 
 
namespace geos
{
namespace finiteElement
{
 
class H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl : public FiniteElementBase_impl< 4, 4, 4 >
{
public:
  using Base = FiniteElementBase_impl< 4, 4, 4 >;
 
  using BASIS = LagrangeBasis1;
 
  using StackVariables = typename Base::StackVariables;
 
  template< typename SubregionType >
  using MeshData = typename Base::template MeshData< SubregionType >;
 
  using Base::numNodes;
 
  using Base::numQuadraturePoints;
 
  using Base::maxSupportPoints;
 
 
  H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl() = default;
  ~H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl() = default;
  H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl( H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl const & ) = default;
  H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl & operator=( H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl const & ) = default;
  H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl( H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl && ) = default;
  H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl & operator=( H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl && ) = default;
 
 
  GEOS_HOST_DEVICE
  static localIndex getNumQuadraturePoints()
  {
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  static localIndex getNumQuadraturePoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  static localIndex getNumSupportPoints()
  {
    return numNodes;
  }
 
  GEOS_HOST_DEVICE
  static localIndex getMaxSupportPoints()
  {
    return maxSupportPoints;
  }
 
  GEOS_HOST_DEVICE
  static localIndex getNumSupportPoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numNodes;
  }
 
  GEOS_HOST_DEVICE
  static void calcN( real64 const (&coords)[2],
                     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 calcBubbleN( real64 const (&pointCoord)[2],
                           real64 (& N)[1] )
  {
    LagrangeBasis1::TensorProduct2D::valueBubble( pointCoord, N );
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcBubbleN( localIndex const q,
                           real64 (& N)[1] )
  {
    real64 const qCoords[2] = { quadratureFactor *parentCoords0( q ),
                                quadratureFactor *parentCoords1( q ) };
 
    calcBubbleN( qCoords, N );
  }
 
  GEOS_HOST_DEVICE
  static real64 transformedQuadratureWeight( localIndex const q,
                                             real64 const (&X)[numNodes][3] );
 
  template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, bool UPPER >
  GEOS_HOST_DEVICE
  inline
  static void addGradGradStabilization( StackVariables const & stack,
                                        real64 ( &matrix )
                                        [maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT]
                                        [maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT],
                                        real64 const & scaleFactor );
 
  GEOS_HOST_DEVICE
  inline
  static void getPermutation( int (& permutation)[numNodes] )
  {
    permutation[0] = 0;
    permutation[1] = 1;
    permutation[2] = 3;
    permutation[3] = 2;
  }
 
private:
  constexpr static real64 parentArea = 4.0;
 
  constexpr static real64 weight = parentArea / numQuadraturePoints;
 
  constexpr static real64 quadratureFactor = 1.0 / 1.732050807568877293528;
 
  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);
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 parentCoords1( localIndex const a )
  {
    return -1.0 + ( a & 2 );
  }
 
};
 
 
template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, bool UPPER >
GEOS_HOST_DEVICE
inline
void H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl::
  addGradGradStabilization( StackVariables const & stack,
                            real64 ( & matrix )
                            [maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT]
                            [maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT],
                            real64 const & scaleFactor )
{
  GEOS_UNUSED_VAR( stack );
  GEOS_UNUSED_VAR( matrix );
  GEOS_UNUSED_VAR( scaleFactor );
}
 
GEOS_HOST_DEVICE
inline
void
H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl::
  calcN( real64 const (&coords)[2],
         real64 (& N)[numNodes] )
{
  for( localIndex a=0; a<numNodes; ++a )
  {
    N[a] = 0.25 *
           ( 1 + quadratureFactor*coords[0]*parentCoords0( a ) ) *
           ( 1 + quadratureFactor*coords[1]*parentCoords1( a ) );
  }
}
GEOS_HOST_DEVICE
inline
void
H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl::
  calcN( localIndex const q,
         real64 (& N)[numNodes] )
{
  for( localIndex a=0; a<numNodes; ++a )
  {
    N[a] = 0.25 *
           ( 1 + quadratureFactor*parentCoords0( q )*parentCoords0( a ) ) *
           ( 1 + quadratureFactor*parentCoords1( q )*parentCoords1( a ) );
  }
}
 
GEOS_HOST_DEVICE
inline
void H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl::
  calcN( localIndex const q,
         StackVariables const & GEOS_UNUSED_PARAM( stack ),
         real64 ( & N )[numNodes] )
{
  return calcN( q, N );
}
 
//*************************************************************************************************
 
GEOS_HOST_DEVICE
inline
real64
H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl::
  transformedQuadratureWeight( localIndex const q,
                               real64 const (&X)[numNodes][3] )
{
  real64 dXdXi[3][2] = {{0}};
 
  real64 const quadratureCoords[2] = { quadratureFactor *parentCoords0( q ),
                                       quadratureFactor *parentCoords1( 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] ) };
  constexpr real64 dpsi[2] = { -0.5, 0.5 };
 
  for( int a=0; a<2; ++a )
  {
    for( int b=0; b<2; ++b )
    {
      real64 const dNdXi[2] = { dpsi[a] * psi1[b],
                                psi0[a] * dpsi[b] };
 
      localIndex const nodeIndex = linearMap( a, b );
 
      for( int i = 0; i < 3; ++i )
      {
        for( int j = 0; j < 2; ++j )
        {
          dXdXi[i][j] = dXdXi[i][j] + dNdXi[ j ] * X[nodeIndex][i];
        }
      }
    }
  }
 
  real64 const detJ = pow( dXdXi[1][0] * dXdXi[2][1] - dXdXi[2][0] * dXdXi[1][1], 2.0 )
                      + pow( dXdXi[2][0] * dXdXi[0][1] - dXdXi[0][0] * dXdXi[2][1], 2.0 )
                      + pow( dXdXi[0][0] * dXdXi[1][1] - dXdXi[1][0] * dXdXi[0][1], 2.0 );
 
  return sqrt( detJ ) * weight;
}
 
 
 
class H1_QuadrilateralFace_Lagrange1_GaussLegendre2 final : public H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl,
  public FiniteElementBase
{
public:
 
  using BASIS = LagrangeBasis1;
 
  using ImplType = H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl;
 
  H1_QuadrilateralFace_Lagrange1_GaussLegendre2():
    FiniteElementBase( ImplType::numNodes,
                       ImplType::maxSupportPoints,
                       ImplType::numQuadraturePoints )
  {}
 
  virtual ~H1_QuadrilateralFace_Lagrange1_GaussLegendre2() override final = default;
 
  H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl * getImpl()
  {
    return static_cast< H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl * >(this);
  }
 
  H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl const * getImpl() const
  {
    return static_cast< H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl const * >(this);
  }
 
 
};
 
}
}
#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1QUADRILATERALFACELAGRANGE1GAUSSLEGENDRE2_HPP_