LagrangeBasis4GL.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) 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_ELEMENTFORMULATIONS_LAGRANGEBASIS4GL_HPP_
#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_ELEMENTFORMULATIONS_LAGRANGEBASIS4GL_HPP_
 
#include "common/DataTypes.hpp"
 
namespace geos
{
namespace finiteElement
{
 
class LagrangeBasis4GL
{
public:
  constexpr static localIndex numSupportPoints = 5;
 
  constexpr static real64 sqrt3_7 = 0.6546536707079771;
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 weight( const int q )
  {
    switch( q )
    {
      case 0:
      case 4:
        return 1.0/10.0;
      case 1:
      case 3:
        return 49.0/90.0;
      default:
        return 32.0/45.0;
    }
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  // depending on the supportPointIndex value
  //Switch case
  constexpr static real64 parentSupportCoord( const localIndex supportPointIndex )
  {
    real64 result=0.0;
 
    switch( supportPointIndex )
    {
      case 0:
        result = -1.0;
        break;
      case 1:
        result = -sqrt3_7;
        break;
      case 2:
        result = 0.0;
        break;
      case 3:
        result = sqrt3_7;
        break;
      case 4:
        result = 1.0;
        break;
      default:
        break;
    }
 
    return result;
  }
 
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 value( const int index,
                                 const real64 xi )
  {
    real64 result=0.0;
 
    switch( index )
    {
      case 0:
        result = LagrangeBasis4GL::value0( xi );
        break;
      case 1:
        result = LagrangeBasis4GL::value1( xi );
        break;
      case 2:
        result = LagrangeBasis4GL::value2( xi );
        break;
      case 3:
        result = LagrangeBasis4GL::value3( xi );
        break;
      case 4:
        result = LagrangeBasis4GL::value4( xi );
        break;
      default:
        break;
    }
 
    return result;
  }
 
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  //MODFI4 : Implemented new base functions and their derivative for Q3
  constexpr static real64 value0( const real64 xi )
  {
    return (1.0/8.0)*(-1.0+xi)*xi*(-3.0+7.0*xi*xi);
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 value1( const real64 xi )
  {
    return (49.0/24.0)*(sqrt3_7-xi)*xi*(-1.0+xi*xi);
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 value2( const real64 xi )
  {
    return (1.0/3.0)*(3.0-10.0*xi*xi+7.0*xi*xi*xi*xi);
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 value3( const real64 xi )
  {
    return -(49.0/24.0)*(sqrt3_7+xi)*xi*(-1.0+xi*xi);
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 value4( const real64 xi )
  {
    return (1.0/8.0)*(1.0+xi)*xi*(-3.0+7.0*xi*xi);
  }
 
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient( const int index,
                                    const real64 xi )
  {
    real64 result=0.0;
 
    switch( index )
    {
      case 0:
        result = LagrangeBasis4GL::gradient0( xi );
        break;
      case 1:
        result = LagrangeBasis4GL::gradient1( xi );
        break;
      case 2:
        result = LagrangeBasis4GL::gradient2( xi );
        break;
      case 3:
        result = LagrangeBasis4GL::gradient3( xi );
        break;
      case 4:
        result = LagrangeBasis4GL::gradient4( xi );
        break;
      default:
        break;
    }
 
    return result;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient0( const real64 xi )
  {
    return (1.0/8.0)*(3.0+xi*(-6.0+7.0*xi*(-3.0+4.0*xi)));
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient1( const real64 xi )
  {
    return (49.0/24.0)*(-sqrt3_7+xi*(2.0+3.0*sqrt3_7*xi-4.0*xi*xi));
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient2( const real64 xi )
  {
    return (4.0/3.0)*xi*(-5.0+7.0*xi*xi);
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient3( const real64 xi )
  {
    return (49.0/24.0)*(sqrt3_7+xi*(2.0-3.0*sqrt3_7*xi-4.0*xi*xi));
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient4( const real64 xi )
  {
    return (1.0/8.0)*(-3.0+xi*(-6.0+7.0*xi*(3.0+4.0*xi)));
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradientAt( const int q,
                                      const int p )
  {
    switch( q )
    {
      case 0:
        switch( p )
        {
          case 0: return -5.0000000000000000000;
          case 1: return -1.2409902530309828578;
          case 2: return 0.37500000000000000000;
        }
        break;
      case 1:
        switch( p )
        {
          case 0: return 6.7565024887242400038;
          case 1: return 0.0;
          case 2: return -1.3365845776954533353;
        }
        break;
      case 2:
        switch( p )
        {
          case 0: return -2.6666666666666666667;
          case 1: return 1.7457431218879390501;
          case 2: return 0.0;
        }
        break;
      case 3:
        switch( p )
        {
          case 0: return 1.4101641779424266628;
          case 1: return -0.7637626158259733344;
          case 2: return 1.3365845776954533353;
        }
        break;
      case 4:
        switch( p )
        {
          case 0: return -0.50000000000000000000;
          case 1: return 0.25900974696901714215;
          case 2: return -0.37500000000000000000;
        }
        break;
    }
    return 0;
  }
 
  struct TensorProduct2D
  {
    constexpr static localIndex numSupportPoints = 25;
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    constexpr static int linearIndex( const int i,
                                      const int j )
    {
      return i + 5 * j;
    }
 
 
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    constexpr static void multiIndex( int const linearIndex,
                                      int & i0,
                                      int & i1 )
    {
 
      i1 = linearIndex/5;
 
      i0 = linearIndex%5;
 
    }
 
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    static void value( const real64 (& coords)[2],
                       real64 (& N)[numSupportPoints] )
    {
      for( int a=0; a<5; ++a )
      {
        for( int b=0; b<5; ++b )
        {
          const int lindex = LagrangeBasis4GL::TensorProduct2D::linearIndex( a, b );
          N[ lindex ] = LagrangeBasis4GL::value( a, coords[0] ) *
                        LagrangeBasis4GL::value( b, coords[1] );
        }
      }
    }
  };
 
  struct TensorProduct3D
  {
    constexpr static localIndex numSupportPoints = 125;
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    constexpr static int linearIndex( const int i,
                                      const int j,
                                      const int k )
    {
      return i + 5 * j + 25 * k;
    }
 
 
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    constexpr static void multiIndex( int const linearIndex,
                                      int & i0,
                                      int & i1,
                                      int & i2 )
    {
 
      i2 = linearIndex/25;
 
      i1 = (linearIndex%25)/5;
 
      i0 = (linearIndex%25)%5;
 
    }
 
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    static void value( const real64 (& coords)[3],
                       real64 (& N)[numSupportPoints] )
    {
      for( int a=0; a<5; ++a )
      {
        for( int b=0; b<5; ++b )
        {
          for( int c=0; c<5; ++c )
          {
            const int lindex = LagrangeBasis4GL::TensorProduct3D::linearIndex( a, b, c );
            N[ lindex ] = LagrangeBasis4GL::value( a, coords[0] ) *
                          LagrangeBasis4GL::value( b, coords[1] ) *
                          LagrangeBasis4GL::value( c, coords[2] );
          }
        }
      }
    }
  };
};
 
 
}
}
 
 
#endif /* GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_ELEMENTFORMULATIONS_LAGRANGEBASIS4GL_HPP_ */