H1_Hexahedron_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_TRILINEARHEXAHEDRON_HPP_
#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_TRILINEARHEXAHEDRON_HPP_
 
#include "FiniteElementBase.hpp"
#include "LagrangeBasis1.hpp"
 
#include <utility>
 
 
 
namespace geos
{
namespace finiteElement
{
 
//constexpr static real64 linearBasisAtQuadrature[2] = { 0.5 + 0.5 * 0.5773502691896257645092,
//                                                       0.5 - 0.5 * 0.5773502691896257645092 };
//__constant__ real64 psiProduct[3] = { 0.5 * linearBasisAtQuadrature[0]*linearBasisAtQuadrature[0],
//                                          0.5 * linearBasisAtQuadrature[0]*linearBasisAtQuadrature[1],
//                                          0.5 * linearBasisAtQuadrature[1]*linearBasisAtQuadrature[1] };
 
//__constant__ real64 psiProduct[3] = {  0.311004233964073108,
//                                       0.083333333333333333,
//                                       0.022329099369260226 };
//
//__constant__ short dpsi[2] = { -1, 1 };
 
class H1_Hexahedron_Lagrange1_GaussLegendre2 final : public FiniteElementBase
{
public:
  constexpr static localIndex numNodes = LagrangeBasis1::TensorProduct3D::numSupportPoints;
 
  constexpr static localIndex numFaces = LagrangeBasis1::TensorProduct3D::numSupportFaces;
 
  constexpr static localIndex maxSupportPoints = numNodes;
 
  constexpr static localIndex numQuadraturePoints = 8;
 
  constexpr static int numSamplingPoints = numSamplingPointsPerDirection * numSamplingPointsPerDirection * numSamplingPointsPerDirection;
 
  USING_FINITEELEMENTBASE
  GEOS_HOST_DEVICE
  virtual ~H1_Hexahedron_Lagrange1_GaussLegendre2() override
  {}
 
  GEOS_HOST_DEVICE
  virtual localIndex getNumQuadraturePoints() const override
  {
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  static localIndex getNumQuadraturePoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numQuadraturePoints;
  }
 
  GEOS_HOST_DEVICE
  virtual localIndex getNumSupportPoints() const override
  {
    return numNodes;
  }
 
  GEOS_HOST_DEVICE
  virtual localIndex getMaxSupportPoints() const override
  {
    return maxSupportPoints;
  }
 
