H1_TriangleFace_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_H1TRIANGLEFACELAGRANGE1GAUSS_HPP_
#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS_HPP_
 
#include "FiniteElementBase.hpp"
#include "LagrangeBasis1.hpp"
 
 
namespace geos
{
namespace finiteElement
{
 
template< typename NUM_Q_POINTS >
class H1_TriangleFace_Lagrange1_Gauss_impl : public FiniteElementBase_impl< 3, 3, NUM_Q_POINTS::value >
{
public:
  using Base = FiniteElementBase_impl< 3, 3, NUM_Q_POINTS::value >;
 
  using StackVariables = typename Base::StackVariables;
 
  template< typename SubregionType >
  using MeshData = typename Base::template MeshData< SubregionType >;
 
  using Base::numNodes;
 
  using Base::numQuadraturePoints;
 
  using Base::maxSupportPoints;
 
  static_assert( ( NUM_Q_POINTS::value == 1 ||
                   NUM_Q_POINTS::value == 4 ||
                   NUM_Q_POINTS::value == 6 ),
                 "NUM_Q_POINTS::value must be 1, 4, or 6!" );
 
  using BASIS = LagrangeBasis1;
 
  H1_TriangleFace_Lagrange1_Gauss_impl() = default;
  ~H1_TriangleFace_Lagrange1_Gauss_impl() = default;
  H1_TriangleFace_Lagrange1_Gauss_impl( H1_TriangleFace_Lagrange1_Gauss_impl const & ) = default;
  H1_TriangleFace_Lagrange1_Gauss_impl & operator=( H1_TriangleFace_Lagrange1_Gauss_impl const & ) = default;
  H1_TriangleFace_Lagrange1_Gauss_impl( H1_TriangleFace_Lagrange1_Gauss_impl && ) = default;
  H1_TriangleFace_Lagrange1_Gauss_impl & operator=( H1_TriangleFace_Lagrange1_Gauss_impl && ) = default;
 
  GEOS_HOST_DEVICE
  localIndex getNumQuadraturePoints()
  {
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  static localIndex getNumQuadraturePoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  localIndex getNumSupportPoints()
  {
    return numNodes;
  }
 
  GEOS_HOST_DEVICE
  static localIndex getNumSupportPoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numNodes;
  }
 
