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 final : public FiniteElementBase
45 {
46 public:
47 
49  constexpr static localIndex num1dNodes = GL_BASIS::numSupportPoints;
50 
52  constexpr static localIndex halfNodes = ( GL_BASIS::numSupportPoints - 1 )/ 2;
53 
55  constexpr static localIndex numNodes = GL_BASIS::TensorProduct3D::numSupportPoints;
56 
58  constexpr static localIndex numNodesPerFace = GL_BASIS::TensorProduct2D::numSupportPoints;
59 
61  constexpr static localIndex maxSupportPoints = numNodes;
62 
64  constexpr static localIndex numQuadraturePoints = numNodes;
65 
73  GEOS_HOST_DEVICE
74  GEOS_FORCE_INLINE
75  constexpr static localIndex linearIndex3DVal( const localIndex qa, localIndex const qb, localIndex const qc )
76  {
77  return qa + qb * num1dNodes + qc * numNodesPerFace;
78  }
79 
85  GEOS_HOST_DEVICE
86  GEOS_FORCE_INLINE
87  constexpr static localIndex meshIndexToLinearIndex3D( localIndex const k )
88  {
89  return linearIndex3DVal( ( num1dNodes - 1 ) * ( k % 2 ),
90  ( num1dNodes - 1 ) * ( ( k % 4 ) / 2 ),
91  ( num1dNodes - 1 ) * ( k / 4 ) );
92  }
93 
94 
101  GEOS_HOST_DEVICE
102  GEOS_FORCE_INLINE
103  constexpr static localIndex linearIndex2DVal( const localIndex qa, const localIndex qb )
104  {
105  return qa + qb * num1dNodes;
106  }
107 
113  GEOS_HOST_DEVICE
114  GEOS_FORCE_INLINE
115  constexpr static localIndex meshIndexToLinearIndex2D( localIndex const k )
116  {
117  return linearIndex2DVal( ( num1dNodes - 1 ) * ( k % 2 ),
118  ( num1dNodes - 1 ) * ( k / 2 ) );
119  }
120 
122  USING_FINITEELEMENTBASE
125  ~Qk_Hexahedron_Lagrange_GaussLobatto() = default;
126 
127  GEOS_HOST_DEVICE
128  virtual localIndex getNumQuadraturePoints() const override
129  {
130  return numQuadraturePoints;
131  }
132 
138  GEOS_HOST_DEVICE
139  GEOS_FORCE_INLINE
140  static localIndex getNumQuadraturePoints( StackVariables const & stack )
141  {
142  GEOS_UNUSED_VAR( stack );
143  return numQuadraturePoints;
144  }
145 
146  GEOS_HOST_DEVICE
147  GEOS_FORCE_INLINE
148  virtual localIndex getNumSupportPoints() const override
149  {
150  return numNodes;
151  }
152 
153  GEOS_HOST_DEVICE
154  GEOS_FORCE_INLINE
155  virtual localIndex getMaxSupportPoints() const override
156  {
157  return maxSupportPoints;
158  }
159 
165  GEOS_HOST_DEVICE
166  GEOS_FORCE_INLINE
167  static localIndex getNumSupportPoints( StackVariables const & stack )
168  {
169  GEOS_UNUSED_VAR( stack );
170  return numNodes;
171  }
172 
173 
179  GEOS_HOST_DEVICE
180  GEOS_FORCE_INLINE
181  static void calcN( real64 const (&coords)[3],
182  real64 (& N)[numNodes] )
183  {
184  GL_BASIS::TensorProduct3D::value( coords, N );
185  }
186 
193  GEOS_HOST_DEVICE
194  GEOS_FORCE_INLINE
195  constexpr static real64 interpolationCoord( const int q, const int k )
196  {
197  const real64 alpha = ( GL_BASIS::parentSupportCoord( q ) + 1.0 ) / 2.0;
198  return k == 0 ? ( 1.0 - alpha ) : alpha;
199  }
200 
201 
208  GEOS_HOST_DEVICE
209  GEOS_FORCE_INLINE
210  constexpr static real64 basisGradientAt( const int q, const int p )
211  {
212  if( p <= halfNodes )
213  {
214  return GL_BASIS::gradientAt( q, p );
215  }
216  else
217  {
218  return -GL_BASIS::gradientAt( GL_BASIS::numSupportPoints - 1 - q, GL_BASIS::numSupportPoints - 1 - p );
219  }
220  }
221 
233  GEOS_HOST_DEVICE
234  GEOS_FORCE_INLINE
235  constexpr static real64 jacobianCoefficient1D( const int q, const int i, const int k, const int dir )
236  {
237  if( i == dir )
238  {
239  return k== 0 ? -1.0/2.0 : 1.0/2.0;
240  }
241  else
242  {
243  return interpolationCoord( q, k );
244  }
245  }
246 
254  GEOS_HOST_DEVICE
255  GEOS_FORCE_INLINE
256  static void calcN( localIndex const q,
257  real64 (& N)[numNodes] )
258  {
259  for( int a=0; a < numNodes; ++a )
260  {
261  N[ a ] = 0;
262  }
263  N[ q ] = 1.0;
264  }
265 
274  GEOS_HOST_DEVICE
275  GEOS_FORCE_INLINE
276  static void calcN( localIndex const q,
277  StackVariables const & stack,
278  real64 ( & N )[numNodes] )
279  {
280  GEOS_UNUSED_VAR( stack );
281  return calcN( q, N );
282  }
283 
284 
294  GEOS_HOST_DEVICE
295  GEOS_FORCE_INLINE
296  static real64 calcGradN( localIndex const q,
297  real64 const (&X)[numNodes][3],
298  real64 ( &gradN )[numNodes][3] );
308  GEOS_HOST_DEVICE
309  GEOS_FORCE_INLINE
310  static real64 calcGradN( real64 const (&coords)[3],
311  real64 const (&X)[numNodes][3],
312  real64 ( &gradN )[numNodes][3] );
313 
324  GEOS_HOST_DEVICE
325  GEOS_FORCE_INLINE
326  static real64 calcGradN( localIndex const q,
327  real64 const (&X)[numNodes][3],
328  StackVariables const & stack,
329  real64 ( &gradN )[numNodes][3] );
339  GEOS_HOST_DEVICE
340  GEOS_FORCE_INLINE
341  static real64 calcGradNWithCorners( localIndex const q,
342  real64 const (&X)[8][3],
343  real64 ( &gradN )[numNodes][3] );
353  GEOS_HOST_DEVICE
354  GEOS_FORCE_INLINE
355  static real64 calcGradNWithCorners( real64 const (&coords)[3],
356  real64 const (&X)[8][3],
357  real64 ( &gradN )[numNodes][3] );
358 
369  GEOS_HOST_DEVICE
370  GEOS_FORCE_INLINE
371  static real64 calcGradNWithCorners( localIndex const q,
372  real64 const (&X)[8][3],
373  StackVariables const & stack,
374  real64 ( &gradN )[numNodes][3] );
375 
383  GEOS_HOST_DEVICE
384  GEOS_FORCE_INLINE
385  static real64 transformedQuadratureWeight( localIndex const q,
386  real64 const (&X)[numNodes][3] );
387 
396  GEOS_HOST_DEVICE
397  GEOS_FORCE_INLINE
398  static void jacobianTransformation2d( int const qa,
399  int const qb,
400  real64 const (&X)[4][3],
401  real64 ( &J )[3][2] );
402 
403 
414  GEOS_HOST_DEVICE
415  GEOS_FORCE_INLINE
416  static real64 invJacobianTransformation( int const qa,
417  int const qb,
418  int const qc,
419  real64 const (&X)[8][3],
420  real64 ( & J )[3][3] )
421  {
422  jacobianTransformation( qa, qb, qc, X, J );
423  return LvArray::tensorOps::invert< 3 >( J );
424  }
425 
434  GEOS_HOST_DEVICE
435  GEOS_FORCE_INLINE
436  static real64 invJacobianTransformation( int const q,
437  real64 const (&X)[8][3],
438  real64 ( & J )[3][3] )
439  {
440  int qa, qb, qc;
441  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
442  return invJacobianTransformation( qa, qb, qc, X, J );
443  }
444 
445 
456  GEOS_HOST_DEVICE
457  GEOS_FORCE_INLINE
458  static void symmetricGradient( int const q,
459  real64 const (&invJ)[3][3],
460  real64 const (&var)[numNodes][3],
461  real64 ( &grad )[6] );
462 
463 
464 
480  GEOS_HOST_DEVICE
481  GEOS_FORCE_INLINE
482  static void gradient( int const q,
483  real64 const (&invJ)[3][3],
484  real64 const (&var)[numNodes][3],
485  real64 ( &grad )[3][3] );
486 
487 
503  GEOS_HOST_DEVICE
504  GEOS_FORCE_INLINE
505  static void plusGradNajAij( int const q,
506  real64 const (&invJ)[3][3],
507  real64 const (&var)[6],
508  real64 ( &R )[numNodes][3] );
509 
510 
511 
521  GEOS_HOST_DEVICE
522  GEOS_FORCE_INLINE
523  static void jacobianTransformation( int const qa,
524  int const qb,
525  int const qc,
526  real64 const (&X)[8][3],
527  real64 ( &J )[3][3] );
528 
536  GEOS_HOST_DEVICE
537  GEOS_FORCE_INLINE
538  static void jacobianTransformation( real64 const (&coords)[3],
539  real64 const (&X)[numNodes][3],
540  real64 ( &J )[3][3] );
541 
551  GEOS_HOST_DEVICE
552  GEOS_FORCE_INLINE
553  static void jacobianTransformationWithCorners( real64 const (&coords)[3],
554  real64 const (&X)[8][3],
555  real64 ( &J )[3][3] );
556 
566  GEOS_HOST_DEVICE
567  GEOS_FORCE_INLINE
568  static void trilinearInterp( real64 const alpha,
569  real64 const beta,
570  real64 const gamma,
571  real64 const (&X)[8][3],
572  real64 ( &coords )[3] );
573 
579  GEOS_HOST_DEVICE
580  GEOS_FORCE_INLINE
581  static void computeLocalCoords( real64 const (&Xmesh)[8][3],
582  real64 const (&X)[numNodes][3] );
583 
591  GEOS_HOST_DEVICE
592  GEOS_FORCE_INLINE
593  static real64 computeMassTerm( localIndex const q,
594  real64 const (&X)[8][3] );
595 
604  GEOS_HOST_DEVICE
605  GEOS_FORCE_INLINE
606  static real64 computeDampingTerm( localIndex const q,
607  real64 const (&X)[4][3] );
608 
619  GEOS_HOST_DEVICE
620  GEOS_FORCE_INLINE
621  static void computeBMatrix( int const qa,
622  int const qb,
623  int const qc,
624  real64 const (&X)[8][3],
625  real64 ( &J )[3][3],
626  real64 ( &B )[6] );
627 
636  template< typename FUNC >
637  GEOS_HOST_DEVICE
638  GEOS_FORCE_INLINE
639  static void computeStiffnessTerm( localIndex const q,
640  real64 const (&X)[8][3],
641  FUNC && func );
642 
653  GEOS_HOST_DEVICE
654  GEOS_FORCE_INLINE
655  static void computeBxyMatrix( int const qa,
656  int const qb,
657  int const qc,
658  real64 const (&X)[8][3],
659  real64 ( &J )[3][3],
660  real64 ( &B )[6] );
661 
670  template< typename FUNC >
671  GEOS_HOST_DEVICE
672  GEOS_FORCE_INLINE
673  static void computeStiffnessxyTerm( localIndex const q,
674  real64 const (&X)[8][3],
675  FUNC && func );
676 
687  GEOS_HOST_DEVICE
688  GEOS_FORCE_INLINE
689  static void computeBzMatrix( int const qa,
690  int const qb,
691  int const qc,
692  real64 const (&X)[8][3],
693  real64 ( &J )[3][3],
694  real64 ( &B )[6] );
695 
704  template< typename FUNC >
705  GEOS_HOST_DEVICE
706  GEOS_FORCE_INLINE
707  static void computeStiffnesszTerm( localIndex const q,
708  real64 const (&X)[8][3],
709  FUNC && func );
710 
720  template< typename FUNC >
721  GEOS_HOST_DEVICE
722  GEOS_FORCE_INLINE
723  static void
724  computeGradPhiBGradPhi( int const qa,
725  int const qb,
726  int const qc,
727  real64 const (&B)[6],
728  FUNC && func );
729 
738  template< typename FUNC >
739  GEOS_HOST_DEVICE
740  GEOS_FORCE_INLINE
741  static void computeFirstOrderStiffnessTermX( localIndex const q,
742  real64 const (&X)[8][3],
743  FUNC && func );
752  template< typename FUNC >
753  GEOS_HOST_DEVICE
754  GEOS_FORCE_INLINE
755  static void computeFirstOrderStiffnessTermY( localIndex const q,
756  real64 const (&X)[8][3],
757  FUNC && func );
766  template< typename FUNC >
767  GEOS_HOST_DEVICE
768  GEOS_FORCE_INLINE
769  static void computeFirstOrderStiffnessTermZ( localIndex const q,
770  real64 const (&X)[8][3],
771  FUNC && func );
781  template< typename FUNC >
782  GEOS_HOST_DEVICE
783  GEOS_FORCE_INLINE
784  static void computeFirstOrderStiffnessTerm( localIndex const q,
785  real64 const (&X)[8][3],
786  FUNC && stiffnessVal );
787 
788 
798  GEOS_HOST_DEVICE
799  GEOS_FORCE_INLINE
800  static void applyTransformationToParentGradients( int const q,
801  real64 const ( &invJ )[3][3],
802  real64 ( &gradN )[numNodes][3] );
803 
813  GEOS_HOST_DEVICE
814  GEOS_FORCE_INLINE
815  static void applyTransformationToParentGradients( real64 const (&coords)[3],
816  real64 const ( &invJ )[3][3],
817  real64 ( &gradN )[numNodes][3] );
818 
819 
820 private:
822  constexpr static real64 parentLength = GL_BASIS::parentSupportCoord( 1 ) - GL_BASIS::parentSupportCoord( 0 );
823 
825  constexpr static real64 parentVolume = parentLength*parentLength*parentLength;
835  template< typename FUNC, typename ... PARAMS >
836  GEOS_HOST_DEVICE
837  GEOS_FORCE_INLINE
838  static void supportLoop( real64 const (&coords)[3],
839  FUNC && func,
840  PARAMS &&... params );
850  template< typename FUNC, typename ... PARAMS >
851  GEOS_HOST_DEVICE
852  GEOS_FORCE_INLINE
853  static void supportLoop( localIndex const q,
854  FUNC && func,
855  PARAMS &&... params );
856 
857 };
858 
860 
861 
862 template< typename GL_BASIS >
863 template< typename FUNC, typename ... PARAMS >
864 GEOS_HOST_DEVICE
865 GEOS_FORCE_INLINE
866 void
867 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::supportLoop( real64 const (&coords)[3],
868  FUNC && func,
869  PARAMS &&... params )
870 {
871  for( int c=0; c<num1dNodes; ++c )
872  {
873  for( int b=0; b<num1dNodes; ++b )
874  {
875  for( int a=0; a<num1dNodes; ++a )
876  {
877  real64 const dNdXi[3] = { GL_BASIS::gradient( a, coords[0] )*
878  GL_BASIS::value( b, coords[1] )*
879  GL_BASIS::value( c, coords[2] ),
880  GL_BASIS::value( a, coords[0] )*
881  GL_BASIS::gradient( b, coords[1] )*
882  GL_BASIS::value( c, coords[2] ),
883  GL_BASIS::value( a, coords[0] )*
884  GL_BASIS::value( b, coords[1] )*
885  GL_BASIS::gradient( c, coords[2] )};
886 
887  localIndex const nodeIndex = GL_BASIS::TensorProduct3D::linearIndex( a, b, c );
888 
889  func( dNdXi, nodeIndex, std::forward< PARAMS >( params )... );
890  }
891  }
892  }
893 }
894 
895 template< typename GL_BASIS >
896 template< typename FUNC, typename ... PARAMS >
897 GEOS_HOST_DEVICE
898 GEOS_FORCE_INLINE
899 void
900 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::supportLoop( localIndex const q,
901  FUNC && func,
902  PARAMS &&... params )
903 {
904  int qa, qb, qc;
905  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
906  for( int c=0; c<num1dNodes; ++c )
907  {
908  for( int b=0; b<num1dNodes; ++b )
909  {
910  for( int a=0; a<num1dNodes; ++a )
911  {
912  real64 const dNdXi[3] = { (b == qb && c == qc ) ? basisGradientAt( a, qa ) : 0,
913  (a == qa && c == qc ) ? basisGradientAt( b, qb ) : 0,
914  (a == qa && b == qb ) ? basisGradientAt( c, qc ) : 0 };
915 
916  localIndex const nodeIndex = GL_BASIS::TensorProduct3D::linearIndex( a, b, c );
917 
918  func( dNdXi, nodeIndex, std::forward< PARAMS >( params )... );
919  }
920  }
921  }
922 }
923 
924 //*************************************************************************************************
925 
926 template< typename GL_BASIS >
927 GEOS_HOST_DEVICE
928 GEOS_FORCE_INLINE
929 real64
930 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::calcGradN( localIndex const q,
931  real64 const (&X)[numNodes][3],
932  real64 (& gradN)[numNodes][3] )
933 {
934  int qa, qb, qc;
935  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
936  real64 Xmesh[8][3] = {{0}};
937  for( int k = 0; k < 8; k++ )
938  {
939  const localIndex nodeIndex = meshIndexToLinearIndex3D( k );
940  for( int i = 0; i < 3; i++ )
941  {
942  Xmesh[ k ][ i ] = X[ nodeIndex ][ i ];
943  }
944  }
945  real64 J[3][3] = {{0}};
946 
947  jacobianTransformation( qa, qb, qc, Xmesh, J );
948 
949  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
950 
951  applyTransformationToParentGradients( q, J, gradN );
952 
953  return detJ;
954 }
955 //*************************************************************************************************
956 template< typename GL_BASIS >
957 GEOS_HOST_DEVICE
958 GEOS_FORCE_INLINE
959 real64
960 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::calcGradN( real64 const (&coords)[3],
961  real64 const (&X)[numNodes][3],
962  real64 (& gradN)[numNodes][3] )
963 {
964  real64 J[3][3] = {{0}};
965 
966  jacobianTransformation( coords, X, J );
967 
968  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
969 
970  applyTransformationToParentGradients( coords, J, gradN );
971 
972  return detJ;
973 }
974 template< typename GL_BASIS >
975 GEOS_HOST_DEVICE
976 GEOS_FORCE_INLINE
977 real64 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
978 calcGradN( localIndex const q,
979  real64 const (&X)[numNodes][3],
980  StackVariables const & GEOS_UNUSED_PARAM( stack ),
981  real64 ( & gradN )[numNodes][3] )
982 {
983  return calcGradN( q, X, gradN );
984 }
985 
986 template< typename GL_BASIS >
987 GEOS_HOST_DEVICE
988 GEOS_FORCE_INLINE
989 real64
990 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::calcGradNWithCorners( localIndex const q,
991  real64 const (&X)[8][3],
992  real64 (& gradN)[numNodes][3] )
993 {
994  int qa, qb, qc;
995  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
996 
997  real64 J[3][3] = {{0}};
998 
999  jacobianTransformation( qa, qb, qc, X, J );
1000 
1001  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
1002 
1003  applyTransformationToParentGradients( q, J, gradN );
1004 
1005  return detJ;
1006 }
1007 //*************************************************************************************************
1008 template< typename GL_BASIS >
1009 GEOS_HOST_DEVICE
1010 GEOS_FORCE_INLINE
1011 real64
1012 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::calcGradNWithCorners( real64 const (&coords)[3],
1013  real64 const (&X)[8][3],
1014  real64 (& gradN)[numNodes][3] )
1015 {
1016  real64 J[3][3] = {{0}};
1017 
1018  jacobianTransformationWithCorners( coords, X, J );
1019 
1020  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
1021 
1022  applyTransformationToParentGradients( coords, J, gradN );
1023 
1024  return detJ;
1025 }
1026 template< typename GL_BASIS >
1027 GEOS_HOST_DEVICE
1028 GEOS_FORCE_INLINE
1029 real64 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1030 calcGradNWithCorners( localIndex const q,
1031  real64 const (&X)[8][3],
1032  StackVariables const & GEOS_UNUSED_PARAM( stack ),
1033  real64 ( & gradN )[numNodes][3] )
1034 {
1035  return calcGradN( q, X, gradN );
1036 }
1037 //*************************************************************************************************
1038 #if __GNUC__
1039 #pragma GCC diagnostic push
1040 #pragma GCC diagnostic ignored "-Wshadow"
1041 #endif
1042 
1043 template< typename GL_BASIS >
1044 GEOS_HOST_DEVICE
1045 GEOS_FORCE_INLINE
1046 void
1047 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1048 jacobianTransformation( int const qa,
1049  int const qb,
1050  int const qc,
1051  real64 const (&X)[8][3],
1052  real64 ( & J )[3][3] )
1053 {
1054  for( int k = 0; k < 8; k++ )
1055  {
1056  const int ka = k % 2;
1057  const int kb = ( k % 4 ) / 2;
1058  const int kc = k / 4;
1059  for( int j = 0; j < 3; j++ )
1060  {
1061  real64 jacCoeff = jacobianCoefficient1D( qa, 0, ka, j ) *
1062  jacobianCoefficient1D( qb, 1, kb, j ) *
1063  jacobianCoefficient1D( qc, 2, kc, j );
1064  for( int i = 0; i < 3; i++ )
1065  {
1066  J[i][j] += jacCoeff * X[k][i];
1067  }
1068  }
1069  }
1070 }
1071 
1072 template< typename GL_BASIS >
1073 GEOS_HOST_DEVICE
1074 GEOS_FORCE_INLINE
1075 void
1076 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1077 jacobianTransformation( real64 const (&coords)[3],
1078  real64 const (&X)[numNodes][3],
1079  real64 ( & J )[3][3] )
1080 {
1081  supportLoop( coords, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1082  int const nodeIndex,
1083  real64 const (&X)[numNodes][3],
1084  real64 (& J)[3][3] )
1085  {
1086  real64 const * const GEOS_RESTRICT Xnode = X[nodeIndex];
1087  for( int i = 0; i < 3; ++i )
1088  {
1089  for( int j = 0; j < 3; ++j )
1090  {
1091  J[i][j] = J[i][j] + dNdXi[ j ] * Xnode[i];
1092  }
1093  }
1094  }, X, J );
1095 }
1096 
1097 template< typename GL_BASIS >
1098 GEOS_HOST_DEVICE
1099 GEOS_FORCE_INLINE
1100 void
1101 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1102 jacobianTransformationWithCorners( real64 const (&coords)[3],
1103  real64 const (&X)[8][3],
1104  real64 ( & J )[3][3] )
1105 {
1106  supportLoop( coords, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1107  int const nodeIndex,
1108  real64 const (&X)[8][3],
1109  real64 (& J)[3][3] )
1110  {
1111  int qa, qb, qc;
1112  GL_BASIS::TensorProduct3D::multiIndex( nodeIndex, qa, qb, qc );
1113  real64 Xnode[3];
1114  real64 alpha = ( GL_BASIS::parentSupportCoord( qa ) + 1.0 ) / 2.0;
1115  real64 beta = ( GL_BASIS::parentSupportCoord( qb ) + 1.0 ) / 2.0;
1116  real64 gamma = ( GL_BASIS::parentSupportCoord( qc ) + 1.0 ) / 2.0;
1117  trilinearInterp( alpha, beta, gamma, X, Xnode );
1118  for( int i = 0; i < 3; ++i )
1119  {
1120  for( int j = 0; j < 3; ++j )
1121  {
1122  J[i][j] = J[i][j] + dNdXi[ j ] * Xnode[i];
1123  }
1124  }
1125  }, X, J );
1126 }
1127 
1128 template< typename GL_BASIS >
1129 GEOS_HOST_DEVICE
1130 GEOS_FORCE_INLINE
1131 void
1132 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1133 trilinearInterp( real64 const alpha,
1134  real64 const beta,
1135  real64 const gamma,
1136  real64 const (&X)[8][3],
1137  real64 (& coords)[3] )
1138 {
1139  for( int i=0; i<3; i++ )
1140  {
1141  coords[i] = X[0][i]*( 1.0-alpha )*( 1.0-beta )*( 1.0-gamma )+
1142  X[1][i]* alpha *( 1.0-beta )*( 1.0-gamma )+
1143  X[2][i]*( 1.0-alpha )* beta *( 1.0-gamma )+
1144  X[3][i]* alpha * beta *( 1.0-gamma )+
1145  X[4][i]*( 1.0-alpha )*( 1.0-beta )* gamma+
1146  X[5][i]* alpha *( 1.0-beta )* gamma+
1147  X[6][i]*( 1.0-alpha )* beta * gamma+
1148  X[7][i]* alpha * beta * gamma;
1149  }
1150 }
1151 
1152 
1153 template< typename GL_BASIS >
1154 GEOS_HOST_DEVICE
1155 GEOS_FORCE_INLINE
1156 void
1157 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1158 computeLocalCoords( real64 const (&Xmesh)[8][3],
1159  real64 const (&X)[numNodes][3] )
1160 {
1161  int qa, qb, qc;
1162  for( int q=0; q<numNodes; q++ )
1163  {
1164  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1165  real64 alpha = ( GL_BASIS::parentSupportCoord( qa ) + 1.0 ) / 2.0;
1166  real64 beta = ( GL_BASIS::parentSupportCoord( qb ) + 1.0 ) / 2.0;
1167  real64 gamma = ( GL_BASIS::parentSupportCoord( qc ) + 1.0 ) / 2.0;
1168  trilinearInterp( alpha, beta, gamma, Xmesh, X[q] );
1169  }
1170 }
1171 
1172 template< typename GL_BASIS >
1173 GEOS_HOST_DEVICE
1174 GEOS_FORCE_INLINE
1175 void
1176 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1177 jacobianTransformation2d( int const qa,
1178  int const qb,
1179  real64 const (&X)[4][3],
1180  real64 ( & J )[3][2] )
1181 {
1182  for( int k = 0; k < 4; k++ )
1183  {
1184  int ka = k % 2;
1185  int kb = k / 2;
1186  for( int j = 0; j < 2; j++ )
1187  {
1188  real64 jacCoeff = jacobianCoefficient1D( qa, 0, ka, j ) *
1189  jacobianCoefficient1D( qb, 1, kb, j );
1190  for( int i = 0; i < 3; i++ )
1191  {
1192  J[i][j] += jacCoeff * X[k][i];
1193  }
1194  }
1195  }
1196 }
1197 
1198 template< typename GL_BASIS >
1199 GEOS_HOST_DEVICE
1200 GEOS_FORCE_INLINE
1201 real64
1202 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1203 computeMassTerm( localIndex const q,
1204  real64 const (&X)[8][3] )
1205 {
1206  int qa, qb, qc;
1207  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1208  const real64 w3D = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1209  real64 J[3][3] = {{0}};
1210  jacobianTransformation( qa, qb, qc, X, J );
1211  return LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( J ) )*w3D;
1212 }
1213 
1214 template< typename GL_BASIS >
1215 GEOS_HOST_DEVICE
1216 GEOS_FORCE_INLINE
1217 real64
1218 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1219 computeDampingTerm( localIndex const q,
1220  real64 const (&X)[4][3] )
1221 {
1222  int qa, qb;
1223  GL_BASIS::TensorProduct2D::multiIndex( q, qa, qb );
1224  const real64 w2D = GL_BASIS::weight( qa )*GL_BASIS::weight( qb );
1225  real64 B[3];
1226  real64 J[3][2] = {{0}};
1227  jacobianTransformation2d( qa, qb, X, J );
1228  // compute J^T.J, using Voigt notation for B
1229  B[0] = J[0][0]*J[0][0]+J[1][0]*J[1][0]+J[2][0]*J[2][0];
1230  B[1] = J[0][1]*J[0][1]+J[1][1]*J[1][1]+J[2][1]*J[2][1];
1231  B[2] = J[0][0]*J[0][1]+J[1][0]*J[1][1]+J[2][0]*J[2][1];
1232  return sqrt( LvArray::math::abs( LvArray::tensorOps::symDeterminant< 2 >( B ) ) )*w2D;
1233 }
1234 
1235 template< typename GL_BASIS >
1236 GEOS_HOST_DEVICE
1237 GEOS_FORCE_INLINE
1238 void
1239 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1240 computeBMatrix( int const qa,
1241  int const qb,
1242  int const qc,
1243  real64 const (&X)[8][3],
1244  real64 (& J)[3][3],
1245  real64 (& B)[6] )
1246 {
1247  jacobianTransformation( qa, qb, qc, X, J );
1248  real64 const detJ = LvArray::tensorOps::determinant< 3 >( J );
1249 
1250  // compute J^T.J/det(J), using Voigt notation for B
1251  B[0] = (J[0][0]*J[0][0]+J[1][0]*J[1][0]+J[2][0]*J[2][0])/detJ;
1252  B[1] = (J[0][1]*J[0][1]+J[1][1]*J[1][1]+J[2][1]*J[2][1])/detJ;
1253  B[2] = (J[0][2]*J[0][2]+J[1][2]*J[1][2]+J[2][2]*J[2][2])/detJ;
1254  B[3] = (J[0][1]*J[0][2]+J[1][1]*J[1][2]+J[2][1]*J[2][2])/detJ;
1255  B[4] = (J[0][0]*J[0][2]+J[1][0]*J[1][2]+J[2][0]*J[2][2])/detJ;
1256  B[5] = (J[0][0]*J[0][1]+J[1][0]*J[1][1]+J[2][0]*J[2][1])/detJ;
1257 
1258  // compute detJ*J^{-1}J^{-T}
1259  LvArray::tensorOps::symInvert< 3 >( B );
1260 }
1261 
1262 template< typename GL_BASIS >
1263 GEOS_HOST_DEVICE
1264 GEOS_FORCE_INLINE
1265 void
1266 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1267 computeBzMatrix( int const qa,
1268  int const qb,
1269  int const qc,
1270  real64 const (&X)[8][3],
1271  real64 (& J)[3][3],
1272  real64 (& B)[6] )
1273 {
1274  jacobianTransformation( qa, qb, qc, X, J );
1275  real64 const detJ = LvArray::tensorOps::determinant< 3 >( J );
1276 
1277  real64 Jinv[3][3] = {{0}};
1278  LvArray::tensorOps::invert< 3 >( Jinv, J );
1279 
1280  // compute det(J)*J^{-1}Az*J^{-T}, using Voigt notation for B
1281  B[0] = detJ*(Jinv[0][2]*Jinv[0][2]);
1282  B[1] = detJ*(Jinv[1][2]*Jinv[1][2]);
1283  B[2] = detJ*(Jinv[2][2]*Jinv[2][2]);
1284  B[3] = detJ*(Jinv[1][2]*Jinv[2][2]);
1285  B[4] = detJ*(Jinv[0][2]*Jinv[2][2]);
1286  B[5] = detJ*(Jinv[0][2]*Jinv[1][2]);
1287 }
1288 
1289 template< typename GL_BASIS >
1290 GEOS_HOST_DEVICE
1291 GEOS_FORCE_INLINE
1292 void
1293 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1294 computeBxyMatrix( int const qa,
1295  int const qb,
1296  int const qc,
1297  real64 const (&X)[8][3],
1298  real64 (& J)[3][3],
1299  real64 (& B)[6] )
1300 {
1301  jacobianTransformation( qa, qb, qc, X, J );
1302  real64 const detJ = LvArray::tensorOps::determinant< 3 >( J );
1303 
1304  real64 Jinv[3][3] = {{0}};
1305  LvArray::tensorOps::invert< 3 >( Jinv, J );
1306 
1307  // compute det(J)*J^{-1}Axy*J^{-T}, using Voigt notation for B
1308  B[0] = detJ*(Jinv[0][0]*Jinv[0][0] + Jinv[0][1]*Jinv[0][1]);
1309  B[1] = detJ*(Jinv[1][1]*Jinv[1][1] + Jinv[1][0]*Jinv[1][0]);
1310  B[2] = detJ*(Jinv[2][0]*Jinv[2][0] + Jinv[2][1]*Jinv[2][1]);
1311  B[3] = detJ*(Jinv[1][0]*Jinv[2][0] + Jinv[1][1]*Jinv[2][1]);
1312  B[4] = detJ*(Jinv[0][0]*Jinv[2][0] + Jinv[0][1]*Jinv[2][1]);
1313  B[5] = detJ*(Jinv[0][0]*Jinv[1][0] + Jinv[0][1]*Jinv[1][1]);
1314 }
1315 
1316 template< typename GL_BASIS >
1317 template< typename FUNC >
1318 GEOS_HOST_DEVICE
1319 GEOS_FORCE_INLINE
1320 void
1321 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1322 computeGradPhiBGradPhi( int const qa,
1323  int const qb,
1324  int const qc,
1325  real64 const (&B)[6],
1326  FUNC && func )
1327 {
1328  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1329  for( int i=0; i<num1dNodes; i++ )
1330  {
1331  const int ibc = GL_BASIS::TensorProduct3D::linearIndex( i, qb, qc );
1332  const int aic = GL_BASIS::TensorProduct3D::linearIndex( qa, i, qc );
1333  const int abi = GL_BASIS::TensorProduct3D::linearIndex( qa, qb, i );
1334  const real64 gia = basisGradientAt( i, qa );
1335  const real64 gib = basisGradientAt( i, qb );
1336  const real64 gic = basisGradientAt( i, qc );
1337  for( int j=0; j<num1dNodes; j++ )
1338  {
1339  const int jbc = GL_BASIS::TensorProduct3D::linearIndex( j, qb, qc );
1340  const int ajc = GL_BASIS::TensorProduct3D::linearIndex( qa, j, qc );
1341  const int abj = GL_BASIS::TensorProduct3D::linearIndex( qa, qb, j );
1342  const real64 gja = basisGradientAt( j, qa );
1343  const real64 gjb = basisGradientAt( j, qb );
1344  const real64 gjc = basisGradientAt( j, qc );
1345  // diagonal terms
1346  const real64 w0 = w * gia * gja;
1347  func( ibc, jbc, w0 * B[0] );
1348  const real64 w1 = w * gib * gjb;
1349  func( aic, ajc, w1 * B[1] );
1350  const real64 w2 = w * gic * gjc;
1351  func( abi, abj, w2 * B[2] );
1352  // off-diagonal terms
1353  const real64 w3 = w * gib * gjc;
1354  func( aic, abj, w3 * B[3] );
1355  func( abj, aic, w3 * B[3] );
1356  const real64 w4 = w * gia * gjc;
1357  func( ibc, abj, w4 * B[4] );
1358  func( abj, ibc, w4 * B[4] );
1359  const real64 w5 = w * gia * gjb;
1360  func( ibc, ajc, w5 * B[5] );
1361  func( ajc, ibc, w5 * B[5] );
1362  }
1363  }
1364 }
1365 
1366 template< typename GL_BASIS >
1367 template< typename FUNC >
1368 GEOS_HOST_DEVICE
1369 GEOS_FORCE_INLINE
1370 void
1371 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1372 computeStiffnessxyTerm( localIndex const q,
1373  real64 const (&X)[8][3],
1374  FUNC && func )
1375 {
1376  int qa, qb, qc;
1377  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1378  real64 B[6] = {0};
1379  real64 J[3][3] = {{0}};
1380  computeBxyMatrix( qa, qb, qc, X, J, B ); // The only change!
1381  computeGradPhiBGradPhi( qa, qb, qc, B, func );
1382 }
1383 
1384 template< typename GL_BASIS >
1385 template< typename FUNC >
1386 GEOS_HOST_DEVICE
1387 GEOS_FORCE_INLINE
1388 void
1389 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1390 computeStiffnesszTerm( localIndex const q,
1391  real64 const (&X)[8][3],
1392  FUNC && func )
1393 {
1394  int qa, qb, qc;
1395  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1396  real64 B[6] = {0};
1397  real64 J[3][3] = {{0}};
1398  computeBzMatrix( qa, qb, qc, X, J, B ); // The only change!
1399  computeGradPhiBGradPhi( qa, qb, qc, B, func );
1400 }
1401 
1402 template< typename GL_BASIS >
1403 template< typename FUNC >
1404 GEOS_HOST_DEVICE
1405 GEOS_FORCE_INLINE
1406 void
1407 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1408 computeStiffnessTerm( localIndex const q,
1409  real64 const (&X)[8][3],
1410  FUNC && func )
1411 {
1412  int qa, qb, qc;
1413  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1414  real64 B[6] = {0};
1415  real64 J[3][3] = {{0}};
1416  computeBMatrix( qa, qb, qc, X, J, B );
1417  computeGradPhiBGradPhi( qa, qb, qc, B, func );
1418 }
1419 
1420 template< typename GL_BASIS >
1421 template< typename FUNC >
1422 GEOS_HOST_DEVICE
1423 GEOS_FORCE_INLINE
1424 void
1425 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1426 computeFirstOrderStiffnessTerm( localIndex const q,
1427  real64 const (&X)[8][3],
1428  FUNC && func )
1429 {
1430  int qa, qb, qc;
1431  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1432  real64 J[3][3] = {{0}};
1433  jacobianTransformation( qa, qb, qc, X, J );
1434  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
1435  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1436  for( int i=0; i<num1dNodes; i++ )
1437  {
1438  const int ibc = GL_BASIS::TensorProduct3D::linearIndex( i, qb, qc );
1439  const int aic = GL_BASIS::TensorProduct3D::linearIndex( qa, i, qc );
1440  const int abi = GL_BASIS::TensorProduct3D::linearIndex( qa, qb, i );
1441  const real64 gia = basisGradientAt( i, qa );
1442  const real64 gib = basisGradientAt( i, qb );
1443  const real64 gic = basisGradientAt( i, qc );
1444  for( int j=0; j<num1dNodes; j++ )
1445  {
1446  const int jbc = GL_BASIS::TensorProduct3D::linearIndex( j, qb, qc );
1447  const int ajc = GL_BASIS::TensorProduct3D::linearIndex( qa, j, qc );
1448  const int abj = GL_BASIS::TensorProduct3D::linearIndex( qa, qb, j );
1449  const real64 gja = basisGradientAt( j, qa );
1450  const real64 gjb = basisGradientAt( j, qb );
1451  const real64 gjc = basisGradientAt( j, qc );
1452  // diagonal terms
1453  const real64 w00 = w * gia * gja;
1454  func( ibc, jbc, w00 * detJ, J, 0, 0 );
1455  const real64 w11 = w * gib * gjb;
1456  func( aic, ajc, w11 * detJ, J, 1, 1 );
1457  const real64 w22 = w * gic * gjc;
1458  func( abi, abj, w22 * detJ, J, 2, 2 );
1459  // off-diagonal terms
1460  const real64 w12 = w * gib * gjc;
1461  func( aic, abj, w12 * detJ, J, 1, 2 );
1462  func( abj, aic, w12 * detJ, J, 2, 1 );
1463  const real64 w02 = w * gia * gjc;
1464  func( ibc, abj, w02 * detJ, J, 0, 2 );
1465  func( abj, ibc, w02 * detJ, J, 2, 0 );
1466  const real64 w01 = w * gia * gjb;
1467  func( ibc, ajc, w01 * detJ, J, 0, 1 );
1468  func( ajc, ibc, w01 * detJ, J, 1, 0 );
1469  }
1470  }
1471 }
1472 
1473 template< typename GL_BASIS >
1474 template< typename FUNC >
1475 GEOS_HOST_DEVICE
1476 GEOS_FORCE_INLINE
1477 void
1478 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1479 computeFirstOrderStiffnessTermX( localIndex const q,
1480  real64 const (&X)[8][3],
1481  FUNC && func )
1482 {
1483  int qa, qb, qc;
1484  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1485  real64 J[3][3] = {{0}};
1486  jacobianTransformation( qa, qb, qc, X, J );
1487  const real64 detJ = LvArray::tensorOps::invert< 3 >( J );
1488  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1489 
1490  for( int i1 = 0; i1 < num1dNodes; ++i1 )
1491  {
1492  auto val = w * basisGradientAt( i1, qa );
1493  func( GL_BASIS::TensorProduct3D::linearIndex( i1, qb, qc ), q, detJ*J[0][0]*val, detJ*J[0][1]*val, detJ*J[0][2]*val );
1494  }
1495 
1496 }
1497 
1498 template< typename GL_BASIS >
1499 template< typename FUNC >
1500 GEOS_HOST_DEVICE
1501 GEOS_FORCE_INLINE
1502 void
1503 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1504 computeFirstOrderStiffnessTermY( localIndex const q,
1505  real64 const (&X)[8][3],
1506  FUNC && func )
1507 {
1508  int qa, qb, qc;
1509  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1510  real64 J[3][3] = {{0}};
1511  jacobianTransformation( qa, qb, qc, X, J );
1512  const real64 detJ = LvArray::tensorOps::invert< 3 >( J );
1513  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1514 
1515  for( int i2 = 0; i2 < num1dNodes; ++i2 )
1516  {
1517  auto val = w * basisGradientAt( i2, qb );
1518  func( GL_BASIS::TensorProduct3D::linearIndex( qa, i2, qc ), q, detJ*J[1][0]*val, detJ*J[1][1]*val, detJ*J[1][2]*val );
1519  }
1520 }
1521 
1522 template< typename GL_BASIS >
1523 template< typename FUNC >
1524 GEOS_HOST_DEVICE
1525 GEOS_FORCE_INLINE
1526 void
1527 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1528 computeFirstOrderStiffnessTermZ( localIndex const q,
1529  real64 const (&X)[8][3],
1530  FUNC && func )
1531 {
1532  int qa, qb, qc;
1533  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1534  real64 J[3][3] = {{0}};
1535  jacobianTransformation( qa, qb, qc, X, J );
1536  const real64 detJ = LvArray::tensorOps::invert< 3 >( J );
1537  const real64 w = GL_BASIS::weight( qa )*GL_BASIS::weight( qb )*GL_BASIS::weight( qc );
1538 
1539  for( int i3 = 0; i3 < num1dNodes; ++i3 )
1540  {
1541  auto val = w * basisGradientAt( i3, qc );
1542  func( GL_BASIS::TensorProduct3D::linearIndex( qa, qb, i3 ), q, detJ*J[2][0]*val, detJ*J[2][1]*val, detJ*J[2][2]*val );
1543  }
1544 }
1545 
1546 //*************************************************************************************************
1547 template< typename GL_BASIS >
1548 GEOS_HOST_DEVICE
1549 GEOS_FORCE_INLINE
1550 void
1551 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1552 applyTransformationToParentGradients( int const q,
1553  real64 const ( &invJ )[3][3],
1554  real64 (& gradN)[numNodes][3] )
1555 {
1556  supportLoop( q, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1557  int const nodeIndex,
1558  real64 const (&invJ)[3][3],
1559  real64 (& gradN)[numNodes][3] )
1560  {
1561  // smaller register footprint by manually unrolling the for loops.
1562  gradN[nodeIndex][0] = dNdXi[0] * invJ[0][0] + dNdXi[1] * invJ[1][0] + dNdXi[2] * invJ[2][0];
1563  gradN[nodeIndex][1] = dNdXi[0] * invJ[0][1] + dNdXi[1] * invJ[1][1] + dNdXi[2] * invJ[2][1];
1564  gradN[nodeIndex][2] = dNdXi[0] * invJ[0][2] + dNdXi[1] * invJ[1][2] + dNdXi[2] * invJ[2][2];
1565 
1566 
1567  }, invJ, gradN );
1568 }
1569 
1570 //*************************************************************************************************
1571 template< typename GL_BASIS >
1572 GEOS_HOST_DEVICE
1573 GEOS_FORCE_INLINE
1574 void
1575 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1576 applyTransformationToParentGradients( real64 const (&coords)[3],
1577  real64 const ( &invJ )[3][3],
1578  real64 (& gradN)[numNodes][3] )
1579 {
1580  supportLoop( coords, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1581  int const nodeIndex,
1582  real64 const (&invJ)[3][3],
1583  real64 (& gradN)[numNodes][3] )
1584  {
1585  gradN[nodeIndex][0] = dNdXi[0] * invJ[0][0] + dNdXi[1] * invJ[1][0] + dNdXi[2] * invJ[2][0];
1586  gradN[nodeIndex][1] = dNdXi[0] * invJ[0][1] + dNdXi[1] * invJ[1][1] + dNdXi[2] * invJ[2][1];
1587  gradN[nodeIndex][2] = dNdXi[0] * invJ[0][2] + dNdXi[1] * invJ[1][2] + dNdXi[2] * invJ[2][2];
1588  }, invJ, gradN );
1589 }
1590 
1591 template< typename GL_BASIS >
1592 GEOS_HOST_DEVICE
1593 GEOS_FORCE_INLINE
1594 real64
1595 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1596 transformedQuadratureWeight( localIndex const q,
1597  real64 const (&X)[numNodes][3] )
1598 {
1599  int qa, qb, qc;
1600  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
1601  real64 J[3][3] = {{0}};
1602 
1603  jacobianTransformation( qa, qb, qc, X, J );
1604 
1605  return LvArray::tensorOps::determinant< 3 >( J );
1606 }
1607 
1608 
1609 
1610 template< typename GL_BASIS >
1611 GEOS_HOST_DEVICE
1612 GEOS_FORCE_INLINE
1613 void Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1614 symmetricGradient( int const q,
1615  real64 const (&invJ)[3][3],
1616  real64 const (&var)[numNodes][3],
1617  real64 (& grad)[6] )
1618 {
1619  supportLoop( q, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1620  int const nodeIndex,
1621  real64 const (&invJ)[3][3],
1622  real64 const (&var)[numNodes][3],
1623  real64 (& grad)[6] )
1624  {
1625 
1626  real64 gradN[3] = {0, 0, 0};
1627  for( int i = 0; i < 3; ++i )
1628  {
1629  for( int j = 0; j < 3; ++j )
1630  {
1631  gradN[i] = gradN[i] + dNdXi[ j ] * invJ[j][i];
1632  }
1633  }
1634 
1635  grad[0] = grad[0] + gradN[0] * var[ nodeIndex ][0];
1636  grad[1] = grad[1] + gradN[1] * var[ nodeIndex ][1];
1637  grad[2] = grad[2] + gradN[2] * var[ nodeIndex ][2];
1638  grad[3] = grad[3] + gradN[2] * var[ nodeIndex ][1] + gradN[1] * var[ nodeIndex ][2];
1639  grad[4] = grad[4] + gradN[2] * var[ nodeIndex ][0] + gradN[0] * var[ nodeIndex ][2];
1640  grad[5] = grad[5] + gradN[1] * var[ nodeIndex ][0] + gradN[0] * var[ nodeIndex ][1];
1641  }, invJ, var, grad );
1642 }
1643 
1644 template< typename GL_BASIS >
1645 GEOS_HOST_DEVICE
1646 GEOS_FORCE_INLINE
1647 void Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1648 plusGradNajAij( int const q,
1649  real64 const (&invJ)[3][3],
1650  real64 const (&var)[6],
1651  real64 (& R)[numNodes][3] )
1652 {
1653  supportLoop( q,
1654  [] GEOS_HOST_DEVICE
1655  ( real64 const (&dNdXi)[3],
1656  int const nodeIndex,
1657  real64 const (&invJ)[3][3],
1658  real64 const (&var)[6],
1659  real64 (& R)[numNodes][3] )
1660  {
1661 
1662  real64 gradN[3] = {0, 0, 0};
1663  for( int i = 0; i < 3; ++i )
1664  {
1665  for( int j = 0; j < 3; ++j )
1666  {
1667  gradN[i] = gradN[i] + dNdXi[ j ] * invJ[j][i];
1668  }
1669  }
1670  R[ nodeIndex ][ 0 ] = R[ nodeIndex ][ 0 ] - var[ 0 ] * gradN[ 0 ] - var[ 5 ] * gradN[ 1 ] - var[ 4 ] * gradN[ 2 ];
1671  R[ nodeIndex ][ 1 ] = R[ nodeIndex ][ 1 ] - var[ 5 ] * gradN[ 0 ] - var[ 1 ] * gradN[ 1 ] - var[ 3 ] * gradN[ 2 ];
1672  R[ nodeIndex ][ 2 ] = R[ nodeIndex ][ 2 ] - var[ 4 ] * gradN[ 0 ] - var[ 3 ] * gradN[ 1 ] - var[ 2 ] * gradN[ 2 ];
1673  }, invJ, var, R );
1674 }
1675 
1676 
1677 
1678 template< typename GL_BASIS >
1679 GEOS_HOST_DEVICE
1680 GEOS_FORCE_INLINE
1681 void Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
1682 gradient( int const q,
1683  real64 const (&invJ)[3][3],
1684  real64 const (&var)[numNodes][3],
1685  real64 (& grad)[3][3] )
1686 {
1687  supportLoop( q, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
1688  int const nodeIndex,
1689  real64 const (&invJ)[3][3],
1690  real64 const (&var)[numNodes][3],
1691  real64 (& grad)[3][3] )
1692  {
1693  for( int i = 0; i < 3; ++i )
1694  {
1695  real64 gradN=0.0;;
1696  for( int j = 0; j < 3; ++j )
1697  {
1698  gradN = gradN + dNdXi[ j ] * invJ[j][i];
1699  }
1700  for( int k = 0; k < 3; ++k )
1701  {
1702  grad[k][i] = grad[k][i] + gradN * var[ nodeIndex ][k];
1703  }
1704  }
1705  }, invJ, var, grad );
1706 }
1733 using Q1_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis1 >;
1776 using Q2_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis2 >;
1820 using Q3_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis3GL >;
1858 using Q4_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis4GL >;
1910 using Q5_Hexahedron_Lagrange_GaussLobatto = Qk_Hexahedron_Lagrange_GaussLobatto< LagrangeBasis5GL >;
1911 
1913 
1914 #if __GNUC__
1915 #pragma GCC diagnostic pop
1916 #endif
1917 #undef PARENT_GRADIENT_METHOD
1918 }
1919 }
1920 
1921 #endif //FINITE_ELEMENT_SHAPE_KERNEL
USING_FINITEELEMENTBASE
#define USING_FINITEELEMENTBASE
Macro to simplify name resolution in derived classes.
Definition: FiniteElementBase.hpp:1076
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
Base class for FEM element implementations.
Definition: FiniteElementBase.hpp:46
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:45
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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:436
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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:115
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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::maxSupportPoints
constexpr static localIndex maxSupportPoints
The maximum number of support points per element.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:61
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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:181
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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:210
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::~Qk_Hexahedron_Lagrange_GaussLobatto
~Qk_Hexahedron_Lagrange_GaussLobatto()=default
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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::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:103
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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::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::halfNodes
constexpr static localIndex halfNodes
Half the number of support points, rounded down. Precomputed for efficiency.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:52
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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:140
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::getMaxSupportPoints
GEOS_HOST_DEVICE virtual GEOS_FORCE_INLINE localIndex getMaxSupportPoints() const override
Get the maximum number of support points for this element.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:155
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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::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:195
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::linearIndex3DVal
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:75
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::num1dNodes
constexpr static localIndex num1dNodes
The number of nodes/support points per element per dimension.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:49
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::getNumQuadraturePoints
virtual GEOS_HOST_DEVICE localIndex getNumQuadraturePoints() const override
Virtual getter for the number of quadrature points per element.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:128
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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::numQuadraturePoints
constexpr static localIndex numQuadraturePoints
The number of quadrature points per element.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:64
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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:276
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::getNumSupportPoints
GEOS_HOST_DEVICE virtual GEOS_FORCE_INLINE localIndex getNumSupportPoints() const override
Virtual getter for the number of support points per element.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:148
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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::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::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::numNodesPerFace
constexpr static localIndex numNodesPerFace
The number of nodes/support points per face.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:58
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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::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:256
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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::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:167
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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:235
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::numNodes
constexpr static localIndex numNodes
The number of nodes/support points per element.
Definition: Qk_Hexahedron_Lagrange_GaussLobatto.hpp:55
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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:416
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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:87
geos::finiteElement::Qk_Hexahedron_Lagrange_GaussLobatto::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::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::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::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::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::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::real64
double real64
64-bit floating point type.
Definition: DataTypes.hpp:99
geos::localIndex
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:85
Generated by doxygen 1.9.1