  GEOS_HOST_DEVICE
  static localIndex getNumSupportPoints( StackVariables const & stack )
  {
    GEOS_UNUSED_VAR( stack );
    return numNodes;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static void getSamplingPointCoordInParentSpace( int const linearIndex,
                                                  real64 (& samplingPointCoord)[3] )
  {
    const int i0 = linearIndex % numSamplingPointsPerDirection;
    const int i1 = ( (linearIndex - i0)/numSamplingPointsPerDirection ) % numSamplingPointsPerDirection;
    const int i2 = ( (linearIndex - i0)/numSamplingPointsPerDirection - i1 ) / numSamplingPointsPerDirection;
 
    constexpr real64 step = 2. / ( numSamplingPointsPerDirection - 1 );
 
    samplingPointCoord[0] = -1 + i0 * step;
    samplingPointCoord[1] = -1 + i1 * step;
    samplingPointCoord[2] = -1 + i2 * step;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  static void calcN( real64 const (&pointCoord)[3],
                     real64 (& N)[numNodes] )
  {
    LagrangeBasis1::TensorProduct3D::value( pointCoord, N );
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcN( localIndex const q,
                     real64 (& N)[numNodes] )
  {
    int qa, qb, qc;
    LagrangeBasis1::TensorProduct3D::multiIndex( q, qa, qb, qc );
    real64 const qCoords[3] = { quadratureFactor * LagrangeBasis1::parentSupportCoord( qa ),
                                quadratureFactor * LagrangeBasis1::parentSupportCoord( qb ),
                                quadratureFactor * LagrangeBasis1::parentSupportCoord( qc ) };
 
    calcN( qCoords, N );
  }
 
  GEOS_HOST_DEVICE
  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 void calcFaceBubbleN( real64 const (&pointCoord)[3],
                               real64 (& N)[numFaces] )
  {
    LagrangeBasis1::TensorProduct3D::valueFaceBubble( pointCoord, N );
  }
 
  GEOS_HOST_DEVICE
  inline
  static void calcFaceBubbleN( localIndex const q,
                               real64 (& N)[numFaces] )
  {
    int qa, qb, qc;
    LagrangeBasis1::TensorProduct3D::multiIndex( q, qa, qb, qc );
    real64 const qCoords[3] = { quadratureFactor * LagrangeBasis1::parentSupportCoord( qa ),
                                quadratureFactor * LagrangeBasis1::parentSupportCoord( qb ),
                                quadratureFactor * LagrangeBasis1::parentSupportCoord( qc ) };
 
    calcFaceBubbleN( qCoords, N );
  }
 
  GEOS_HOST_DEVICE
  static real64 calcGradN( localIndex const q,
                           real64 const (&X)[numNodes][3],
                           real64 ( &gradN )[numNodes][3] );
 
  GEOS_HOST_DEVICE
  inline
  static real64 calcGradN( localIndex const q,
                           real64 const (&X)[numNodes][3],
                           StackVariables const & stack,
                           real64 ( &gradN )[numNodes][3] );
 
  GEOS_HOST_DEVICE
  static real64 calcGradFaceBubbleN( localIndex const q,
                                     real64 const (&X)[numNodes][3],
                                     real64 ( &gradN )[numFaces][3] );
 
  GEOS_HOST_DEVICE
  static real64 transformedQuadratureWeight( localIndex const q,
                                             real64 const (&X)[numNodes][3] );
 
 
  GEOS_HOST_DEVICE
  static real64 transformedQuadratureWeight( localIndex const q,
                                             real64 const (&X)[numNodes][3],
                                             StackVariables const & stack )
  { GEOS_UNUSED_VAR( stack ); return transformedQuadratureWeight( q, X ); }
 
 
  GEOS_HOST_DEVICE
  static real64 invJacobianTransformation( int const q,
                                           real64 const (&X)[numNodes][3],
                                           real64 ( & J )[3][3] )
  {
    int qa, qb, qc;
    LagrangeBasis1::TensorProduct3D::multiIndex( q, qa, qb, qc );
    jacobianTransformation( qa, qb, qc, X, J );
    return LvArray::tensorOps::invert< 3 >( J );
  }
 
 
  GEOS_HOST_DEVICE
  static void symmetricGradient( int const q,
                                 real64 const (&invJ)[3][3],
                                 real64 const (&var)[numNodes][3],
                                 real64 ( &grad )[6] );
 
 
 
  GEOS_HOST_DEVICE
  static void gradient( int const q,
                        real64 const (&invJ)[3][3],
                        real64 const (&var)[numNodes][3],
                        real64 ( &grad )[3][3] );
 
 
  GEOS_HOST_DEVICE
  static void plusGradNajAij( int const q,
                              real64 const (&invJ)[3][3],
                              real64 const (&var)[6],
                              real64 ( &R )[numNodes][3] );
 
 
 
  GEOS_HOST_DEVICE
  static void jacobianTransformation( int const qa,
                                      int const qb,
                                      int const qc,
                                      real64 const (&X)[numNodes][3],
                                      real64 ( &J )[3][3] );
 
 
  GEOS_HOST_DEVICE
  static void
    applyTransformationToParentGradients( int const qa,
                                          int const qb,
                                          int const qc,
                                          real64 const ( &invJ )[3][3],
                                          real64 ( &gradN )[numNodes][3] );
 
 
private:
  constexpr static real64 parentLength = LagrangeBasis1::parentSupportCoord( 1 ) - LagrangeBasis1::parentSupportCoord( 0 );
 
  constexpr static real64 parentVolume = parentLength*parentLength*parentLength;
 
  constexpr static real64 weight = parentVolume / numQuadraturePoints;
 
  constexpr static real64 quadratureFactor = 1.0 / 1.732050807568877293528;
 
//  constexpr static real64 psiProduct0 = 0.311004233964073108;
//  constexpr static real64 psiProduct1 = 0.083333333333333333;
//  constexpr static real64 psiProduct2 = 0.022329099369260226;
//
//  constexpr static short dpsi0 = -1;
//  constexpr static short dpsi1 = 1;
 
 
  template< typename FUNC, typename ... PARAMS >
  GEOS_HOST_DEVICE
  static void supportLoop( int const qa,
                           int const qb,
                           int const qc,
                           FUNC && func,
                           PARAMS &&... params );
 
};
 
//GEOS_HOST_DEVICE inline real64
//psiProductFunc( int const n )
//{
//  // factor = 1./(2 + Sqrt[3]) = 0.267949192431122706
//
//  real64 rval = 0.311004233964073108;
//  for( int a=0; a<n; ++a )
//  {
//    rval *= 0.267949192431122706;
//  }
//  return rval;
//}
 
 
template< typename FUNC, typename ... PARAMS >
GEOS_HOST_DEVICE inline void
H1_Hexahedron_Lagrange1_GaussLegendre2::supportLoop( int const qa,
                                                     int const qb,
                                                     int const qc,
                                                     FUNC && func,
                                                     PARAMS &&... params )
{
 
  #define PARENT_GRADIENT_METHOD 2
#if PARENT_GRADIENT_METHOD == 1
  // This option calculates the basis values at the quadrature point for each
  // linear basis index.
 
  real64 const quadratureCoords[3] = { -quadratureFactor + 1.154700538379252 * qa,
                                       -quadratureFactor + 1.154700538379252 * qb,
                                       -quadratureFactor + 1.154700538379252 * qc };
 
  real64 const psi0[2] = { 0.5 - 0.5 * quadratureCoords[0],
                           0.5 + 0.5 * quadratureCoords[0] };
  real64 const psi1[2] = { 0.5 - 0.5 * quadratureCoords[1],
                           0.5 + 0.5 * quadratureCoords[1] };
  real64 const psi2[2] = { 0.5 - 0.5 * quadratureCoords[2],
                           0.5 + 0.5 * quadratureCoords[2] };
  constexpr real64 dpsi[2] = { -0.5, 0.5 };
#elif PARENT_GRADIENT_METHOD == 2
  // This option calculates the product of linear basis prior to use.
  // The tensor product basis gradient may be expressed as a permutation of the
  // product between the two possible linear basis gradients. Thus the values
  // in the basis loop of qaa/qbb/qcc indicate which permutation to choose.
  // The quantities qaa, qbb, qcc are the difference in index between the
  // quadrature point and the basis. This is possible because there are 8
  // basis, and 8 quadrature points, which have correlated indices. So the
  // values of qaa/qbb/qcc are the "distance" from the quadrature point index
  // and the support point index.
  // THIS approach uses about 10 less registers than option 1, with no apparent
  // cost.
 
//  constexpr static real64 linearBasisAtQuadrature[2] = { 0.5 + 0.5 * quadratureFactor,
//                                                         0.5 - 0.5 * quadratureFactor };
//  constexpr static real64 psiProduct[3] = { 0.5 * linearBasisAtQuadrature[0]*linearBasisAtQuadrature[0],
//                                            0.5 * linearBasisAtQuadrature[0]*linearBasisAtQuadrature[1],
//                                            0.5 * linearBasisAtQuadrature[1]*linearBasisAtQuadrature[1] };
 
  constexpr static real64 psiProduct[3] = { 0.311004233964073108, 0.083333333333333333, 0.022329099369260226};
  constexpr static int dpsi[2] = { -1, 1 };
 
//  constexpr static real64 psiProduct[3] = { psiProduct0, psiProduct1, psiProduct2 };
//  constexpr short dpsi[2] = { dpsi0, dpsi1 };
#endif
 
  // Loop over the linear basis indices in each direction.
  for( int a=0; a<2; ++a )
  {
#if PARENT_GRADIENT_METHOD == 2
    int const qaa = ( a^qa ); // abs(a-qa)
#endif
    for( int b=0; b<2; ++b )
    {
#if PARENT_GRADIENT_METHOD == 2
      int const qbb = ( b^qb );
#endif
      for( int c=0; c<2; ++c )
      {
#if PARENT_GRADIENT_METHOD == 2
        int const qcc = ( c^qc );
#endif
 
#if PARENT_GRADIENT_METHOD == 1
        real64 const dNdXi[3] = { dpsi[a] * psi1[b] * psi2[c],
                                  psi0[a] * dpsi[b] * psi2[c],
                                  psi0[a] * psi1[b] * dpsi[c] };
#elif PARENT_GRADIENT_METHOD == 2
 
 
//        const real64 dNdXi[3] = { dpsi[a] * psiProductFunc( qbb + qcc ),
//                                  dpsi[b] * psiProductFunc( qaa + qcc ),
//                                  dpsi[c] * psiProductFunc( qaa + qbb ) };
 
        real64 const dNdXi[3] = { dpsi[a] * psiProduct[ qbb + qcc ],
                                  dpsi[b] * psiProduct[ qaa + qcc ],
                                  dpsi[c] * psiProduct[ qaa + qbb ] };
#endif
        localIndex const nodeIndex = LagrangeBasis1::TensorProduct3D::linearIndex( a, b, c );
 
        func( dNdXi, nodeIndex, std::forward< PARAMS >( params )... );
      }
    }
  }
}
 
//*************************************************************************************************
GEOS_HOST_DEVICE
inline
real64
H1_Hexahedron_Lagrange1_GaussLegendre2::calcGradN( localIndex const q,
                                                   real64 const (&X)[numNodes][3],
                                                   real64 (& gradN)[numNodes][3] )
{
  real64 J[3][3] = {{0}};
 
 
  int qa, qb, qc;
  LagrangeBasis1::TensorProduct3D::multiIndex( q, qa, qb, qc );
 
  jacobianTransformation( qa, qb, qc, X, J );
 
  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
 
  applyTransformationToParentGradients( qa, qb, qc, J, gradN );
 
  return detJ * weight;
}
 
GEOS_HOST_DEVICE
inline
real64 H1_Hexahedron_Lagrange1_GaussLegendre2::
  calcGradN( localIndex const q,
             real64 const (&X)[numNodes][3],
             StackVariables const & GEOS_UNUSED_PARAM( stack ),
             real64 ( & gradN )[numNodes][3] )
{
  return calcGradN( q, X, gradN );
}
 
GEOS_HOST_DEVICE
inline
real64
H1_Hexahedron_Lagrange1_GaussLegendre2::calcGradFaceBubbleN( localIndex const q,
                                                             real64 const (&X)[numNodes][3],
                                                             real64 (& gradN)[numFaces][3] )
{
  real64 J[3][3] = {{0}};
 
 
  int qa, qb, qc;
  LagrangeBasis1::TensorProduct3D::multiIndex( q, qa, qb, qc );
 
  jacobianTransformation( qa, qb, qc, X, J );
 
  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
 
  real64 dNdXi[numFaces][3] = {{0}};
 
  real64 const qCoords[3] = { quadratureFactor * LagrangeBasis1::parentSupportCoord( qa ),
                              quadratureFactor * LagrangeBasis1::parentSupportCoord( qb ),
                              quadratureFactor * LagrangeBasis1::parentSupportCoord( qc ) };
 
  LagrangeBasis1::TensorProduct3D::gradientFaceBubble( qCoords, dNdXi );
 
  for( int fi=0; fi<numFaces; ++fi )
  {
    gradN[fi][0] = dNdXi[fi][0] * J[0][0] + dNdXi[fi][1] * J[1][0] + dNdXi[fi][2] * J[2][0];
    gradN[fi][1] = dNdXi[fi][0] * J[0][1] + dNdXi[fi][1] * J[1][1] + dNdXi[fi][2] * J[2][1];
    gradN[fi][2] = dNdXi[fi][0] * J[0][2] + dNdXi[fi][1] * J[1][2] + dNdXi[fi][2] * J[2][2];
  }
 
  return detJ * weight;
}
 
//*************************************************************************************************
#if __GNUC__
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wshadow"
#endif
 
GEOS_HOST_DEVICE
inline
void
H1_Hexahedron_Lagrange1_GaussLegendre2::
  jacobianTransformation( int const qa,
                          int const qb,
                          int const qc,
                          real64 const (&X)[numNodes][3],
                          real64 ( & J )[3][3] )
{
  supportLoop( qa, qb, qc, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
                                                 int const nodeIndex,
                                                 real64 const (&X)[numNodes][3],
                                                 real64 (& J)[3][3] )
  {
    real64 const * const GEOS_RESTRICT Xnode = X[nodeIndex];
    for( int i = 0; i < 3; ++i )
    {
      for( int j = 0; j < 3; ++j )
      {
        J[i][j] = J[i][j] + dNdXi[ j ] * Xnode[i];
      }
    }
 
//    J[0][0] = J[0][0] + dNdXi[0] * Xnode[0];
//    J[0][1] = J[0][1] + dNdXi[1] * Xnode[0];
//    J[0][2] = J[0][2] + dNdXi[2] * Xnode[0];
//    J[1][0] = J[1][0] + dNdXi[0] * Xnode[1];
//    J[1][1] = J[1][1] + dNdXi[1] * Xnode[1];
//    J[1][2] = J[1][2] + dNdXi[2] * Xnode[1];
//    J[2][0] = J[2][0] + dNdXi[0] * Xnode[2];
//    J[2][1] = J[2][1] + dNdXi[1] * Xnode[2];
//    J[2][2] = J[2][2] + dNdXi[2] * Xnode[2];
 
  }, X, J );
}
 
 
//*************************************************************************************************
GEOS_HOST_DEVICE
inline
void
H1_Hexahedron_Lagrange1_GaussLegendre2::
  applyTransformationToParentGradients( int const qa,
                                        int const qb,
                                        int const qc,
                                        real64 const ( &invJ )[3][3],
                                        real64 (& gradN)[numNodes][3] )
{
  supportLoop( qa, qb, qc, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
                                                 int const nodeIndex,
                                                 real64 const (&invJ)[3][3],
                                                 real64 (& gradN)[numNodes][3] )
  {
//    for( int i = 0; i < 3; ++i )
//    {
//      gradN[nodeIndex][i] = 0.0;
//      for( int j = 0; j < 3; ++j )
//      {
//        gradN[nodeIndex][i] = gradN[nodeIndex][i] + dNdXi[ j ] * invJ[j][i];
//      }
//    }
    // smaller register footprint by manually unrolling the for loops.
    gradN[nodeIndex][0] = dNdXi[0] * invJ[0][0] + dNdXi[1] * invJ[1][0] + dNdXi[2] * invJ[2][0];
    gradN[nodeIndex][1] = dNdXi[0] * invJ[0][1] + dNdXi[1] * invJ[1][1] + dNdXi[2] * invJ[2][1];
    gradN[nodeIndex][2] = dNdXi[0] * invJ[0][2] + dNdXi[1] * invJ[1][2] + dNdXi[2] * invJ[2][2];
 
 
  }, invJ, gradN );
}
 
 
GEOS_HOST_DEVICE
inline
real64
H1_Hexahedron_Lagrange1_GaussLegendre2::
  transformedQuadratureWeight( localIndex const q,
                               real64 const (&X)[numNodes][3] )
{
  real64 J[3][3] = {{0}};
 
  int qa, qb, qc;
  LagrangeBasis1::TensorProduct3D::multiIndex( q, qa, qb, qc );
 
  jacobianTransformation( qa, qb, qc, X, J );
 
  return LvArray::tensorOps::determinant< 3 >( J );
}
 
 
 
GEOS_HOST_DEVICE
inline
void H1_Hexahedron_Lagrange1_GaussLegendre2::symmetricGradient( int const q,
                                                                real64 const (&invJ)[3][3],
                                                                real64 const (&var)[numNodes][3],
                                                                real64 (& grad)[6] )
{
  int qa, qb, qc;
  LagrangeBasis1::TensorProduct3D::multiIndex( q, qa, qb, qc );
 
  supportLoop( qa, qb, qc, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
                                                 int const nodeIndex,
                                                 real64 const (&invJ)[3][3],
                                                 real64 const (&var)[numNodes][3],
                                                 real64 (& grad)[6] )
  {
 
    real64 gradN[3] = {0, 0, 0};
    for( int i = 0; i < 3; ++i )
    {
      for( int j = 0; j < 3; ++j )
      {
        gradN[i] = gradN[i] + dNdXi[ j ] * invJ[j][i];
      }
    }
 
    grad[0] = grad[0] + gradN[0] * var[ nodeIndex ][0];
    grad[1] = grad[1] + gradN[1] * var[ nodeIndex ][1];
    grad[2] = grad[2] + gradN[2] * var[ nodeIndex ][2];
    grad[3] = grad[3] + gradN[2] * var[ nodeIndex ][1] + gradN[1] * var[ nodeIndex ][2];
    grad[4] = grad[4] + gradN[2] * var[ nodeIndex ][0] + gradN[0] * var[ nodeIndex ][2];
    grad[5] = grad[5] + gradN[1] * var[ nodeIndex ][0] + gradN[0] * var[ nodeIndex ][1];
  }, invJ, var, grad );
}
 
GEOS_HOST_DEVICE
inline
void H1_Hexahedron_Lagrange1_GaussLegendre2::plusGradNajAij( int const q,
                                                             real64 const (&invJ)[3][3],
                                                             real64 const (&var)[6],
                                                             real64 (& R)[numNodes][3] )
{
  int qa, qb, qc;
  LagrangeBasis1::TensorProduct3D::multiIndex( q, qa, qb, qc );
 
  supportLoop( qa, qb, qc,
               [] GEOS_HOST_DEVICE
                 ( real64 const (&dNdXi)[3],
                 int const nodeIndex,
                 real64 const (&invJ)[3][3],
                 real64 const (&var)[6],
                 real64 (& R)[numNodes][3] )
  {
 
    real64 gradN[3] = {0, 0, 0};
    for( int i = 0; i < 3; ++i )
    {
      for( int j = 0; j < 3; ++j )
      {
        gradN[i] = gradN[i] + dNdXi[ j ] * invJ[j][i];
      }
    }
    R[ nodeIndex ][ 0 ] = R[ nodeIndex ][ 0 ] - var[ 0 ] * gradN[ 0 ] - var[ 5 ] * gradN[ 1 ] - var[ 4 ] * gradN[ 2 ];
    R[ nodeIndex ][ 1 ] = R[ nodeIndex ][ 1 ] - var[ 5 ] * gradN[ 0 ] - var[ 1 ] * gradN[ 1 ] - var[ 3 ] * gradN[ 2 ];
    R[ nodeIndex ][ 2 ] = R[ nodeIndex ][ 2 ] - var[ 4 ] * gradN[ 0 ] - var[ 3 ] * gradN[ 1 ] - var[ 2 ] * gradN[ 2 ];
  }, invJ, var, R );
}
 
 
 
GEOS_HOST_DEVICE
inline
void H1_Hexahedron_Lagrange1_GaussLegendre2::gradient( int const q,
                                                       real64 const (&invJ)[3][3],
                                                       real64 const (&var)[numNodes][3],
                                                       real64 (& grad)[3][3] )
{
  int qa, qb, qc;
  LagrangeBasis1::TensorProduct3D::multiIndex( q, qa, qb, qc );
 
  supportLoop( qa, qb, qc, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
                                                 int const nodeIndex,
                                                 real64 const (&invJ)[3][3],
                                                 real64 const (&var)[numNodes][3],
                                                 real64 (& grad)[3][3] )
  {
    for( int i = 0; i < 3; ++i )
    {
      real64 gradN=0.0;;
      for( int j = 0; j < 3; ++j )
      {
        gradN = gradN + dNdXi[ j ] * invJ[j][i];
      }
      for( int k = 0; k < 3; ++k )
      {
        grad[k][i] = grad[k][i] + gradN * var[ nodeIndex ][k];
      }
    }
  }, invJ, var, grad );
}
 
 
#if __GNUC__
#pragma GCC diagnostic pop
#endif
 
#undef PARENT_GRADIENT_METHOD
}
}
 
#endif //FINITE_ELEMENT_SHAPE_KERNEL