  GEOS_HOST_DEVICE
  localIndex getMaxSupportPoints()
  {
    return maxSupportPoints;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  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] )
  {
 
    real64 const r  = pointCoord[0];
    real64 const s  = pointCoord[1];
 
    N[0] = (1.0 - r - s) * r * s;
 
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcBubbleN( localIndex const q,
                           real64 (& N)[1] )
  {
 
    real64 const qCoords[2] = {quadratureParentCoords0( q ), quadratureParentCoords1( 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] = 2;
  }
 
private:
  constexpr static real64 parentArea = 0.5;
 
  //constexpr static real64 weight = parentArea / numQuadraturePoints;
  constexpr static real64 weight = parentArea;
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureWeight( localIndex const q )
  {
 
    if constexpr (numQuadraturePoints == 1)
    {
      constexpr real64 w[numQuadraturePoints] = { 1.0 };
      return w[q];
    }
    else if constexpr (numQuadraturePoints == 4)
    {
      constexpr real64 w[numQuadraturePoints] = {-0.562500000000000,
                                                 0.520833333333333,
                                                 0.520833333333333,
                                                 0.520833333333333 };
      return w[q];
    }
    else if constexpr (numQuadraturePoints == 6)
    {
      real64 const w[numQuadraturePoints] = { 0.166666666666666,
                                              0.166666666666666,
                                              0.166666666666666,
                                              0.166666666666666,
                                              0.166666666666666,
                                              0.166666666666666 };
      return w[q];
    }
 
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords0( localIndex const q )
  {
 
    if constexpr (numQuadraturePoints == 1)
    {
      constexpr real64 qCoords[numQuadraturePoints] = { 1.0/3.0 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 4)
    {
      constexpr real64 qCoords[numQuadraturePoints] = { 0.333333333333333,
                                                        0.600000000000000,
                                                        0.200000000000000,
                                                        0.200000000000000 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 6)
    {
      constexpr real64 qCoords[numQuadraturePoints] = { 0.659027622374092,
                                                        0.109039009072877,
                                                        0.231933368553031,
                                                        0.659027622374092,
                                                        0.109039009072877,
                                                        0.231933368553031 };
      return qCoords[q];
    }
 
  }
 
  GEOS_HOST_DEVICE
  inline
  constexpr static real64 quadratureParentCoords1( localIndex const q )
  {
 
    if constexpr (numQuadraturePoints == 1)
    {
      constexpr real64 qCoords[numQuadraturePoints] = { 1.0/3.0 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 4)
    {
      constexpr real64 qCoords[numQuadraturePoints] = { 0.333333333333333,
                                                        0.200000000000000,
                                                        0.600000000000000,
                                                        0.200000000000000 };
      return qCoords[q];
    }
    else if constexpr (numQuadraturePoints == 6)
    {
      constexpr real64 qCoords[numQuadraturePoints] = { 0.231933368553031,
                                                        0.659027622374092,
                                                        0.109039009072877,
                                                        0.109039009072877,
                                                        0.231933368553031,
                                                        0.659027622374092 };
      return qCoords[q];
    }
 
  }
 
};
 
 
 
template< typename NUM_Q_POINTS >
template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, bool UPPER >
GEOS_HOST_DEVICE
inline
void H1_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS >::
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 );
}
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
void
H1_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS >::
calcN( real64 const (&coords)[2],
       real64 ( & N )[numNodes] )
{
  real64 const r  = coords[0];
  real64 const s  = coords[1];
 
  N[0] = 1.0 - r - s;
  N[1] = r;
  N[2] = s;
}
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
void
H1_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS >::
calcN( localIndex const q,
       real64 (& N)[numNodes] )
{
 
  real64 const qCoords[2] = {quadratureParentCoords0( q ), quadratureParentCoords1( q ) };
 
  calcN( qCoords, N );
 
}
 
template< typename NUM_Q_POINTS >
GEOS_HOST_DEVICE
inline
void H1_TriangleFace_Lagrange1_Gauss_impl< 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_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS >::
transformedQuadratureWeight( localIndex const q,
                             real64 const (&X)[numNodes][3] )
{
  //GEOS_UNUSED_VAR( q );
  real64 n[3] = { ( X[1][1] - X[0][1] ) * ( X[2][2] - X[0][2] ) - ( X[2][1] - X[0][1] ) * ( X[1][2] - X[0][2] ),
                  ( X[2][0] - X[0][0] ) * ( X[1][2] - X[0][2] ) - ( X[1][0] - X[0][0] ) * ( X[2][2] - X[0][2] ),
                  ( X[1][0] - X[0][0] ) * ( X[2][1] - X[0][1] ) - ( X[2][0] - X[0][0] ) * ( X[1][1] - X[0][1] )};
 
  return sqrt( n[0] * n[0] + n[1] * n[1] + n[2] * n[2] ) * weight * quadratureWeight( q );
}
 
 
 
template< typename NUM_Q_POINTS >
class H1_TriangleFace_Lagrange1_Gauss final : public H1_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS >,
  public FiniteElementBase
{
public:
 
  using BASIS = LagrangeBasis1;
 
  using ImplType = H1_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS >;
 
  H1_TriangleFace_Lagrange1_Gauss():
    FiniteElementBase( ImplType::numNodes,
                       ImplType::maxSupportPoints,
                       ImplType::numQuadraturePoints )
  {}
 
  virtual ~H1_TriangleFace_Lagrange1_Gauss() override final = default;
 
  H1_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS > * getImpl()
  {
    return static_cast< H1_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS > * >(this);
  }
 
  const H1_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS > * getImpl() const
  {
    return static_cast< H1_TriangleFace_Lagrange1_Gauss_impl< NUM_Q_POINTS > const * >(this);
  }
 
};
 
using H1_TriangleFace_Lagrange1_Gauss1 = H1_TriangleFace_Lagrange1_Gauss< std::integral_constant< int, 1 > >;
using H1_TriangleFace_Lagrange1_Gauss4 = H1_TriangleFace_Lagrange1_Gauss< std::integral_constant< int, 4 > >;
using H1_TriangleFace_Lagrange1_Gauss6 = H1_TriangleFace_Lagrange1_Gauss< std::integral_constant< int, 6 > >;
 
using H1_TriangleFace_Lagrange1_Gauss1_impl = H1_TriangleFace_Lagrange1_Gauss_impl< std::integral_constant< int, 1 > >;
using H1_TriangleFace_Lagrange1_Gauss4_impl = H1_TriangleFace_Lagrange1_Gauss_impl< std::integral_constant< int, 4 > >;
using H1_TriangleFace_Lagrange1_Gauss6_impl = H1_TriangleFace_Lagrange1_Gauss_impl< std::integral_constant< int, 6 > >;
 
}
}
#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS_HPP_