LagrangeBasis2.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_LAGRANGEBASIS2_HPP_
#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_ELEMENTFORMULATIONS_LAGRANGEBASIS2_HPP_
#include "common/DataTypes.hpp"
 
namespace geos
{
namespace finiteElement
{
 
class LagrangeBasis2
{
public:
  constexpr static localIndex numSupportPoints = 3;
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 weight( const int q )
  {
    switch( q )
    {
      case 0:
      case 2:
        return 1.0/3.0;
      default:
        return 4.0/3.0;
    }
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 parentSupportCoord( const localIndex supportPointIndex )
  {
    switch( supportPointIndex )
    {
      case 0:
        return -1.0;
        break;
      case 2:
        return 1.0;
      case 1:
      default:
        return 0.0;
    }
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 value( const int index,
                                 const real64 xi )
  {
 
    switch( index )
    {
      case 0:
        return value0( xi );
      case 2:
        return value2( xi );
      case 1:
      default:
        return value1( xi );
    }
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 value0( const real64 xi )
  {
    const real64 xi_div2 = 0.5 * xi;
    return -xi_div2 + xi_div2 * xi;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 value1( const real64 xi )
  {
    return 1.0 - xi * xi;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 value2( const real64 xi )
  {
    const real64 xi_div2 = 0.5 * xi;
    return xi_div2 + xi_div2 * xi;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient0( const real64 xi )
  {
    return -0.5 + xi;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient1( const real64 xi )
  {
    return -2 * xi;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient2( const real64 xi )
  {
    return 0.5 + xi;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradient( const int index,
                                    const real64 xi )
  {
    switch( index )
    {
      case 0:
        return gradient0( xi );
      case 2:
        return gradient2( xi );
      case 1:
      default:
        return gradient1( xi );
    }
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr static real64 gradientAt( const int q,
                                      const int p )
  {
    switch( q )
    {
      case 0:
        return p == 0 ? -1.5 : -0.5;
      case 1:
        return p == 0 ? 2.0 : 0.0;
      case 2:
        return p == 0 ? -0.5 : 0.5;
      default:
        return 0;
    }
  }
 
  struct TensorProduct2D
  {
 
    constexpr static localIndex numSupportPoints = 9;
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    constexpr static int linearIndex( const int i,
                                      const int j )
    {
      return i + 3 * j;
    }
 
 
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    constexpr static void multiIndex( const int linearIndex,
                                      int & i0,
                                      int & i1 )
    {
 
      i1 = ( ( linearIndex * 22 ) >> 6 );
      //i1 = a/3;
 
      i0 = linearIndex - i1 * 3;
    }
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    static void value( real64 const (&coords)[2],
                       real64 (& N)[numSupportPoints] )
    {
      for( int a=0; a<3; ++a )
      {
        for( int b=0; b<3; ++b )
        {
          const int lindex = LagrangeBasis2::TensorProduct2D::linearIndex( a, b );
          N[ lindex ] = LagrangeBasis2::value( a, coords[0] ) *
                        LagrangeBasis2::value( b, coords[1] );
        }
      }
    }
  };
 
  struct TensorProduct3D
  {
 
    constexpr static localIndex numSupportPoints = 27;
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    constexpr static int linearIndex( const int i,
                                      const int j,
                                      const int k )
    {
      return i + 3 * j + 9 * k;
    }
 
 
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    constexpr static void multiIndex( const int linearIndex,
                                      int & i0,
                                      int & i1,
                                      int & i2 )
    {
      i2 = ( linearIndex * 29 ) >> 8;
      //i2 = a/9;
 
      i1 = ( ( linearIndex * 22 ) >> 6 ) - i2 * 3;
      //i1 = a/3 - i2 * 3;
 
      i0 = linearIndex - i1 * 3 - i2 * 9;
    }
 
    GEOS_HOST_DEVICE
    GEOS_FORCE_INLINE
    static void value( const real64 (& coords)[3],
                       real64 (& N)[numSupportPoints] )
    {
      for( int a=0; a<3; ++a )
      {
        for( int b=0; b<3; ++b )
        {
          for( int c=0; c<3; ++c )
          {
            const int lindex = LagrangeBasis2::TensorProduct3D::linearIndex( a, b, c );
            N[ lindex ] = LagrangeBasis2::value( a, coords[0] ) *
                          LagrangeBasis2::value( b, coords[1] ) *
                          LagrangeBasis2::value( c, coords[2] );
          }
        }
      }
    }
  };
};
 
}
}
 
#endif /* GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_ELEMENTFORMULATIONS_LAGRANGEBASIS2_HPP_ */