GEOS
Qk_Hexahedron_Lagrange_GaussLobatto.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 TotalEnergies
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
20 #ifndef GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_Q1HEXAHEDRON_HPP_
21 #define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_Q1HEXAHEDRON_HPP_
22 
23 #include "FiniteElementBase.hpp"
24 #include "LagrangeBasis1.hpp"
25 #include "LagrangeBasis2.hpp"
26 #include "LagrangeBasis3GL.hpp"
27 #include "LagrangeBasis4GL.hpp"
28 #include "LagrangeBasis5GL.hpp"
29 #include <utility>
30 
31 
32 
33 namespace geos
34 {
35 namespace finiteElement
36 {
37 
43 template< typename GL_BASIS >
44 class Qk_Hexahedron_Lagrange_GaussLobatto_impl : public FiniteElementBase_impl< GL_BASIS::TensorProduct3D::numSupportPoints,
45  6,
46  GL_BASIS::TensorProduct3D::numSupportPoints >
47 {
48 public:
50  using Base = FiniteElementBase_impl< GL_BASIS::TensorProduct3D::numSupportPoints,
51  6,
52  GL_BASIS::TensorProduct3D::numSupportPoints >;
53 
55  using Base::numNodes;
56 
58  using Base::numQuadraturePoints;
59 
61  using Base::maxSupportPoints;
62 
64  constexpr static localIndex num1dNodes = GL_BASIS::numSupportPoints;
65 
67  constexpr static localIndex halfNodes = ( GL_BASIS::numSupportPoints - 1 )/ 2;
68 
70  constexpr static localIndex numNodesPerFace = GL_BASIS::TensorProduct2D::numSupportPoints;
71 
73  using StackVariables = typename Base::StackVariables;
74 
76  template< typename SubregionType >
77  using MeshData = typename Base::template MeshData< SubregionType >;
78 
80  Qk_Hexahedron_Lagrange_GaussLobatto_impl() = default;
81  ~Qk_Hexahedron_Lagrange_GaussLobatto_impl() = default;
82  Qk_Hexahedron_Lagrange_GaussLobatto_impl( Qk_Hexahedron_Lagrange_GaussLobatto_impl const & ) = default;
83  Qk_Hexahedron_Lagrange_GaussLobatto_impl & operator=( Qk_Hexahedron_Lagrange_GaussLobatto_impl const & ) = default;
84  Qk_Hexahedron_Lagrange_GaussLobatto_impl( Qk_Hexahedron_Lagrange_GaussLobatto_impl && ) = default;
85  Qk_Hexahedron_Lagrange_GaussLobatto_impl & operator=( Qk_Hexahedron_Lagrange_GaussLobatto_impl && ) = default;
87 
95  GEOS_HOST_DEVICE
96  GEOS_FORCE_INLINE
97  constexpr static localIndex linearIndex3DVal( const localIndex qa, localIndex const qb, localIndex const qc )
98  {
99  return qa + qb * num1dNodes + qc * numNodesPerFace;
100  }
101 
107  GEOS_HOST_DEVICE
108  GEOS_FORCE_INLINE
109  constexpr static localIndex meshIndexToLinearIndex3D( localIndex const k )
110  {
111  return linearIndex3DVal( ( num1dNodes - 1 ) * ( k % 2 ),
112  ( num1dNodes - 1 ) * ( ( k % 4 ) / 2 ),
113  ( num1dNodes - 1 ) * ( k / 4 ) );
114  }
115 
116 
123  GEOS_HOST_DEVICE
124  GEOS_FORCE_INLINE
125  constexpr static localIndex linearIndex2DVal( const localIndex qa, const localIndex qb )
126  {
127  return qa + qb * num1dNodes;
128  }
129 
135  GEOS_HOST_DEVICE
136  GEOS_FORCE_INLINE
137  constexpr static localIndex meshIndexToLinearIndex2D( localIndex const k )
138  {
139  return linearIndex2DVal( ( num1dNodes - 1 ) * ( k % 2 ),
140  ( num1dNodes - 1 ) * ( k / 2 ) );
141  }
142 
148  GEOS_HOST_DEVICE
149  GEOS_FORCE_INLINE
150  static localIndex getNumQuadraturePoints( StackVariables const & stack )
151  {
152  GEOS_UNUSED_VAR( stack );
153  return numQuadraturePoints;
154  }
155 
161  GEOS_HOST_DEVICE
162  GEOS_FORCE_INLINE
163  static localIndex getNumSupportPoints( StackVariables const & stack )
164  {
165  GEOS_UNUSED_VAR( stack );
166  return numNodes;
167  }
168 
169 
175  GEOS_HOST_DEVICE
176  GEOS_FORCE_INLINE
177  static void calcN( real64 const (&coords)[3],
178  real64 (& N)[numNodes] )
179  {
180  GL_BASIS::TensorProduct3D::value( coords, N );
181  }
182 
189  GEOS_HOST_DEVICE
190  GEOS_FORCE_INLINE
191  constexpr static real64 interpolationCoord( const int q, const int k )
192  {
193  const real64 alpha = ( GL_BASIS::parentSupportCoord( q ) + 1.0 ) / 2.0;
194  return k == 0 ? ( 1.0 - alpha ) : alpha;
195  }
196 
197 
204  GEOS_HOST_DEVICE
205  GEOS_FORCE_INLINE
206  constexpr static real64 basisGradientAt( const int q, const int p )
207  {
208  if( p <= halfNodes )
209  {
210  return GL_BASIS::gradientAt( q, p );
211  }
212  else
213  {
214  return -GL_BASIS::gradientAt( GL_BASIS::numSupportPoints - 1 - q, GL_BASIS::numSupportPoints - 1 - p );
215  }
216  }
217 
229  GEOS_HOST_DEVICE
230  GEOS_FORCE_INLINE
231  constexpr static real64 jacobianCoefficient1D( const int q, const int i, const int k, const int dir )
232  {
233  if( i == dir )
234  {
235  return k== 0 ? -1.0/2.0 : 1.0/2.0;
236  }
237  else
238  {
239  return interpolationCoord( q, k );
240  }
241  }
242 
250  GEOS_HOST_DEVICE
251  GEOS_FORCE_INLINE
252  static void calcN( localIndex const q,
253  real64 (& N)[numNodes] )
254  {
255  for( int a=0; a < numNodes; ++a )
256  {
257  N[ a ] = 0;
258  }
259  N[ q ] = 1.0;
260  }
261 
270  GEOS_HOST_DEVICE
271  GEOS_FORCE_INLINE
272  static void calcN( localIndex const q,
273  StackVariables const & stack,
274  real64 ( & N )[numNodes] )
275  {
276  GEOS_UNUSED_VAR( stack );
277  return calcN( q, N );
278  }
279 
287  GEOS_HOST_DEVICE
288  GEOS_FORCE_INLINE
289  static
290  real64 calcJacobian( localIndex const q,
291  real64 const (&X)[numNodes][3],
292  real64 ( &J )[3][3] );
293 
302  GEOS_HOST_DEVICE
303  GEOS_FORCE_INLINE
304  static
305  real64 calcJacobian( localIndex const q,
306  real64 const (&X)[numNodes][3],
307  StackVariables const &stack,
308  real64 ( &J )[3][3] );
309 
319  GEOS_HOST_DEVICE
320  GEOS_FORCE_INLINE
321  static real64 calcGradN( localIndex const q,
322  real64 const (&X)[numNodes][3],
323  real64 ( &gradN )[numNodes][3] );
333  GEOS_HOST_DEVICE
334  GEOS_FORCE_INLINE
335  static real64 calcGradN( real64 const (&coords)[3],
336  real64 const (&X)[numNodes][3],
337  real64 ( &gradN )[numNodes][3] );
338 
349  GEOS_HOST_DEVICE
350  GEOS_FORCE_INLINE
351  static real64 calcGradN( localIndex const q,
352  real64 const (&X)[numNodes][3],
353  StackVariables const & stack,
354  real64 ( &gradN )[numNodes][3] );
364  GEOS_HOST_DEVICE
365  GEOS_FORCE_INLINE
366  static real64 calcGradNWithCorners( localIndex const q,
367  real64 const (&X)[8][3],
368  real64 ( &gradN )[numNodes][3] );
378  GEOS_HOST_DEVICE
379  GEOS_FORCE_INLINE
380  static real64 calcGradNWithCorners( real64 const (&coords)[3],
381  real64 const (&X)[8][3],
382  real64 ( &gradN )[numNodes][3] );
383 
394  GEOS_HOST_DEVICE
395  GEOS_FORCE_INLINE
396  static real64 calcGradNWithCorners( localIndex const q,
397  real64 const (&X)[8][3],
398  StackVariables const & stack,
399  real64 ( &gradN )[numNodes][3] );
400 
408  GEOS_HOST_DEVICE
409  GEOS_FORCE_INLINE
410  static real64 transformedQuadratureWeight( localIndex const q,
411  real64 const (&X)[numNodes][3] );
412 
421  GEOS_HOST_DEVICE
422  GEOS_FORCE_INLINE
423  static void jacobianTransformation2d( int const qa,
424  int const qb,
425  real64 const (&X)[4][3],
426  real64 ( &J )[3][2] );
427 
428 
439  GEOS_HOST_DEVICE
440  GEOS_FORCE_INLINE
441  static real64 invJacobianTransformation( int const qa,
442  int const qb,
443  int const qc,
444  real64 const (&X)[8][3],
445  real64 ( & J )[3][3] )
446  {
447  jacobianTransformation( qa, qb, qc, X, J );
448  return LvArray::tensorOps::invert< 3 >( J );
449  }
450 
459  GEOS_HOST_DEVICE
460  GEOS_FORCE_INLINE
461  static real64 invJacobianTransformation( int const q,
462  real64 const (&X)[8][3],
463  real64 ( & J )[3][3] )
464  {
465  int qa, qb, qc;
466  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
467  return invJacobianTransformation( qa, qb, qc, X, J );
468  }
469 
470 
481  GEOS_HOST_DEVICE
482  GEOS_FORCE_INLINE
483  static void symmetricGradient( int const q,
484  real64 const (&invJ)[3][3],
485  real64 const (&var)[numNodes][3],
486  real64 ( &grad )[6] );
487 
488 
489 
505  GEOS_HOST_DEVICE
506  GEOS_FORCE_INLINE
507  static void gradient( int const q,
508  real64 const (&invJ)[3][3],
509  real64 const (&var)[numNodes][3],
510  real64 ( &grad )[3][3] );
511 
512 
528  GEOS_HOST_DEVICE
529  GEOS_FORCE_INLINE
530  static void plusGradNajAij( int const q,
531  real64 const (&invJ)[3][3],
532  real64 const (&var)[6],
533  real64 ( &R )[numNodes][3] );
534 
535 
536 
546  GEOS_HOST_DEVICE
547  GEOS_FORCE_INLINE
548  static void jacobianTransformation( int const qa,
549  int const qb,
550  int const qc,
551  real64 const (&X)[8][3],
552  real64 ( &J )[3][3] );
553 
561  GEOS_HOST_DEVICE
562  GEOS_FORCE_INLINE
563  static void jacobianTransformation( real64 const (&coords)[3],
564  real64 const (&X)[numNodes][3],
565  real64 ( &J )[3][3] );
566 
576  GEOS_HOST_DEVICE
577  GEOS_FORCE_INLINE
578  static void jacobianTransformationWithCorners( real64 const (&coords)[3],
579  real64 const (&X)[8][3],
580  real64 ( &J )[3][3] );
581 
591  GEOS_HOST_DEVICE
592  GEOS_FORCE_INLINE
593  static void trilinearInterp( real64 const alpha,
594  real64 const beta,
595  real64 const gamma,
596  real64 const (&X)[8][3],
597  real64 ( &coords )[3] );
598 
604  GEOS_HOST_DEVICE
605  GEOS_FORCE_INLINE
606  static void computeLocalCoords( real64 const (&Xmesh)[8][3],
607  real64 const (&X)[numNodes][3] );
608 
616  GEOS_HOST_DEVICE
617  GEOS_FORCE_INLINE
618  static real64 computeMassTerm( localIndex const q,
619  real64 const (&X)[8][3] );
620 
629  GEOS_HOST_DEVICE
630  GEOS_FORCE_INLINE
631  static real64 computeDampingTerm( localIndex const q,
632  real64 const (&X)[4][3] );
633 
644  GEOS_HOST_DEVICE
645  GEOS_FORCE_INLINE
646  static void computeBMatrix( int const qa,
647  int const qb,
648  int const qc,
649  real64 const (&X)[8][3],
650  real64 ( &J )[3][3],
651  real64 ( &B )[6] );
652 
661  template< typename FUNC >
662  GEOS_HOST_DEVICE
663  GEOS_FORCE_INLINE
664  static void computeStiffnessTerm( localIndex const q,
665  real64 const (&X)[8][3],
666  FUNC && func );
667 
678  GEOS_HOST_DEVICE
679  GEOS_FORCE_INLINE
680  static void computeBxyMatrix( int const qa,
681  int const qb,
682  int const qc,
683  real64 const (&X)[8][3],
684  real64 ( &J )[3][3],
685  real64 ( &B )[6] );
686 
695  template< typename FUNC >
696  GEOS_HOST_DEVICE
697  GEOS_FORCE_INLINE
698  static void computeStiffnessxyTerm( localIndex const q,
699  real64 const (&X)[8][3],
700  FUNC && func );
701 
712  GEOS_HOST_DEVICE
713  GEOS_FORCE_INLINE
714  static void computeBzMatrix( int const qa,
715  int const qb,
716  int const qc,
717  real64 const (&X)[8][3],
718  real64 ( &J )[3][3],
719  real64 ( &B )[6] );
720 
729  template< typename FUNC >
730  GEOS_HOST_DEVICE
731  GEOS_FORCE_INLINE
732  static void computeStiffnesszTerm( localIndex const q,
733  real64 const (&X)[8][3],
734  FUNC && func );
735 
745  template< typename FUNC >
746  GEOS_HOST_DEVICE
747  GEOS_FORCE_INLINE
748  static void
749  computeGradPhiBGradPhi( int const qa,
750  int const qb,
751  int const qc,
752  real64 const (&B)[6],
753  FUNC && func );
754 
763  template< typename FUNC >
764  GEOS_HOST_DEVICE
765  GEOS_FORCE_INLINE
766  static void computeFirstOrderStiffnessTermX( localIndex const q,
767  real64 const (&X)[8][3],
768  FUNC && func );
777  template< typename FUNC >
778  GEOS_HOST_DEVICE
779  GEOS_FORCE_INLINE
780  static void computeFirstOrderStiffnessTermY( localIndex const q,
781  real64 const (&X)[8][3],
782  FUNC && func );
791  template< typename FUNC >
792  GEOS_HOST_DEVICE
793  GEOS_FORCE_INLINE
794  static void computeFirstOrderStiffnessTermZ( localIndex const q,
795  real64 const (&X)[8][3],
796  FUNC && func );
806  template< typename FUNC >
807  GEOS_HOST_DEVICE
808  GEOS_FORCE_INLINE
809  static void computeFirstOrderStiffnessTerm( localIndex const q,
810  real64 const (&X)[8][3],
811  FUNC && stiffnessVal );
812 
813 
823  GEOS_HOST_DEVICE
824  GEOS_FORCE_INLINE
825  static void applyTransformationToParentGradients( int const q,
826  real64 const ( &invJ )[3][3],
827  real64 ( &gradN )[numNodes][3] );
828 
838  GEOS_HOST_DEVICE
839  GEOS_FORCE_INLINE
840  static void applyTransformationToParentGradients( real64 const (&coords)[3],
841  real64 const ( &invJ )[3][3],
842  real64 ( &gradN )[numNodes][3] );
843 
844 
845 private:
847  constexpr static real64 parentLength = GL_BASIS::parentSupportCoord( 1 ) - GL_BASIS::parentSupportCoord( 0 );
848 
850  constexpr static real64 parentVolume = parentLength*parentLength*parentLength;
860  template< typename FUNC, typename ... PARAMS >
861  GEOS_HOST_DEVICE
862  GEOS_FORCE_INLINE
863  static void supportLoop( real64 const (&coords)[3],
864  FUNC && func,
865  PARAMS &&... params );
875  template< typename FUNC, typename ... PARAMS >
876  GEOS_HOST_DEVICE
877  GEOS_FORCE_INLINE
878  static void supportLoop( localIndex const q,
879  FUNC && func,
880  PARAMS &&... params );
881 
882 };
883 
885 
886 
887 template< typename GL_BASIS >
888 template< typename FUNC, typename ... PARAMS >
889 GEOS_HOST_DEVICE
890 GEOS_FORCE_INLINE
891 void
892 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::supportLoop( real64 const (&coords)[3],
893  FUNC && func,
894  PARAMS &&... params )
895 {
896  for( int c=0; c<num1dNodes; ++c )
897  {
898  for( int b=0; b<num1dNodes; ++b )
899  {
900  for( int a=0; a<num1dNodes; ++a )
901  {
902  real64 const dNdXi[3] = { GL_BASIS::gradient( a, coords[0] )*
903  GL_BASIS::value( b, coords[1] )*
904  GL_BASIS::value( c, coords[2] ),
905  GL_BASIS::value( a, coords[0] )*
906  GL_BASIS::gradient( b, coords[1] )*
907  GL_BASIS::value( c, coords[2] ),
908  GL_BASIS::value( a, coords[0] )*
909  GL_BASIS::value( b, coords[1] )*
910  GL_BASIS::gradient( c, coords[2] )};
911 
912  localIndex const nodeIndex = GL_BASIS::TensorProduct3D::linearIndex( a, b, c );
913 
914  func( dNdXi, nodeIndex, std::forward< PARAMS >( params )... );
915  }
916  }
917  }
918 }
919 
920 template< typename GL_BASIS >
921 template< typename FUNC, typename ... PARAMS >
922 GEOS_HOST_DEVICE
923 GEOS_FORCE_INLINE
924 void
925 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::supportLoop( localIndex const q,
926  FUNC && func,
927  PARAMS &&... params )
928 {
929  int qa, qb, qc;
930  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
931  for( int c=0; c<num1dNodes; ++c )
932  {
933  for( int b=0; b<num1dNodes; ++b )
934  {
935  for( int a=0; a<num1dNodes; ++a )
936  {
937  real64 const dNdXi[3] = { (b == qb && c == qc ) ? basisGradientAt( a, qa ) : 0,
938  (a == qa && c == qc ) ? basisGradientAt( b, qb ) : 0,
939  (a == qa && b == qb ) ? basisGradientAt( c, qc ) : 0 };
940 
941  localIndex const nodeIndex = GL_BASIS::TensorProduct3D::linearIndex( a, b, c );
942 
943  func( dNdXi, nodeIndex, std::forward< PARAMS >( params )... );
944  }
945  }
946  }
947 }
948 
949 //*************************************************************************************************
950 
951 template< typename GL_BASIS >
952 GEOS_HOST_DEVICE
953 GEOS_FORCE_INLINE
954 real64
955 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::calcGradN( localIndex const q,
956  real64 const (&X)[numNodes][3],
957  real64 (& gradN)[numNodes][3] )
958 {
959  int qa, qb, qc;
960  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
961  real64 Xmesh[8][3] = {{0}};
962  for( int k = 0; k < 8; k++ )
963  {
964  const localIndex nodeIndex = meshIndexToLinearIndex3D( k );
965  for( int i = 0; i < 3; i++ )
966  {
967  Xmesh[ k ][ i ] = X[ nodeIndex ][ i ];
968  }
969  }
970  real64 J[3][3] = {{0}};
971 
972  jacobianTransformation( qa, qb, qc, Xmesh, J );
973 
974  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
975 
976  applyTransformationToParentGradients( q, J, gradN );
977 
978  return detJ * GL_BASIS::weight( qa ) * GL_BASIS::weight( qb ) * GL_BASIS::weight( qc );
979 }
980 //*************************************************************************************************
981 template< typename GL_BASIS >
982 GEOS_HOST_DEVICE
983 GEOS_FORCE_INLINE
984 real64
985 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::calcJacobian( localIndex const q,
986  real64 const (&X)[numNodes][3],
987  real64 (& J)[3][3] )
988 {
989  for( int i = 0; i < 3; ++i )
990  {
991  for( int j = 0; j < 3; ++j )
992  {
993  J[i][j] = 0;
994  }
995  }
996  int qa, qb, qc;
997  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
998  real64 Xmesh[8][3] = {{0}};
999  for( int k = 0; k < 8; k++ )
1000  {
1001  const localIndex nodeIndex = meshIndexToLinearIndex3D( k );
1002  for( int i = 0; i < 3; i++ )
1003  {
1004  Xmesh[k][i] = X[nodeIndex][i];
1005  }
1006  }
1007  jacobianTransformation( qa, qb, qc, Xmesh, J );
1008  return LvArray::tensorOps::determinant< 3 >( J );
1009 }
1010 
1011 template< typename GL_BASIS >
1012 GEOS_HOST_DEVICE
1013 GEOS_FORCE_INLINE
1014 real64
1015 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::calcJacobian( localIndex const q,
1016  real64 const (&X)[numNodes][3],
1017  StackVariables const & GEOS_UNUSED_PARAM( stack ),
1018  real64 (& J)[3][3] )
1019 {
1020  return calcJacobian( q, X, J );
1021 }
1022 
1023 template< typename GL_BASIS >
1024 GEOS_HOST_DEVICE
1025 GEOS_FORCE_INLINE
1026 real64
1027 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::calcGradN( real64 const (&coords)[3],
1028  real64 const (&X)[numNodes][3],
1029  real64 (& gradN)[numNodes][3] )
1030 {
1031  real64 J[3][3] = {{0}};
1032 
1033  jacobianTransformation( coords, X, J );
1034 
1035  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
1036 
1037  applyTransformationToParentGradients( coords, J, gradN );
1038 
1039  return detJ;
1040 }
1041 template< typename GL_BASIS >
1042 GEOS_HOST_DEVICE
1043 GEOS_FORCE_INLINE
1044 real64 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1045 calcGradN( localIndex const q,
1046  real64 const (&X)[numNodes][3],
1047  StackVariables const & GEOS_UNUSED_PARAM( stack ),
1048  real64 ( & gradN )[numNodes][3] )
1049 {
1050  return calcGradN( q, X, gradN );
1051 }
1052 
1053 template< typename GL_BASIS >
1054 GEOS_HOST_DEVICE
1055 GEOS_FORCE_INLINE
1056 real64
1057 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::calcGradNWithCorners( localIndex const q,
1058  real64 const (&X)[8][3],
1059  real64 (& gradN)[numNodes][3] )
1060 {
1061  int qa, qb, qc;
1062  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1063 
1064  real64 J[3][3] = {{0}};
1065 
1066  jacobianTransformation( qa, qb, qc, X, J );
1067 
1068  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
1069 
1070  applyTransformationToParentGradients( q, J, gradN );
1071 
1072  return detJ * GL_BASIS::weight( qa ) * GL_BASIS::weight( qb ) * GL_BASIS::weight( qc );
1073 }
1074 //*************************************************************************************************
1075 template< typename GL_BASIS >
1076 GEOS_HOST_DEVICE
1077 GEOS_FORCE_INLINE
1078 real64
1079 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::calcGradNWithCorners( real64 const (&coords)[3],
1080  real64 const (&X)[8][3],
1081  real64 (& gradN)[numNodes][3] )
1082 {
1083  real64 J[3][3] = {{0}};
1084 
1085  jacobianTransformationWithCorners( coords, X, J );
1086 
1087  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
1088 
1089  applyTransformationToParentGradients( coords, J, gradN );
1090 
1091  return detJ;
1092 }
1093 template< typename GL_BASIS >
1094 GEOS_HOST_DEVICE
1095 GEOS_FORCE_INLINE
1096 real64 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1097 calcGradNWithCorners( localIndex const q,
1098  real64 const (&X)[8][3],
1099  StackVariables const & GEOS_UNUSED_PARAM( stack ),
1100  real64 ( & gradN )[numNodes][3] )
1101 {
1102  return calcGradN( q, X, gradN );
1103 }
1104 //*************************************************************************************************
1105 #if __GNUC__
1106 #pragma GCC diagnostic push
1107 #pragma GCC diagnostic ignored "-Wshadow"
1108 #endif
1109 
1110 template< typename GL_BASIS >
1111 GEOS_HOST_DEVICE
1112 GEOS_FORCE_INLINE
1113 void
1114 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1115 jacobianTransformation( int const qa,
1116  int const qb,
1117  int const qc,
1118  real64 const (&X)[8][3],
1119  real64 ( & J )[3][3] )
1120 {
1121  for( int k = 0; k < 8; k++ )
1122  {
1123  const int ka = k % 2;
1124  const int kb = ( k % 4 ) / 2;
1125  const int kc = k / 4;
1126  for( int j = 0; j < 3; j++ )
1127  {
1128  real64 jacCoeff = jacobianCoefficient1D( qa, 0, ka, j ) *
1129  jacobianCoefficient1D( qb, 1, kb, j ) *
1130  jacobianCoefficient1D( qc, 2, kc, j );
1131  for( int i = 0; i < 3; i++ )
1132  {
1133  J[i][j] += jacCoeff * X[k][i];
1134  }
1135  }
1136  }
1137 }
1138 
1139 template< typename GL_BASIS >
1140 GEOS_HOST_DEVICE
1141 GEOS_FORCE_INLINE
1142 void
1143 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1144 jacobianTransformation( real64 const (&coords)[3],
1145  real64 const (&X)[numNodes][3],
1146  real64 ( & J )[3][3] )
1147 {
1148  supportLoop( coords, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1149  int const nodeIndex,
1150  real64 const (&X)[numNodes][3],
1151  real64 (& J)[3][3] )
1152  {
1153  real64 const * const GEOS_RESTRICT Xnode = X[nodeIndex];
1154  for( int i = 0; i < 3; ++i )
1155  {
1156  for( int j = 0; j < 3; ++j )
1157  {
1158  J[i][j] = J[i][j] + dNdXi[ j ] * Xnode[i];
1159  }
1160  }
1161  }, X, J );
1162 }
1163 
1164 template< typename GL_BASIS >
1165 GEOS_HOST_DEVICE
1166 GEOS_FORCE_INLINE
1167 void
1168 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1169 jacobianTransformationWithCorners( real64 const (&coords)[3],
1170  real64 const (&X)[8][3],
1171  real64 ( & J )[3][3] )
1172 {
1173  supportLoop( coords, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1174  int const nodeIndex,
1175  real64 const (&X)[8][3],
1176  real64 (& J)[3][3] )
1177  {
1178  int qa, qb, qc;
1179  GL_BASIS::TensorProduct3D::multiIndex( nodeIndex, qa, qb, qc );
1180  real64 Xnode[3];
1181  real64 alpha = ( GL_BASIS::parentSupportCoord( qa ) + 1.0 ) / 2.0;
1182  real64 beta = ( GL_BASIS::parentSupportCoord( qb ) + 1.0 ) / 2.0;
1183  real64 gamma = ( GL_BASIS::parentSupportCoord( qc ) + 1.0 ) / 2.0;
1184  trilinearInterp( alpha, beta, gamma, X, Xnode );
1185  for( int i = 0; i < 3; ++i )
1186  {
1187  for( int j = 0; j < 3; ++j )
1188  {
1189  J[i][j] = J[i][j] + dNdXi[ j ] * Xnode[i];
1190  }
1191  }
1192  }, X, J );
1193 }
1194 
1195 template< typename GL_BASIS >
1196 GEOS_HOST_DEVICE
1197 GEOS_FORCE_INLINE
1198 void
1199 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1200 trilinearInterp( real64 const alpha,
1201  real64 const beta,
1202  real64 const gamma,
1203  real64 const (&X)[8][3],
1204  real64 (& coords)[3] )
1205 {
1206  for( int i=0; i<3; i++ )
1207  {
1208  coords[i] = X[0][i]*( 1.0-alpha )*( 1.0-beta )*( 1.0-gamma )+
1209  X[1][i]* alpha *( 1.0-beta )*( 1.0-gamma )+
1210  X[2][i]*( 1.0-alpha )* beta *( 1.0-gamma )+
1211  X[3][i]* alpha * beta *( 1.0-gamma )+
1212  X[4][i]*( 1.0-alpha )*( 1.0-beta )* gamma+
1213  X[5][i]* alpha *( 1.0-beta )* gamma+
1214  X[6][i]*( 1.0-alpha )* beta * gamma+
1215  X[7][i]* alpha * beta * gamma;
1216  }
1217 }
1218 
1219 
1220 template< typename GL_BASIS >
1221 GEOS_HOST_DEVICE
1222 GEOS_FORCE_INLINE
1223 void
1224 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1225 computeLocalCoords( real64 const (&Xmesh)[8][3],
1226  real64 const (&X)[numNodes][3] )
1227 {
1228  int qa, qb, qc;
1229  for( int q=0; q<numNodes; q++ )
1230  {
1231  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1232  real64 alpha = ( GL_BASIS::parentSupportCoord( qa ) + 1.0 ) / 2.0;
1233  real64 beta = ( GL_BASIS::parentSupportCoord( qb ) + 1.0 ) / 2.0;
1234  real64 gamma = ( GL_BASIS::parentSupportCoord( qc ) + 1.0 ) / 2.0;
1235  trilinearInterp( alpha, beta, gamma, Xmesh, X[q] );
1236  }
1237 }
1238 
1239 template< typename GL_BASIS >
1240 GEOS_HOST_DEVICE
1241 GEOS_FORCE_INLINE
1242 void
1243 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1244 jacobianTransformation2d( int const qa,
1245  int const qb,
1246  real64 const (&X)[4][3],
1247  real64 ( & J )[3][2] )
1248 {
1249  for( int k = 0; k < 4; k++ )
1250  {
1251  int ka = k % 2;
1252  int kb = k / 2;
1253  for( int j = 0; j < 2; j++ )
1254  {
1255  real64 jacCoeff = jacobianCoefficient1D( qa, 0, ka, j ) *
1256  jacobianCoefficient1D( qb, 1, kb, j );
1257  for( int i = 0; i < 3; i++ )
1258  {
1259  J[i][j] += jacCoeff * X[k][i];
1260  }
1261  }
1262  }
1263 }
1264 
1265 template< typename GL_BASIS >
1266 GEOS_HOST_DEVICE
1267 GEOS_FORCE_INLINE
1268 real64
1269 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1270 computeMassTerm( localIndex const q,
1271  real64 const (&X)[8][3] )
1272 {
1273  int qa, qb, qc;
1274  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1275  const real64 w3D = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1276  real64 J[3][3] = {{0}};
1277  jacobianTransformation( qa, qb, qc, X, J );
1278  return LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( J ) )*w3D;
1279 }
1280 
1281 template< typename GL_BASIS >
1282 GEOS_HOST_DEVICE
1283 GEOS_FORCE_INLINE
1284 real64
1285 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1286 computeDampingTerm( localIndex const q,
1287  real64 const (&X)[4][3] )
1288 {
1289  int qa, qb;
1290  GL_BASIS::TensorProduct2D::multiIndex( q, qa, qb );
1291  const real64 w2D = GL_BASIS::weight( qa )*GL_BASIS::weight( qb );
1292  real64 B[3];
1293  real64 J[3][2] = {{0}};
1294  jacobianTransformation2d( qa, qb, X, J );
1295  // compute J^T.J, using Voigt notation for B
1296  B[0] = J[0][0]*J[0][0]+J[1][0]*J[1][0]+J[2][0]*J[2][0];
1297  B[1] = J[0][1]*J[0][1]+J[1][1]*J[1][1]+J[2][1]*J[2][1];
1298  B[2] = J[0][0]*J[0][1]+J[1][0]*J[1][1]+J[2][0]*J[2][1];
1299  return sqrt( LvArray::math::abs( LvArray::tensorOps::symDeterminant< 2 >( B ) ) )*w2D;
1300 }
1301 
1302 template< typename GL_BASIS >
1303 GEOS_HOST_DEVICE
1304 GEOS_FORCE_INLINE
1305 void
1306 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1307 computeBMatrix( int const qa,
1308  int const qb,
1309  int const qc,
1310  real64 const (&X)[8][3],
1311  real64 (& J)[3][3],
1312  real64 (& B)[6] )
1313 {
1314  jacobianTransformation( qa, qb, qc, X, J );
1315  real64 const detJ = LvArray::tensorOps::determinant< 3 >( J );
1316 
1317  // compute J^T.J/det(J), using Voigt notation for B
1318  B[0] = (J[0][0]*J[0][0]+J[1][0]*J[1][0]+J[2][0]*J[2][0])/detJ;
1319  B[1] = (J[0][1]*J[0][1]+J[1][1]*J[1][1]+J[2][1]*J[2][1])/detJ;
1320  B[2] = (J[0][2]*J[0][2]+J[1][2]*J[1][2]+J[2][2]*J[2][2])/detJ;
1321  B[3] = (J[0][1]*J[0][2]+J[1][1]*J[1][2]+J[2][1]*J[2][2])/detJ;
1322  B[4] = (J[0][0]*J[0][2]+J[1][0]*J[1][2]+J[2][0]*J[2][2])/detJ;
1323  B[5] = (J[0][0]*J[0][1]+J[1][0]*J[1][1]+J[2][0]*J[2][1])/detJ;
1324 
1325  // compute detJ*J^{-1}J^{-T}
1326  LvArray::tensorOps::symInvert< 3 >( B );
1327 }
1328 
1329 template< typename GL_BASIS >
1330 GEOS_HOST_DEVICE
1331 GEOS_FORCE_INLINE
1332 void
1333 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1334 computeBzMatrix( int const qa,
1335  int const qb,
1336  int const qc,
1337  real64 const (&X)[8][3],
1338  real64 (& J)[3][3],
1339  real64 (& B)[6] )
1340 {
1341  jacobianTransformation( qa, qb, qc, X, J );
1342  real64 const detJ = LvArray::tensorOps::determinant< 3 >( J );
1343 
1344  real64 Jinv[3][3] = {{0}};
1345  LvArray::tensorOps::invert< 3 >( Jinv, J );
1346 
1347  // compute det(J)*J^{-1}Az*J^{-T}, using Voigt notation for B
1348  B[0] = detJ*(Jinv[0][2]*Jinv[0][2]);
1349  B[1] = detJ*(Jinv[1][2]*Jinv[1][2]);
1350  B[2] = detJ*(Jinv[2][2]*Jinv[2][2]);
1351  B[3] = detJ*(Jinv[1][2]*Jinv[2][2]);
1352  B[4] = detJ*(Jinv[0][2]*Jinv[2][2]);
1353  B[5] = detJ*(Jinv[0][2]*Jinv[1][2]);
1354 }
1355 
1356 template< typename GL_BASIS >
1357 GEOS_HOST_DEVICE
1358 GEOS_FORCE_INLINE
1359 void
1360 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1361 computeBxyMatrix( int const qa,
1362  int const qb,
1363  int const qc,
1364  real64 const (&X)[8][3],
1365  real64 (& J)[3][3],
1366  real64 (& B)[6] )
1367 {
1368  jacobianTransformation( qa, qb, qc, X, J );
1369  real64 const detJ = LvArray::tensorOps::determinant< 3 >( J );
1370 
1371  real64 Jinv[3][3] = {{0}};
1372  LvArray::tensorOps::invert< 3 >( Jinv, J );
1373 
1374  // compute det(J)*J^{-1}Axy*J^{-T}, using Voigt notation for B
1375  B[0] = detJ*(Jinv[0][0]*Jinv[0][0] + Jinv[0][1]*Jinv[0][1]);
1376  B[1] = detJ*(Jinv[1][1]*Jinv[1][1] + Jinv[1][0]*Jinv[1][0]);
1377  B[2] = detJ*(Jinv[2][0]*Jinv[2][0] + Jinv[2][1]*Jinv[2][1]);
1378  B[3] = detJ*(Jinv[1][0]*Jinv[2][0] + Jinv[1][1]*Jinv[2][1]);
1379  B[4] = detJ*(Jinv[0][0]*Jinv[2][0] + Jinv[0][1]*Jinv[2][1]);
1380  B[5] = detJ*(Jinv[0][0]*Jinv[1][0] + Jinv[0][1]*Jinv[1][1]);
1381 }
1382 
1383 template< typename GL_BASIS >
1384 template< typename FUNC >
1385 GEOS_HOST_DEVICE
1386 GEOS_FORCE_INLINE
1387 void
1388 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1389 computeGradPhiBGradPhi( int const qa,
1390  int const qb,
1391  int const qc,
1392  real64 const (&B)[6],
1393  FUNC && func )
1394 {
1395  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1396  for( int i=0; i<num1dNodes; i++ )
1397  {
1398  const int ibc = GL_BASIS::TensorProduct3D::linearIndex( i, qb, qc );
1399  const int aic = GL_BASIS::TensorProduct3D::linearIndex( qa, i, qc );
1400  const int abi = GL_BASIS::TensorProduct3D::linearIndex( qa, qb, i );
1401  const real64 gia = basisGradientAt( i, qa );
1402  const real64 gib = basisGradientAt( i, qb );
1403  const real64 gic = basisGradientAt( i, qc );
1404  for( int j=0; j<num1dNodes; j++ )
1405  {
1406  const int jbc = GL_BASIS::TensorProduct3D::linearIndex( j, qb, qc );
1407  const int ajc = GL_BASIS::TensorProduct3D::linearIndex( qa, j, qc );
1408  const int abj = GL_BASIS::TensorProduct3D::linearIndex( qa, qb, j );
1409  const real64 gja = basisGradientAt( j, qa );
1410  const real64 gjb = basisGradientAt( j, qb );
1411  const real64 gjc = basisGradientAt( j, qc );
1412  // diagonal terms
1413  const real64 w0 = w * gia * gja;
1414  func( ibc, jbc, w0 * B[0] );
1415  const real64 w1 = w * gib * gjb;
1416  func( aic, ajc, w1 * B[1] );
1417  const real64 w2 = w * gic * gjc;
1418  func( abi, abj, w2 * B[2] );
1419  // off-diagonal terms
1420  const real64 w3 = w * gib * gjc;
1421  func( aic, abj, w3 * B[3] );
1422  func( abj, aic, w3 * B[3] );
1423  const real64 w4 = w * gia * gjc;
1424  func( ibc, abj, w4 * B[4] );
1425  func( abj, ibc, w4 * B[4] );
1426  const real64 w5 = w * gia * gjb;
1427  func( ibc, ajc, w5 * B[5] );
1428  func( ajc, ibc, w5 * B[5] );
1429  }
1430  }
1431 }
1432 
1433 template< typename GL_BASIS >
1434 template< typename FUNC >
1435 GEOS_HOST_DEVICE
1436 GEOS_FORCE_INLINE
1437 void
1438 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1439 computeStiffnessxyTerm( localIndex const q,
1440  real64 const (&X)[8][3],
1441  FUNC && func )
1442 {
1443  int qa, qb, qc;
1444  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1445  real64 B[6] = {0};
1446  real64 J[3][3] = {{0}};
1447  computeBxyMatrix( qa, qb, qc, X, J, B ); // The only change!
1448  computeGradPhiBGradPhi( qa, qb, qc, B, func );
1449 }
1450 
1451 template< typename GL_BASIS >
1452 template< typename FUNC >
1453 GEOS_HOST_DEVICE
1454 GEOS_FORCE_INLINE
1455 void
1456 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1457 computeStiffnesszTerm( localIndex const q,
1458  real64 const (&X)[8][3],
1459  FUNC && func )
1460 {
1461  int qa, qb, qc;
1462  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1463  real64 B[6] = {0};
1464  real64 J[3][3] = {{0}};
1465  computeBzMatrix( qa, qb, qc, X, J, B ); // The only change!
1466  computeGradPhiBGradPhi( qa, qb, qc, B, func );
1467 }
1468 
1469 template< typename GL_BASIS >
1470 template< typename FUNC >
1471 GEOS_HOST_DEVICE
1472 GEOS_FORCE_INLINE
1473 void
1474 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1475 computeStiffnessTerm( localIndex const q,
1476  real64 const (&X)[8][3],
1477  FUNC && func )
1478 {
1479  int qa, qb, qc;
1480  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1481  real64 B[6] = {0};
1482  real64 J[3][3] = {{0}};
1483  computeBMatrix( qa, qb, qc, X, J, B );
1484  computeGradPhiBGradPhi( qa, qb, qc, B, func );
1485 }
1486 
1487 template< typename GL_BASIS >
1488 template< typename FUNC >
1489 GEOS_HOST_DEVICE
1490 GEOS_FORCE_INLINE
1491 void
1492 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1493 computeFirstOrderStiffnessTerm( localIndex const q,
1494  real64 const (&X)[8][3],
1495  FUNC && func )
1496 {
1497  int qa, qb, qc;
1498  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1499  real64 J[3][3] = {{0}};
1500  jacobianTransformation( qa, qb, qc, X, J );
1501  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
1502  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1503  for( int i=0; i<num1dNodes; i++ )
1504  {
1505  const int ibc = GL_BASIS::TensorProduct3D::linearIndex( i, qb, qc );
1506  const int aic = GL_BASIS::TensorProduct3D::linearIndex( qa, i, qc );
1507  const int abi = GL_BASIS::TensorProduct3D::linearIndex( qa, qb, i );
1508  const real64 gia = basisGradientAt( i, qa );
1509  const real64 gib = basisGradientAt( i, qb );
1510  const real64 gic = basisGradientAt( i, qc );
1511  for( int j=0; j<num1dNodes; j++ )
1512  {
1513  const int jbc = GL_BASIS::TensorProduct3D::linearIndex( j, qb, qc );
1514  const int ajc = GL_BASIS::TensorProduct3D::linearIndex( qa, j, qc );
1515  const int abj = GL_BASIS::TensorProduct3D::linearIndex( qa, qb, j );
1516  const real64 gja = basisGradientAt( j, qa );
1517  const real64 gjb = basisGradientAt( j, qb );
1518  const real64 gjc = basisGradientAt( j, qc );
1519  // diagonal terms
1520  const real64 w00 = w * gia * gja;
1521  func( ibc, jbc, w00 * detJ, J, 0, 0 );
1522  const real64 w11 = w * gib * gjb;
1523  func( aic, ajc, w11 * detJ, J, 1, 1 );
1524  const real64 w22 = w * gic * gjc;
1525  func( abi, abj, w22 * detJ, J, 2, 2 );
1526  // off-diagonal terms
1527  const real64 w12 = w * gib * gjc;
1528  func( aic, abj, w12 * detJ, J, 1, 2 );
1529  func( abj, aic, w12 * detJ, J, 2, 1 );
1530  const real64 w02 = w * gia * gjc;
1531  func( ibc, abj, w02 * detJ, J, 0, 2 );
1532  func( abj, ibc, w02 * detJ, J, 2, 0 );
1533  const real64 w01 = w * gia * gjb;
1534  func( ibc, ajc, w01 * detJ, J, 0, 1 );
1535  func( ajc, ibc, w01 * detJ, J, 1, 0 );
1536  }
1537  }
1538 }
1539 
1540 template< typename GL_BASIS >
1541 template< typename FUNC >
1542 GEOS_HOST_DEVICE
1543 GEOS_FORCE_INLINE
1544 void
1545 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1546 computeFirstOrderStiffnessTermX( localIndex const q,
1547  real64 const (&X)[8][3],
1548  FUNC && func )
1549 {
1550  int qa, qb, qc;
1551  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1552  real64 J[3][3] = {{0}};
1553  jacobianTransformation( qa, qb, qc, X, J );
1554  const real64 detJ = LvArray::tensorOps::invert< 3 >( J );
1555  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1556 
1557  for( int i1 = 0; i1 < num1dNodes; ++i1 )
1558  {
1559  auto val = w * basisGradientAt( i1, qa );
1560  func( GL_BASIS::TensorProduct3D::linearIndex( i1, qb, qc ), q, detJ*J[0][0]*val, detJ*J[0][1]*val, detJ*J[0][2]*val );
1561  }
1562 
1563 }
1564 
1565 template< typename GL_BASIS >
1566 template< typename FUNC >
1567 GEOS_HOST_DEVICE
1568 GEOS_FORCE_INLINE
1569 void
1570 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1571 computeFirstOrderStiffnessTermY( localIndex const q,
1572  real64 const (&X)[8][3],
1573  FUNC && func )
1574 {
1575  int qa, qb, qc;
1576  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1577  real64 J[3][3] = {{0}};
1578  jacobianTransformation( qa, qb, qc, X, J );
1579  const real64 detJ = LvArray::tensorOps::invert< 3 >( J );
1580  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1581 
1582  for( int i2 = 0; i2 < num1dNodes; ++i2 )
1583  {
1584  auto val = w * basisGradientAt( i2, qb );
1585  func( GL_BASIS::TensorProduct3D::linearIndex( qa, i2, qc ), q, detJ*J[1][0]*val, detJ*J[1][1]*val, detJ*J[1][2]*val );
1586  }
1587 }
1588 
1589 template< typename GL_BASIS >
1590 template< typename FUNC >
1591 GEOS_HOST_DEVICE
1592 GEOS_FORCE_INLINE
1593 void
1594 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1595 computeFirstOrderStiffnessTermZ( localIndex const q,
1596  real64 const (&X)[8][3],
1597  FUNC && func )
1598 {
1599  int qa, qb, qc;
1600  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1601  real64 J[3][3] = {{0}};
1602  jacobianTransformation( qa, qb, qc, X, J );
1603  const real64 detJ = LvArray::tensorOps::invert< 3 >( J );
1604  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1605 
1606  for( int i3 = 0; i3 < num1dNodes; ++i3 )
1607  {
1608  auto val = w * basisGradientAt( i3, qc );
1609  func( GL_BASIS::TensorProduct3D::linearIndex( qa, qb, i3 ), q, detJ*J[2][0]*val, detJ*J[2][1]*val, detJ*J[2][2]*val );
1610  }
1611 }
1612 
1613 //*************************************************************************************************
1614 template< typename GL_BASIS >
1615 GEOS_HOST_DEVICE
1616 GEOS_FORCE_INLINE
1617 void
1618 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1619 applyTransformationToParentGradients( int const q,
1620  real64 const ( &invJ )[3][3],
1621  real64 (& gradN)[numNodes][3] )
1622 {
1623  supportLoop( q, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1624  int const nodeIndex,
1625  real64 const (&invJ)[3][3],
1626  real64 (& gradN)[numNodes][3] )
1627  {
1628  // smaller register footprint by manually unrolling the for loops.
1629  gradN[nodeIndex][0] = dNdXi[0] * invJ[0][0] + dNdXi[1] * invJ[1][0] + dNdXi[2] * invJ[2][0];
1630  gradN[nodeIndex][1] = dNdXi[0] * invJ[0][1] + dNdXi[1] * invJ[1][1] + dNdXi[2] * invJ[2][1];
1631  gradN[nodeIndex][2] = dNdXi[0] * invJ[0][2] + dNdXi[1] * invJ[1][2] + dNdXi[2] * invJ[2][2];
1632 
1633 
1634  }, invJ, gradN );
1635 }
1636 
1637 //*************************************************************************************************
1638 template< typename GL_BASIS >
1639 GEOS_HOST_DEVICE
1640 GEOS_FORCE_INLINE
1641 void
1642 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1643 applyTransformationToParentGradients( real64 const (&coords)[3],
1644  real64 const ( &invJ )[3][3],
1645  real64 (& gradN)[numNodes][3] )
1646 {
1647  supportLoop( coords, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1648  int const nodeIndex,
1649  real64 const (&invJ)[3][3],
1650  real64 (& gradN)[numNodes][3] )
1651  {
1652  gradN[nodeIndex][0] = dNdXi[0] * invJ[0][0] + dNdXi[1] * invJ[1][0] + dNdXi[2] * invJ[2][0];
1653  gradN[nodeIndex][1] = dNdXi[0] * invJ[0][1] + dNdXi[1] * invJ[1][1] + dNdXi[2] * invJ[2][1];
1654  gradN[nodeIndex][2] = dNdXi[0] * invJ[0][2] + dNdXi[1] * invJ[1][2] + dNdXi[2] * invJ[2][2];
1655  }, invJ, gradN );
1656 }
1657 
1658 template< typename GL_BASIS >
1659 GEOS_HOST_DEVICE
1660 GEOS_FORCE_INLINE
1661 real64
1662 Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1663 transformedQuadratureWeight( localIndex const q,
1664  real64 const (&X)[numNodes][3] )
1665 {
1666  int qa, qb, qc;
1667  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1668  real64 J[3][3] = {{0}};
1669 
1670  jacobianTransformation( qa, qb, qc, X, J );
1671 
1672  return LvArray::tensorOps::determinant< 3 >( J );
1673 }
1674 
1675 
1676 
1677 template< typename GL_BASIS >
1678 GEOS_HOST_DEVICE
1679 GEOS_FORCE_INLINE
1680 void Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1681 symmetricGradient( int const q,
1682  real64 const (&invJ)[3][3],
1683  real64 const (&var)[numNodes][3],
1684  real64 (& grad)[6] )
1685 {
1686  supportLoop( q, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1687  int const nodeIndex,
1688  real64 const (&invJ)[3][3],
1689  real64 const (&var)[numNodes][3],
1690  real64 (& grad)[6] )
1691  {
1692 
1693  real64 gradN[3] = {0, 0, 0};
1694  for( int i = 0; i < 3; ++i )
1695  {
1696  for( int j = 0; j < 3; ++j )
1697  {
1698  gradN[i] = gradN[i] + dNdXi[ j ] * invJ[j][i];
1699  }
1700  }
1701 
1702  grad[0] = grad[0] + gradN[0] * var[ nodeIndex ][0];
1703  grad[1] = grad[1] + gradN[1] * var[ nodeIndex ][1];
1704  grad[2] = grad[2] + gradN[2] * var[ nodeIndex ][2];
1705  grad[3] = grad[3] + gradN[2] * var[ nodeIndex ][1] + gradN[1] * var[ nodeIndex ][2];
1706  grad[4] = grad[4] + gradN[2] * var[ nodeIndex ][0] + gradN[0] * var[ nodeIndex ][2];
1707  grad[5] = grad[5] + gradN[1] * var[ nodeIndex ][0] + gradN[0] * var[ nodeIndex ][1];
1708  }, invJ, var, grad );
1709 }
1710 
1711 template< typename GL_BASIS >
1712 GEOS_HOST_DEVICE
1713 GEOS_FORCE_INLINE
1714 void Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1715 plusGradNajAij( int const q,
1716  real64 const (&invJ)[3][3],
1717  real64 const (&var)[6],
1718  real64 (& R)[numNodes][3] )
1719 {
1720  supportLoop( q,
1721  [] GEOS_HOST_DEVICE
1722  ( real64 const (&dNdXi)[3],
1723  int const nodeIndex,
1724  real64 const (&invJ)[3][3],
1725  real64 const (&var)[6],
1726  real64 (& R)[numNodes][3] )
1727  {
1728 
1729  real64 gradN[3] = {0, 0, 0};
1730  for( int i = 0; i < 3; ++i )
1731  {
1732  for( int j = 0; j < 3; ++j )
1733  {
1734  gradN[i] = gradN[i] + dNdXi[ j ] * invJ[j][i];
1735  }
1736  }
1737  R[ nodeIndex ][ 0 ] = R[ nodeIndex ][ 0 ] - var[ 0 ] * gradN[ 0 ] - var[ 5 ] * gradN[ 1 ] - var[ 4 ] * gradN[ 2 ];
1738  R[ nodeIndex ][ 1 ] = R[ nodeIndex ][ 1 ] - var[ 5 ] * gradN[ 0 ] - var[ 1 ] * gradN[ 1 ] - var[ 3 ] * gradN[ 2 ];
1739  R[ nodeIndex ][ 2 ] = R[ nodeIndex ][ 2 ] - var[ 4 ] * gradN[ 0 ] - var[ 3 ] * gradN[ 1 ] - var[ 2 ] * gradN[ 2 ];
1740  }, invJ, var, R );
1741 }
1742 
1743 
1744 
1745 template< typename GL_BASIS >
1746 GEOS_HOST_DEVICE
1747 GEOS_FORCE_INLINE
1748 void Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >::
1749 gradient( int const q,
1750  real64 const (&invJ)[3][3],
1751  real64 const (&var)[numNodes][3],
1752  real64 (& grad)[3][3] )
1753 {
1754  supportLoop( q, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1755  int const nodeIndex,
1756  real64 const (&invJ)[3][3],
1757  real64 const (&var)[numNodes][3],
1758  real64 (& grad)[3][3] )
1759  {
1760  for( int i = 0; i < 3; ++i )
1761  {
1762  real64 gradN=0.0;;
1763  for( int j = 0; j < 3; ++j )
1764  {
1765  gradN = gradN + dNdXi[ j ] * invJ[j][i];
1766  }
1767  for( int k = 0; k < 3; ++k )
1768  {
1769  grad[k][i] = grad[k][i] + gradN * var[ nodeIndex ][k];
1770  }
1771  }
1772  }, invJ, var, grad );
1773 }
1774 
1775 
1776 template< typename GL_BASIS >
1777 class Qk_Hexahedron_Lagrange_GaussLobatto final : public Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >,
1778  public FiniteElementBase
1779 {
1780 public:
1781 
1783  using BASIS = GL_BASIS;
1784 
1785  using ImplType = Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS >;
1786 
1787  Qk_Hexahedron_Lagrange_GaussLobatto(): FiniteElementBase( ImplType::numNodes,
1788  ImplType::maxSupportPoints,
1789  ImplType::numQuadraturePoints )
1790  { }
1791 
1796  Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS > * getImpl()
1797  {
1798  return static_cast< Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS > * >(this);
1799  }
1800 
1805  const Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS > * getImpl() const
1806  {
1807  return static_cast< const Qk_Hexahedron_Lagrange_GaussLobatto_impl< GL_BASIS > * >(this);
1808  }
1809 
1810  virtual ~Qk_Hexahedron_Lagrange_GaussLobatto() override final = default;
1811 };
1812 
1813 
1840 using Q1_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis1 >;
1883 using Q2_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis2 >;
1927 using Q3_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis3GL >;
1965 using Q4_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis4GL >;
2017 using Q5_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis5GL >;
2018 
2019 using Q1_Hexahedron_Lagrange_GaussLobatto_impl = Qk_Hexahedron_Lagrange_GaussLobatto_impl< LagrangeBasis1 >;
2020 using Q2_Hexahedron_Lagrange_GaussLobatto_impl = Qk_Hexahedron_Lagrange_GaussLobatto_impl< LagrangeBasis2 >;
2021 using Q3_Hexahedron_Lagrange_GaussLobatto_impl = Qk_Hexahedron_Lagrange_GaussLobatto_impl< LagrangeBasis3GL >;
2022 using Q4_Hexahedron_Lagrange_GaussLobatto_impl = Qk_Hexahedron_Lagrange_GaussLobatto_impl< LagrangeBasis4GL >;
2023 using Q5_Hexahedron_Lagrange_GaussLobatto_impl = Qk_Hexahedron_Lagrange_GaussLobatto_impl< LagrangeBasis5GL >;
2024 
2026 
2027 #if __GNUC__
2028 #pragma GCC diagnostic pop
2029 #endif
2030 #undef PARENT_GRADIENT_METHOD
2031 }
2032 }
2033 
2034 #endif //FINITE_ELEMENT_SHAPE_KERNEL
GEOS_HOST_DEVICE
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
GEOS_UNUSED_VAR
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
GEOS_FORCE_INLINE
#define GEOS_FORCE_INLINE
Marks a function or lambda for inlining.
Definition: GeosxMacros.hpp:51
GEOS_RESTRICT
#define GEOS_RESTRICT
preprocessor variable for the C99 restrict keyword for use with pointers
Definition: GeosxMacros.hpp:101
GEOS_UNUSED_PARAM
#define GEOS_UNUSED_PARAM(X)
Mark an unused argument and silence compiler warnings.
Definition: GeosxMacros.hpp:72
geos::finiteElement::FiniteElementBase_impl
Device-compatible (non virtual) Base class for all finite element formulations.
Definition: FiniteElementBase.hpp:42
geos::finiteElement::FiniteElementBase_impl< GL_BASIS::TensorProduct3D::numSupportPoints, 6, GL_BASIS::TensorProduct3D::numSupportPoints >::numQuadraturePoints
constexpr static localIndex numQuadraturePoints
The number of quadrature points per element.
Definition: FiniteElementBase.hpp:58
geos::finiteElement::FiniteElementBase_impl< GL_BASIS::TensorProduct3D::numSupportPoints, 6, GL_BASIS::TensorProduct3D::numSupportPoints >::maxSupportPoints
constexpr static localIndex maxSupportPoints
The maximum number of support points per element.
Definition: FiniteElementBase.hpp:52
geos::finiteElement::FiniteElementBase_impl< GL_BASIS::TensorProduct3D::numSupportPoints, 6, GL_BASIS::TensorProduct3D::numSupportPoints >::numNodes
constexpr static localIndex numNodes
The number of nodes per element.
Definition: FiniteElementBase.hpp:46
geos::finiteElement::FiniteElementBase_impl< GL_BASIS::TensorProduct3D::numSupportPoints, 6, GL_BASIS::TensorProduct3D::numSupportPoints >::operator=
GEOS_HOST_DEVICE FiniteElementBase_impl & operator=(FiniteElementBase_impl const &)=default
Default copy assignment operator.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:47
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeDampingTerm
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 computeDampingTerm(localIndex const q, real64 const (&X)[4][3])
computes the non-zero contributions of the d.o.f. indexd by q to the damping matrix M,...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeBMatrix
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeBMatrix(int const qa, int const qb, int const qc, real64 const (&X)[8][3], real64(&J)[3][3], real64(&B)[6])
computes the matrix B, defined as J^{-T}J^{-1}/det(J), where J is the Jacobian matrix,...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::jacobianCoefficient1D
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 jacobianCoefficient1D(const int q, const int i, const int k, const int dir)
Compute the 1D factor of the coefficient of the jacobian on the q-th quadrature point,...
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:231
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::linearIndex2DVal
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE localIndex linearIndex2DVal(const localIndex qa, const localIndex qb)
The linear index associated to the given one-dimensional indices in the two directions.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:125
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeBzMatrix
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeBzMatrix(int const qa, int const qb, int const qc, real64 const (&X)[8][3], real64(&J)[3][3], real64(&B)[6])
computes the matrix B in the case of quasi-stiffness (e.g. for pseudo-acoustic case),...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeMassTerm
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 computeMassTerm(localIndex const q, real64 const (&X)[8][3])
computes the non-zero contributions of the d.o.f. indexd by q to the mass matrix M,...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcN(localIndex const q, real64(&N)[numNodes])
Calculate shape functions values for each support point at a quadrature point.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:252
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::transformedQuadratureWeight
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 transformedQuadratureWeight(localIndex const q, real64 const (&X)[numNodes][3])
Calculate the integration weights for a quadrature point.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeLocalCoords
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeLocalCoords(real64 const (&Xmesh)[8][3], real64 const (&X)[numNodes][3])
computes the real-world coordinates of the support nodes
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeStiffnessxyTerm
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeStiffnessxyTerm(localIndex const q, real64 const (&X)[8][3], FUNC &&func)
computes the non-zero contributions of the d.o.f. indexed by q to the partial-stiffness matrix R,...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcGradNWithCorners
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcGradNWithCorners(localIndex const q, real64 const (&X)[8][3], real64(&gradN)[numNodes][3])
Calculate the shape functions derivatives wrt the physical coordinates.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::meshIndexToLinearIndex2D
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE localIndex meshIndexToLinearIndex2D(localIndex const k)
Converts from the index of the point in the mesh and the linear 2D index of the corresponding dof.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:137
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::interpolationCoord
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 interpolationCoord(const int q, const int k)
Compute the interpolation coefficients of the q-th quadrature point in a given direction.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:191
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeStiffnessTerm
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeStiffnessTerm(localIndex const q, real64 const (&X)[8][3], FUNC &&func)
computes the non-zero contributions of the d.o.f. indexed by q to the stiffness matrix R,...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcN(real64 const (&coords)[3], real64(&N)[numNodes])
Calculate shape functions values at a single point.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:177
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::halfNodes
constexpr static localIndex halfNodes
Half the number of support points, rounded down. Precomputed for efficiency.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:67
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::getNumQuadraturePoints
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE localIndex getNumQuadraturePoints(StackVariables const &stack)
Get the number of quadrature points.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:150
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcGradN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcGradN(localIndex const q, real64 const (&X)[numNodes][3], real64(&gradN)[numNodes][3])
Calculate the shape functions derivatives wrt the physical coordinates.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeGradPhiBGradPhi
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeGradPhiBGradPhi(int const qa, int const qb, int const qc, real64 const (&B)[6], FUNC &&func)
Computes the "Grad(Phi)*B*Grad(Phi)" coefficient of the stiffness term. The matrix B must be provided...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcGradNWithCorners
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcGradNWithCorners(localIndex const q, real64 const (&X)[8][3], StackVariables const &stack, real64(&gradN)[numNodes][3])
Calculate the shape functions derivatives wrt the physical coordinates.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcJacobian
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcJacobian(localIndex const q, real64 const (&X)[numNodes][3], StackVariables const &stack, real64(&J)[3][3])
Calculate the Jacobian matrix at a quadrature point.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeFirstOrderStiffnessTermZ
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeFirstOrderStiffnessTermZ(localIndex const q, real64 const (&X)[8][3], FUNC &&func)
computes the non-zero contributions of the d.o.f. indexd by q to the z-part of the first order stiffn...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::meshIndexToLinearIndex3D
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE localIndex meshIndexToLinearIndex3D(localIndex const k)
Converts from the index of the point in the mesh and the linear 3D index of the corresponding dof.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:109
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::basisGradientAt
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 basisGradientAt(const int q, const int p)
Compute the 1st derivative of the q-th 1D basis function at quadrature point p.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:206
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcJacobian
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcJacobian(localIndex const q, real64 const (&X)[numNodes][3], real64(&J)[3][3])
Calculate the Jacobian matrix at a quadrature point.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::trilinearInterp
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void trilinearInterp(real64 const alpha, real64 const beta, real64 const gamma, real64 const (&X)[8][3], real64(&coords)[3])
performs a trilinear interpolation to determine the real-world coordinates of a vertex
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::jacobianTransformation2d
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void jacobianTransformation2d(int const qa, int const qb, real64 const (&X)[4][3], real64(&J)[3][2])
Calculates the isoparametric "Jacobian" transformation matrix/mapping from the parent space to the ph...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::numNodesPerFace
constexpr static localIndex numNodesPerFace
The number of nodes/support points per face.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:70
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeFirstOrderStiffnessTermX
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeFirstOrderStiffnessTermX(localIndex const q, real64 const (&X)[8][3], FUNC &&func)
computes the non-zero contributions of the d.o.f. indexd by q to the x-part of the first order stiffn...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeFirstOrderStiffnessTerm
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeFirstOrderStiffnessTerm(localIndex const q, real64 const (&X)[8][3], FUNC &&stiffnessVal)
computes the non-zero contributions of the d.o.f. indexd by q to the stiffness matrix R for the elast...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::applyTransformationToParentGradients
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void applyTransformationToParentGradients(int const q, real64 const (&invJ)[3][3], real64(&gradN)[numNodes][3])
Apply a Jacobian transformation matrix from the parent space to the physical space on the parent shap...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcN(localIndex const q, StackVariables const &stack, real64(&N)[numNodes])
Calculate shape functions values for each support point at a quadrature point.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:272
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::jacobianTransformation
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void jacobianTransformation(int const qa, int const qb, int const qc, real64 const (&X)[8][3], real64(&J)[3][3])
Calculates the isoparametric "Jacobian" transformation matrix/mapping from the parent space to the ph...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::invJacobianTransformation
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 invJacobianTransformation(int const qa, int const qb, int const qc, real64 const (&X)[8][3], real64(&J)[3][3])
Calculates the isoparametric "Jacobian" transformation matrix/mapping from the parent space to the ph...
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:441
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::jacobianTransformationWithCorners
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void jacobianTransformationWithCorners(real64 const (&coords)[3], real64 const (&X)[8][3], real64(&J)[3][3])
Calculates the isoparametric "Jacobian" transformation matrix/mapping from the parent space to the ph...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::applyTransformationToParentGradients
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void applyTransformationToParentGradients(real64 const (&coords)[3], real64 const (&invJ)[3][3], real64(&gradN)[numNodes][3])
Apply a Jacobian transformation matrix from the parent space to the physical space on the parent shap...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::invJacobianTransformation
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 invJacobianTransformation(int const q, real64 const (&X)[8][3], real64(&J)[3][3])
Calculates the isoparametric "Jacobian" transformation matrix/mapping from the parent space to the ph...
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:461
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::plusGradNajAij
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void plusGradNajAij(int const q, real64 const (&invJ)[3][3], real64 const (&var)[6], real64(&R)[numNodes][3])
Inner product of all basis function gradients and a rank-2 symmetric tensor evaluated at a quadrature...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeBxyMatrix
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeBxyMatrix(int const qa, int const qb, int const qc, real64 const (&X)[8][3], real64(&J)[3][3], real64(&B)[6])
computes the matrix B in the case of quasi-stiffness (e.g. for pseudo-acoustic case),...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::symmetricGradient
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void symmetricGradient(int const q, real64 const (&invJ)[3][3], real64 const (&var)[numNodes][3], real64(&grad)[6])
Calculate the symmetric gradient of a vector valued support field at a quadrature point using the sto...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeStiffnesszTerm
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeStiffnesszTerm(localIndex const q, real64 const (&X)[8][3], FUNC &&func)
computes the non-zero contributions of the d.o.f. indexed by q to the partial-stiffness matrix R,...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::num1dNodes
constexpr static localIndex num1dNodes
The number of nodes/support points per element per dimension.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:64
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::jacobianTransformation
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void jacobianTransformation(real64 const (&coords)[3], real64 const (&X)[numNodes][3], real64(&J)[3][3])
Calculates the isoparametric "Jacobian" transformation matrix/mapping from the parent space to the ph...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcGradN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcGradN(localIndex const q, real64 const (&X)[numNodes][3], StackVariables const &stack, real64(&gradN)[numNodes][3])
Calculate the shape functions derivatives wrt the physical coordinates.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcGradN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcGradN(real64 const (&coords)[3], real64 const (&X)[numNodes][3], real64(&gradN)[numNodes][3])
Calculate the shape functions derivatives wrt the physical coordinates at a single point.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::computeFirstOrderStiffnessTermY
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeFirstOrderStiffnessTermY(localIndex const q, real64 const (&X)[8][3], FUNC &&func)
computes the non-zero contributions of the d.o.f. indexd by q to the y-part of the first order stiffn...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::linearIndex3DVal
DO_NOT_DOCUMENT GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE localIndex linearIndex3DVal(const localIndex qa, localIndex const qb, localIndex const qc)
The linear index associated to the given one-dimensional indices in the three directions.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:97
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::MeshData
typename Base::template MeshData< SubregionType > MeshData
Mesh data structure for the element.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:77
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::calcGradNWithCorners
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcGradNWithCorners(real64 const (&coords)[3], real64 const (&X)[8][3], real64(&gradN)[numNodes][3])
Calculate the shape functions derivatives wrt the physical coordinates at a single point.
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::gradient
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void gradient(int const q, real64 const (&invJ)[3][3], real64 const (&var)[numNodes][3], real64(&grad)[3][3])
Calculate the gradient of a vector valued support field at a point using the stored basis function gr...
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto_impl::getNumSupportPoints
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE localIndex getNumSupportPoints(StackVariables const &stack)
Get the number of support points.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:163
geos::real64
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
geos::localIndex
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
StackVariables
Kernel variables (dof numbers, jacobian and residual) located on the stack.
Generated by doxygen 1.9.1