GEOS
FiniteElementBase.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 #if defined(GEOS_USE_DEVICE)
21 #define CALC_FEM_SHAPE_IN_KERNEL
22 #endif
23 
24 
25 
26 #ifndef GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_FINITEELEMENTBASE_HPP_
27 #define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_FINITEELEMENTBASE_HPP_
28 
29 #include "common/DataTypes.hpp"
30 #include "common/GeosxMacros.hpp"
31 #include "finiteElement/PDEUtilities.hpp"
32 #include "LvArray/src/tensorOps.hpp"
33 #include "mesh/NodeManager.hpp"
34 #include "mesh/EdgeManager.hpp"
35 #include "mesh/FaceManager.hpp"
36 
37 namespace geos
38 {
39 namespace finiteElement
40 {
41 
45 class FiniteElementBase
46 {
47 public:
48 
50  FiniteElementBase() = default;
51 
52 
54  constexpr static int numSamplingPointsPerDirection = 10;
55 
60  GEOS_HOST_DEVICE
61  FiniteElementBase( FiniteElementBase const & source ):
62 #ifdef CALC_FEM_SHAPE_IN_KERNEL
63  m_viewGradN(),
64  m_viewDetJ()
65 #else
66  m_viewGradN( source.m_viewGradN ),
67  m_viewDetJ( source.m_viewDetJ )
68 #endif
69  {
70  GEOS_UNUSED_VAR( source ); // suppress warning when CALC_FEM_SHAPE_IN_KERNEL is defined
71  }
72 
74  FiniteElementBase( FiniteElementBase && ) = default;
75 
80  FiniteElementBase & operator=( FiniteElementBase const & ) = delete;
81 
86  FiniteElementBase & operator=( FiniteElementBase && ) = delete;
87 
91  GEOS_HOST_DEVICE
92  virtual ~FiniteElementBase()
93  {}
94 
102  struct StackVariables
103  {};
104 
109  template< typename SUBREGION_TYPE >
110  struct MeshData
111  {};
112 
113 
122  template< typename SUBREGION_TYPE >
123  static void fillMeshData( NodeManager const & nodeManager,
124  EdgeManager const & edgeManager,
125  FaceManager const & faceManager,
126  SUBREGION_TYPE const & cellSubRegion,
127  MeshData< SUBREGION_TYPE > & meshData )
128  {
129  GEOS_UNUSED_VAR( nodeManager,
130  edgeManager,
131  faceManager,
132  cellSubRegion,
133  meshData );
134  }
135 
146  template< typename LEAF, typename SUBREGION_TYPE >
147  static void initialize( NodeManager const & nodeManager,
148  EdgeManager const & edgeManager,
149  FaceManager const & faceManager,
150  SUBREGION_TYPE const & cellSubRegion,
151  typename LEAF::template MeshData< SUBREGION_TYPE > & meshData
152  )
153  {
154  LEAF::template fillMeshData< SUBREGION_TYPE >( nodeManager,
155  edgeManager,
156  faceManager,
157  cellSubRegion,
158  meshData );
159  }
160 
161 
168  template< typename SUBREGION_TYPE >
169  GEOS_HOST_DEVICE
170  GEOS_FORCE_INLINE
171  static void setupStack( localIndex const & cellIndex,
172  MeshData< SUBREGION_TYPE > const & meshData,
173  StackVariables & stack )
174  {
175  GEOS_UNUSED_VAR( cellIndex,
176  meshData,
177  stack );
178 
179  }
180 
189  template< typename LEAF, typename SUBREGION_TYPE >
190  GEOS_HOST_DEVICE
191  void setup( localIndex const & cellIndex,
192  typename LEAF::template MeshData< SUBREGION_TYPE > const & meshData,
193  typename LEAF::StackVariables & stack ) const
194  {
195  LEAF::setupStack( cellIndex, meshData, stack );
196  }
197 
202  GEOS_HOST_DEVICE
203  virtual localIndex getNumQuadraturePoints() const = 0;
204 
209  GEOS_HOST_DEVICE
210  virtual localIndex getNumSupportPoints() const = 0;
211 
217  template< int N >
218  struct FunctionSpaceHelper
219  {};
220 
227  template< int N >
228  GEOS_HOST_DEVICE
229  constexpr static PDEUtilities::FunctionSpace getFunctionSpace();
230 
237  template< typename LEAF >
238  GEOS_HOST_DEVICE
239  localIndex numSupportPoints( typename LEAF::StackVariables const & stack ) const
240  {
241  return LEAF::getNumSupportPoints( stack );
242  }
243 
250  GEOS_HOST_DEVICE
251  virtual localIndex getMaxSupportPoints() const = 0;
252 
264  template< typename LEAF >
265  GEOS_HOST_DEVICE
266  real64 getGradN( localIndex const k,
267  localIndex const q,
268  real64 const (&X)[LEAF::maxSupportPoints][3],
269  real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
270 
283  template< typename LEAF >
284  GEOS_HOST_DEVICE
285  real64 getGradN( localIndex const k,
286  localIndex const q,
287  real64 const (&X)[LEAF::maxSupportPoints][3],
288  typename LEAF::StackVariables const & stack,
289  real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
290 
302  template< typename LEAF >
303  GEOS_HOST_DEVICE
304  real64 getGradN( localIndex const k,
305  localIndex const q,
306  int const X,
307  real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
320  template< typename LEAF >
321  GEOS_HOST_DEVICE
322  real64 getGradN( localIndex const k,
323  localIndex const q,
324  int const X,
325  typename LEAF::StackVariables const & stack,
326  real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
327 
328 
339  template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, localIndex MAXSUPPORTPOINTS, bool UPPER >
340  GEOS_HOST_DEVICE
341  GEOS_FORCE_INLINE
342  static void addGradGradStabilization( StackVariables const & stack,
343  real64 ( & matrix )
344  [MAXSUPPORTPOINTS * NUMDOFSPERTRIALSUPPORTPOINT]
345  [MAXSUPPORTPOINTS * NUMDOFSPERTRIALSUPPORTPOINT],
346  real64 const & scaleFactor )
347  {
348  GEOS_UNUSED_VAR( stack,
349  matrix,
350  scaleFactor );
351  }
352 
353 
363  template< typename LEAF, localIndex NUMDOFSPERTRIALSUPPORTPOINT, bool UPPER = false >
364  GEOS_HOST_DEVICE
365  void addGradGradStabilizationMatrix( typename LEAF::StackVariables const & stack,
366  real64 ( & matrix )
367  [LEAF::maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT]
368  [LEAF::maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT],
369  real64 const scaleFactor = 1.0 ) const
370  {
371  LEAF::template addGradGradStabilization< NUMDOFSPERTRIALSUPPORTPOINT,
372  LEAF::maxSupportPoints,
373  UPPER >( stack,
374  matrix,
375  scaleFactor );
376  }
377 
391  template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, localIndex MAXSUPPORTPOINTS >
392  GEOS_HOST_DEVICE
393  GEOS_FORCE_INLINE
394  static void addEvaluatedGradGradStabilization( StackVariables const & stack,
395  real64 const ( &dofs )[MAXSUPPORTPOINTS][NUMDOFSPERTRIALSUPPORTPOINT],
396  real64 ( & targetVector )[MAXSUPPORTPOINTS][NUMDOFSPERTRIALSUPPORTPOINT],
397  real64 const scaleFactor )
398  {
399  GEOS_UNUSED_VAR( stack );
400  GEOS_UNUSED_VAR( dofs );
401  GEOS_UNUSED_VAR( targetVector );
402  GEOS_UNUSED_VAR( scaleFactor );
403  }
404 
418  template< typename LEAF, localIndex NUMDOFSPERTRIALSUPPORTPOINT >
419  GEOS_HOST_DEVICE
420  GEOS_FORCE_INLINE
421  void
422  addEvaluatedGradGradStabilizationVector( typename LEAF::StackVariables const & stack,
423  real64 const ( &dofs )[LEAF::maxSupportPoints]
424  [NUMDOFSPERTRIALSUPPORTPOINT],
425  real64 ( & targetVector )[LEAF::maxSupportPoints]
426  [NUMDOFSPERTRIALSUPPORTPOINT],
427  real64 const scaleFactor = 1.0 ) const
428  {
429  LEAF::template
430  addEvaluatedGradGradStabilization< NUMDOFSPERTRIALSUPPORTPOINT, LEAF::maxSupportPoints >( stack,
431  dofs,
432  targetVector,
433  scaleFactor );
434  }
435 
440 
460  template< int NUM_SUPPORT_POINTS >
461  GEOS_HOST_DEVICE
462  static
463  void value( real64 const (&N)[NUM_SUPPORT_POINTS],
464  real64 const (&var)[NUM_SUPPORT_POINTS],
465  real64 & value );
466 
472  template< int NUM_SUPPORT_POINTS,
473  int NUM_COMPONENTS >
474  GEOS_HOST_DEVICE
475  static
476  void value( real64 const (&N)[NUM_SUPPORT_POINTS],
477  real64 const (&var)[NUM_SUPPORT_POINTS][NUM_COMPONENTS],
478  real64 ( &value )[NUM_COMPONENTS] );
479 
481 
486 
505  template< int NUM_SUPPORT_POINTS,
506  typename GRADIENT_TYPE >
507  GEOS_HOST_DEVICE
508  static void symmetricGradient( GRADIENT_TYPE const & gradN,
509  real64 const (&var)[NUM_SUPPORT_POINTS][3],
510  real64 ( &gradVar )[6] );
511 
524  template< int NUM_SUPPORT_POINTS,
525  typename GRADIENT_TYPE >
526  GEOS_HOST_DEVICE
527  static real64 symmetricGradientTrace( GRADIENT_TYPE const & gradN,
528  real64 const (&var)[NUM_SUPPORT_POINTS][3] );
529 
546  template< int NUM_SUPPORT_POINTS,
547  typename GRADIENT_TYPE >
548  GEOS_HOST_DEVICE
549  static void gradient( GRADIENT_TYPE const & gradN,
550  real64 const (&var)[NUM_SUPPORT_POINTS],
551  real64 ( &gradVar )[3] );
552 
563  template< int NUM_SUPPORT_POINTS,
564  typename GRADIENT_TYPE >
565  GEOS_HOST_DEVICE
566  static void gradient( GRADIENT_TYPE const & gradN,
567  real64 const (&var)[NUM_SUPPORT_POINTS][3],
568  real64 ( &gradVar )[3][3] );
570 
575 
589  template< int NUM_SUPPORT_POINTS,
590  typename GRADIENT_TYPE >
591  GEOS_HOST_DEVICE
592  static void valueAndGradient( real64 const (&N)[NUM_SUPPORT_POINTS],
593  GRADIENT_TYPE const & gradN,
594  real64 const (&var)[NUM_SUPPORT_POINTS],
595  real64 & value,
596  real64 ( &gradVar )[3] );
597 
599 
607 
628  template< int NUM_SUPPORT_POINTS,
629  typename GRADIENT_TYPE >
630  GEOS_HOST_DEVICE
631  static void plusGradNajAij( GRADIENT_TYPE const & gradN,
632  real64 const (&var_detJxW)[6],
633  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
634 
640  template< int NUM_SUPPORT_POINTS,
641  typename GRADIENT_TYPE >
642  GEOS_HOST_DEVICE
643  static void plusGradNajAij( GRADIENT_TYPE const & gradN,
644  real64 const (&var_detJxW)[3][3],
645  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
646 
655  template< int NUM_SUPPORT_POINTS >
656  GEOS_HOST_DEVICE
657  static void plusNaFi( real64 const (&N)[NUM_SUPPORT_POINTS],
658  real64 const (&forcingTerm_detJxW)[3],
659  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
660 
678  template< int NUM_SUPPORT_POINTS,
679  typename GRADIENT_TYPE >
680  GEOS_HOST_DEVICE
681  static void plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
682  real64 const (&var_detJxW)[3][3],
683  real64 const (&N)[NUM_SUPPORT_POINTS],
684  real64 const (&forcingTerm_detJxW)[3],
685  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
686 
692  template< int NUM_SUPPORT_POINTS,
693  typename GRADIENT_TYPE >
694  GEOS_HOST_DEVICE
695  static void plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
696  real64 const (&var_detJxW)[6],
697  real64 const (&N)[NUM_SUPPORT_POINTS],
698  real64 const (&forcingTerm_detJxW)[3],
699  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
700 
701 
706  void setGradNView( arrayView4d< real64 const > const & source )
707  {
708  GEOS_ERROR_IF_NE_MSG( source.size( 1 ),
709  getNumQuadraturePoints(),
710  "2nd-dimension of gradN array does not match number of quadrature points" );
711  GEOS_ERROR_IF_NE_MSG( source.size( 2 ),
712  getMaxSupportPoints(),
713  "3rd-dimension of gradN array does not match number of support points" );
714  GEOS_ERROR_IF_NE_MSG( source.size( 3 ),
715  3,
716  "4th-dimension of gradN array does not match 3" );
717 
718  m_viewGradN = source;
719  }
720 
725  void setDetJView( arrayView2d< real64 const > const & source )
726  {
727  GEOS_ERROR_IF_NE_MSG( source.size( 1 ),
728  getNumQuadraturePoints(),
729  "2nd-dimension of gradN array does not match number of quadrature points" );
730  m_viewDetJ = source;
731  }
732 
737  arrayView4d< real64 const > getGradNView() const
738  {
739  return m_viewGradN;
740  }
741 
746  arrayView2d< real64 const > getDetJView() const
747  {
748  return m_viewDetJ;
749  }
750 
751 
752 protected:
754  arrayView4d< real64 const > m_viewGradN;
755 
758  arrayView2d< real64 const > m_viewDetJ;
759 };
760 
762 
763 //*************************************************************************************************
764 //***** Definitions *******************************************************************************
765 //*************************************************************************************************
766 
767 template<>
768 struct FiniteElementBase::FunctionSpaceHelper< 1 >
769 {
770  GEOS_HOST_DEVICE
771  constexpr static PDEUtilities::FunctionSpace getFunctionSpace()
772  {
773  return PDEUtilities::FunctionSpace::H1;
774  }
775 };
776 
777 template<>
778 struct FiniteElementBase::FunctionSpaceHelper< 3 >
779 {
780  GEOS_HOST_DEVICE
781  constexpr static PDEUtilities::FunctionSpace getFunctionSpace()
782  {
783  return PDEUtilities::FunctionSpace::H1vector;
784  }
785 };
786 
787 template< int N >
788 GEOS_HOST_DEVICE
789 constexpr PDEUtilities::FunctionSpace FiniteElementBase::getFunctionSpace()
790 {
791  return FunctionSpaceHelper< N >::getFunctionSpace();
792 }
793 
794 template< typename LEAF >
795 GEOS_HOST_DEVICE
796 GEOS_FORCE_INLINE
797 real64 FiniteElementBase::getGradN( localIndex const k,
798  localIndex const q,
799  real64 const (&X)[LEAF::maxSupportPoints][3],
800  real64 (& gradN)[LEAF::maxSupportPoints][3] ) const
801 {
802  GEOS_UNUSED_VAR( k );
803  return LEAF::calcGradN( q, X, gradN );
804 }
805 
806 template< typename LEAF >
807 GEOS_HOST_DEVICE
808 GEOS_FORCE_INLINE
809 real64 FiniteElementBase::getGradN( localIndex const k,
810  localIndex const q,
811  real64 const (&X)[LEAF::maxSupportPoints][3],
812  typename LEAF::StackVariables const & stack,
813  real64 ( & gradN )[LEAF::maxSupportPoints][3] ) const
814 {
815  GEOS_UNUSED_VAR( k );
816  return LEAF::calcGradN( q, X, stack, gradN );
817 }
818 
819 template< typename LEAF >
820 GEOS_HOST_DEVICE
821 GEOS_FORCE_INLINE
822 real64 FiniteElementBase::getGradN( localIndex const k,
823  localIndex const q,
824  int const X,
825  real64 (& gradN)[LEAF::maxSupportPoints][3] ) const
826 {
827  GEOS_UNUSED_VAR( X );
828 
829  LvArray::tensorOps::copy< LEAF::maxSupportPoints, 3 >( gradN, m_viewGradN[ k ][ q ] );
830 
831  return m_viewDetJ( k, q );
832 }
833 
834 template< typename LEAF >
835 GEOS_HOST_DEVICE
836 GEOS_FORCE_INLINE
837 real64 FiniteElementBase::getGradN( localIndex const k,
838  localIndex const q,
839  int const X,
840  typename LEAF::StackVariables const & stack,
841  real64 (& gradN)[LEAF::maxSupportPoints][3] ) const
842 {
843  GEOS_UNUSED_VAR( X );
844  GEOS_UNUSED_VAR( stack );
845 
846  LvArray::tensorOps::copy< LEAF::maxSupportPoints, 3 >( gradN, m_viewGradN[ k ][ q ] );
847 
848  return m_viewDetJ( k, q );
849 }
850 
851 //*************************************************************************************************
852 //***** Interpolated Value Functions **************************************************************
853 //*************************************************************************************************
854 
855 template< int NUM_SUPPORT_POINTS >
856 GEOS_HOST_DEVICE
857 GEOS_FORCE_INLINE
858 void FiniteElementBase::value( real64 const (&N)[NUM_SUPPORT_POINTS],
859  real64 const (&var)[NUM_SUPPORT_POINTS],
860  real64 & value )
861 {
862  value = LvArray::tensorOps::AiBi< NUM_SUPPORT_POINTS >( N, var );
863 }
864 
865 template< int NUM_SUPPORT_POINTS,
866  int NUM_COMPONENTS >
867 GEOS_HOST_DEVICE
868 GEOS_FORCE_INLINE
869 void FiniteElementBase::value( real64 const (&N)[NUM_SUPPORT_POINTS],
870  real64 const (&var)[NUM_SUPPORT_POINTS][NUM_COMPONENTS],
871  real64 (& value)[NUM_COMPONENTS] )
872 {
873 
874  LvArray::tensorOps::Ri_eq_AjiBj< 3, NUM_SUPPORT_POINTS >( value, var, N );
875 }
876 
877 
878 //*************************************************************************************************
879 //***** Variable Gradient Functions ***************************************************************
880 //*************************************************************************************************
881 
882 template< int NUM_SUPPORT_POINTS,
883  typename GRADIENT_TYPE >
884 GEOS_HOST_DEVICE
885 GEOS_FORCE_INLINE
886 void FiniteElementBase::symmetricGradient( GRADIENT_TYPE const & gradN,
887  real64 const (&var)[NUM_SUPPORT_POINTS][3],
888  real64 (& gradVar)[6] )
889 {
890  gradVar[0] = gradN[0][0] * var[0][0];
891  gradVar[1] = gradN[0][1] * var[0][1];
892  gradVar[2] = gradN[0][2] * var[0][2];
893  gradVar[3] = gradN[0][2] * var[0][1] + gradN[0][1] * var[0][2];
894  gradVar[4] = gradN[0][2] * var[0][0] + gradN[0][0] * var[0][2];
895  gradVar[5] = gradN[0][1] * var[0][0] + gradN[0][0] * var[0][1];
896 
897  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
898  {
899  gradVar[0] = gradVar[0] + gradN[a][0] * var[ a ][0];
900  gradVar[1] = gradVar[1] + gradN[a][1] * var[ a ][1];
901  gradVar[2] = gradVar[2] + gradN[a][2] * var[ a ][2];
902  gradVar[3] = gradVar[3] + gradN[a][2] * var[ a ][1] + gradN[a][1] * var[ a ][2];
903  gradVar[4] = gradVar[4] + gradN[a][2] * var[ a ][0] + gradN[a][0] * var[ a ][2];
904  gradVar[5] = gradVar[5] + gradN[a][1] * var[ a ][0] + gradN[a][0] * var[ a ][1];
905  }
906 }
907 
908 template< int NUM_SUPPORT_POINTS,
909  typename GRADIENT_TYPE >
910 GEOS_HOST_DEVICE
911 GEOS_FORCE_INLINE
912 real64 FiniteElementBase::symmetricGradientTrace( GRADIENT_TYPE const & gradN,
913  real64 const (&var)[NUM_SUPPORT_POINTS][3] )
914 {
915  real64 result = gradN[0][0] * var[0][0] + gradN[0][1] * var[0][1] + gradN[0][2] * var[0][2];
916 
917  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
918  {
919  result = result + gradN[a][0] * var[a][0] + gradN[a][1] * var[a][1] + gradN[a][2] * var[a][2];
920  }
921  return result;
922 }
923 
924 template< int NUM_SUPPORT_POINTS,
925  typename GRADIENT_TYPE >
926 GEOS_HOST_DEVICE
927 GEOS_FORCE_INLINE
928 void FiniteElementBase::gradient( GRADIENT_TYPE const & gradN,
929  real64 const (&var)[NUM_SUPPORT_POINTS],
930  real64 (& gradVar)[3] )
931 {
932  LvArray::tensorOps::Ri_eq_AjiBj< 3, NUM_SUPPORT_POINTS >( gradVar, gradN, var );
933 }
934 
935 template< int NUM_SUPPORT_POINTS,
936  typename GRADIENT_TYPE >
937 GEOS_HOST_DEVICE
938 GEOS_FORCE_INLINE
939 void FiniteElementBase::gradient( GRADIENT_TYPE const & gradN,
940  real64 const (&var)[NUM_SUPPORT_POINTS][3],
941  real64 (& gradVar)[3][3] )
942 {
943  LvArray::tensorOps::Rij_eq_AkiBkj< 3, 3, NUM_SUPPORT_POINTS >( gradVar, var, gradN );
944 }
945 
946 
947 
948 template< int NUM_SUPPORT_POINTS,
949  typename GRADIENT_TYPE >
950 GEOS_HOST_DEVICE
951 GEOS_FORCE_INLINE
952 void FiniteElementBase::valueAndGradient( real64 const (&N)[NUM_SUPPORT_POINTS],
953  GRADIENT_TYPE const & gradN,
954  real64 const (&var)[NUM_SUPPORT_POINTS],
955  real64 & value,
956  real64 (& gradVar)[3] )
957 {
958  value = N[0] * var[0];
959  for( int i = 0; i < 3; ++i )
960  {
961  gradVar[i] = var[0] * gradN[0][i];
962  }
963 
964  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
965  {
966  value = value + N[a] * var[a];
967  for( int i = 0; i < 3; ++i )
968  {
969  gradVar[i] = gradVar[i] + var[ a ] * gradN[a][i];
970  }
971  }
972 }
973 
974 
975 
976 template< int NUM_SUPPORT_POINTS,
977  typename GRADIENT_TYPE >
978 GEOS_HOST_DEVICE
979 GEOS_FORCE_INLINE
980 void FiniteElementBase::plusGradNajAij( GRADIENT_TYPE const & gradN,
981  real64 const (&var_detJxW)[6],
982  real64 (& R)[NUM_SUPPORT_POINTS][3] )
983 {
984  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
985  {
986  R[a][0] = R[a][0] + var_detJxW[0] * gradN[a][0] + var_detJxW[5] * gradN[a][1] + var_detJxW[4] * gradN[a][2];
987  R[a][1] = R[a][1] + var_detJxW[5] * gradN[a][0] + var_detJxW[1] * gradN[a][1] + var_detJxW[3] * gradN[a][2];
988  R[a][2] = R[a][2] + var_detJxW[4] * gradN[a][0] + var_detJxW[3] * gradN[a][1] + var_detJxW[2] * gradN[a][2];
989  }
990 }
991 
992 
993 template< int NUM_SUPPORT_POINTS,
994  typename GRADIENT_TYPE >
995 GEOS_HOST_DEVICE
996 GEOS_FORCE_INLINE
997 void FiniteElementBase::plusGradNajAij( GRADIENT_TYPE const & gradN,
998  real64 const (&var_detJxW)[3][3],
999  real64 (& R)[NUM_SUPPORT_POINTS][3] )
1000 {
1001  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
1002  {
1003  LvArray::tensorOps::Ri_add_AijBj< 3, 3 >( R[a], var_detJxW, gradN[a] );
1004  }
1005 }
1006 
1007 template< int NUM_SUPPORT_POINTS >
1008 GEOS_HOST_DEVICE
1009 GEOS_FORCE_INLINE
1010 void FiniteElementBase::plusNaFi( real64 const (&N)[NUM_SUPPORT_POINTS],
1011  real64 const (&var_detJxW)[3],
1012  real64 ( & R )[NUM_SUPPORT_POINTS][3] )
1013 {
1014  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
1015  {
1016  LvArray::tensorOps::scaledAdd< 3 >( R[a], var_detJxW, N[a] );
1017  }
1018 }
1019 
1020 
1021 template< int NUM_SUPPORT_POINTS,
1022  typename GRADIENT_TYPE >
1023 GEOS_HOST_DEVICE
1024 GEOS_FORCE_INLINE
1025 void FiniteElementBase::plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
1026  real64 const (&var_detJxW)[6],
1027  real64 const (&N)[NUM_SUPPORT_POINTS],
1028  real64 const (&forcingTerm_detJxW)[3],
1029  real64 (& R)[NUM_SUPPORT_POINTS][3] )
1030 {
1031  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
1032  {
1033  R[a][0] = R[a][0] + var_detJxW[0] * gradN[a][0] + var_detJxW[5] * gradN[a][1] + var_detJxW[4] * gradN[a][2] + forcingTerm_detJxW[0] * N[a];
1034  R[a][1] = R[a][1] + var_detJxW[5] * gradN[a][0] + var_detJxW[1] * gradN[a][1] + var_detJxW[3] * gradN[a][2] + forcingTerm_detJxW[1] * N[a];
1035  R[a][2] = R[a][2] + var_detJxW[4] * gradN[a][0] + var_detJxW[3] * gradN[a][1] + var_detJxW[2] * gradN[a][2] + forcingTerm_detJxW[2] * N[a];
1036  }
1037 }
1038 
1039 template< int NUM_SUPPORT_POINTS,
1040  typename GRADIENT_TYPE >
1041 GEOS_HOST_DEVICE
1042 GEOS_FORCE_INLINE
1043 void FiniteElementBase::plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
1044  real64 const (&var_detJxW)[3][3],
1045  real64 const (&N)[NUM_SUPPORT_POINTS],
1046  real64 const (&forcingTerm_detJxW)[3],
1047  real64 (& R)[NUM_SUPPORT_POINTS][3] )
1048 {
1049  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
1050  {
1051  R[a][0] = R[a][0] + var_detJxW[0][0] * gradN[a][0] + var_detJxW[0][1] * gradN[a][1] + var_detJxW[0][2] * gradN[a][2] + forcingTerm_detJxW[0] * N[a];
1052  R[a][1] = R[a][1] + var_detJxW[1][0] * gradN[a][0] + var_detJxW[1][1] * gradN[a][1] + var_detJxW[1][2] * gradN[a][2] + forcingTerm_detJxW[1] * N[a];
1053  R[a][2] = R[a][2] + var_detJxW[2][0] * gradN[a][0] + var_detJxW[2][1] * gradN[a][1] + var_detJxW[2][2] * gradN[a][2] + forcingTerm_detJxW[2] * N[a];
1054  }
1055 }
1057 
1058 }
1059 }
1060 
1061 
1063 #define USING_FINITEELEMENTBASE \
1064  using FiniteElementBase::value; \
1065  using FiniteElementBase::symmetricGradient; \
1066  using FiniteElementBase::gradient; \
1067  using FiniteElementBase::valueAndGradient; \
1068  using FiniteElementBase::plusGradNajAij; \
1069  using FiniteElementBase::plusNaFi; \
1070  using FiniteElementBase::plusGradNajAijPlusNaFi;
1071 
1072 #endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_FINITEELEMENTBASE_HPP_
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_ERROR_IF_NE_MSG
#define GEOS_ERROR_IF_NE_MSG(lhs, rhs, msg)
Raise a hard error if two values are not equal.
Definition: Logger.hpp:243
geos::EdgeManager
This class provides an interface to ObjectManagerBase in order to manage edge data.
Definition: EdgeManager.hpp:43
geos::FaceManager
The FaceManager class provides an interface to ObjectManagerBase in order to manage face data.
Definition: FaceManager.hpp:44
geos::NodeManager
The NodeManager class provides an interface to ObjectManagerBase in order to manage node data.
Definition: NodeManager.hpp:46
geos::finiteElement::FiniteElementBase
Base class for FEM element implementations.
Definition: FiniteElementBase.hpp:46
geos::finiteElement::FiniteElementBase::value
static GEOS_HOST_DEVICE void value(real64 const (&N)[NUM_SUPPORT_POINTS], real64 const (&var)[NUM_SUPPORT_POINTS][NUM_COMPONENTS], real64(&value)[NUM_COMPONENTS])
Compute the interpolated value of a vector variable.
geos::finiteElement::FiniteElementBase::getNumSupportPoints
virtual GEOS_HOST_DEVICE localIndex getNumSupportPoints() const =0
Virtual getter for the number of support points per element.
geos::finiteElement::FiniteElementBase::plusGradNajAijPlusNaFi
static GEOS_HOST_DEVICE void plusGradNajAijPlusNaFi(GRADIENT_TYPE const &gradN, real64 const (&var_detJxW)[3][3], real64 const (&N)[NUM_SUPPORT_POINTS], real64 const (&forcingTerm_detJxW)[3], real64(&R)[NUM_SUPPORT_POINTS][3])
Inner product of each basis function gradient with a rank-2 symmetric tensor added to the product eac...
geos::finiteElement::FiniteElementBase::getMaxSupportPoints
virtual GEOS_HOST_DEVICE localIndex getMaxSupportPoints() const =0
Get the maximum number of support points for this element.
geos::finiteElement::FiniteElementBase::fillMeshData
static void fillMeshData(NodeManager const &nodeManager, EdgeManager const &edgeManager, FaceManager const &faceManager, SUBREGION_TYPE const &cellSubRegion, MeshData< SUBREGION_TYPE > &meshData)
Method to fill a MeshData object.
Definition: FiniteElementBase.hpp:123
geos::finiteElement::FiniteElementBase::getDetJView
arrayView2d< real64 const > getDetJView() const
Getter for m_viewDetJ.
Definition: FiniteElementBase.hpp:746
geos::finiteElement::FiniteElementBase::setupStack
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void setupStack(localIndex const &cellIndex, MeshData< SUBREGION_TYPE > const &meshData, StackVariables &stack)
Empty setup method.
Definition: FiniteElementBase.hpp:171
geos::finiteElement::FiniteElementBase::plusGradNajAijPlusNaFi
static GEOS_HOST_DEVICE void plusGradNajAijPlusNaFi(GRADIENT_TYPE const &gradN, real64 const (&var_detJxW)[6], real64 const (&N)[NUM_SUPPORT_POINTS], real64 const (&forcingTerm_detJxW)[3], real64(&R)[NUM_SUPPORT_POINTS][3])
Inner product of each basis function gradient with a rank-2 tensor added to the product each shape fu...
geos::finiteElement::FiniteElementBase::~FiniteElementBase
virtual GEOS_HOST_DEVICE ~FiniteElementBase()
Destructor.
Definition: FiniteElementBase.hpp:92
geos::finiteElement::FiniteElementBase::gradient
static GEOS_HOST_DEVICE void gradient(GRADIENT_TYPE const &gradN, real64 const (&var)[NUM_SUPPORT_POINTS], real64(&gradVar)[3])
Calculate the gradient of a scalar valued support field at a point using the input basis function gra...
geos::finiteElement::FiniteElementBase::plusGradNajAij
static GEOS_HOST_DEVICE void plusGradNajAij(GRADIENT_TYPE const &gradN, real64 const (&var_detJxW)[3][3], real64(&R)[NUM_SUPPORT_POINTS][3])
Inner product of each basis function gradient with a rank-2 symmetric tensor.
geos::finiteElement::FiniteElementBase::operator=
FiniteElementBase & operator=(FiniteElementBase const &)=delete
Deleted copy assignment operator.
geos::finiteElement::FiniteElementBase::setup
GEOS_HOST_DEVICE void setup(localIndex const &cellIndex, typename LEAF::template MeshData< SUBREGION_TYPE > const &meshData, typename LEAF::StackVariables &stack) const
Abstract setup method, possibly computing cell-dependent properties.
Definition: FiniteElementBase.hpp:191
geos::finiteElement::FiniteElementBase::m_viewGradN
arrayView4d< real64 const > m_viewGradN
View to potentially hold pre-calculated shape function gradients.
Definition: FiniteElementBase.hpp:754
geos::finiteElement::FiniteElementBase::FiniteElementBase
GEOS_HOST_DEVICE FiniteElementBase(FiniteElementBase const &source)
Copy Constructor.
Definition: FiniteElementBase.hpp:61
geos::finiteElement::FiniteElementBase::addGradGradStabilization
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void addGradGradStabilization(StackVariables const &stack, real64(&matrix)[MAXSUPPORTPOINTS *NUMDOFSPERTRIALSUPPORTPOINT][MAXSUPPORTPOINTS *NUMDOFSPERTRIALSUPPORTPOINT], real64 const &scaleFactor)
Empty method, here for compatibility with methods that require a stabilization of the grad-grad bilin...
Definition: FiniteElementBase.hpp:342
geos::finiteElement::FiniteElementBase::getGradN
GEOS_HOST_DEVICE real64 getGradN(localIndex const k, localIndex const q, real64 const (&X)[LEAF::maxSupportPoints][3], typename LEAF::StackVariables const &stack, real64(&gradN)[LEAF::maxSupportPoints][3]) const
Get the shape function gradients.
geos::finiteElement::FiniteElementBase::setDetJView
void setDetJView(arrayView2d< real64 const > const &source)
Sets m_viewDetJ equal to an input view.
Definition: FiniteElementBase.hpp:725
geos::finiteElement::FiniteElementBase::plusGradNajAij
static GEOS_HOST_DEVICE void plusGradNajAij(GRADIENT_TYPE const &gradN, real64 const (&var_detJxW)[6], real64(&R)[NUM_SUPPORT_POINTS][3])
Inner product of each basis function gradient with a rank-2 symmetric tensor.
geos::finiteElement::FiniteElementBase::FiniteElementBase
FiniteElementBase()=default
Default Constructor.
geos::finiteElement::FiniteElementBase::getGradN
GEOS_HOST_DEVICE real64 getGradN(localIndex const k, localIndex const q, int const X, real64(&gradN)[LEAF::maxSupportPoints][3]) const
Get the shape function gradients.
geos::finiteElement::FiniteElementBase::numSupportPoints
GEOS_HOST_DEVICE localIndex numSupportPoints(typename LEAF::StackVariables const &stack) const
Getter for the number of support points per element.
Definition: FiniteElementBase.hpp:239
geos::finiteElement::FiniteElementBase::symmetricGradient
static GEOS_HOST_DEVICE void symmetricGradient(GRADIENT_TYPE const &gradN, real64 const (&var)[NUM_SUPPORT_POINTS][3], real64(&gradVar)[6])
Calculate the symmetric gradient of a vector valued support field at a point using the stored basis f...
geos::finiteElement::FiniteElementBase::getNumQuadraturePoints
virtual GEOS_HOST_DEVICE localIndex getNumQuadraturePoints() const =0
Virtual getter for the number of quadrature points per element.
geos::finiteElement::FiniteElementBase::initialize
static void initialize(NodeManager const &nodeManager, EdgeManager const &edgeManager, FaceManager const &faceManager, SUBREGION_TYPE const &cellSubRegion, typename LEAF::template MeshData< SUBREGION_TYPE > &meshData)
Abstract initialization method.
Definition: FiniteElementBase.hpp:147
geos::finiteElement::FiniteElementBase::addEvaluatedGradGradStabilization
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void addEvaluatedGradGradStabilization(StackVariables const &stack, real64 const (&dofs)[MAXSUPPORTPOINTS][NUMDOFSPERTRIALSUPPORTPOINT], real64(&targetVector)[MAXSUPPORTPOINTS][NUMDOFSPERTRIALSUPPORTPOINT], real64 const scaleFactor)
Empty method, here for compatibility with methods that require a stabilization of the grad-grad bilin...
Definition: FiniteElementBase.hpp:394
geos::finiteElement::FiniteElementBase::operator=
FiniteElementBase & operator=(FiniteElementBase &&)=delete
Deleted move assignment operator.
geos::finiteElement::FiniteElementBase::symmetricGradientTrace
static GEOS_HOST_DEVICE real64 symmetricGradientTrace(GRADIENT_TYPE const &gradN, real64 const (&var)[NUM_SUPPORT_POINTS][3])
Calculate the trace of the symmetric gradient of a vector valued support field (i....
geos::finiteElement::FiniteElementBase::addGradGradStabilizationMatrix
GEOS_HOST_DEVICE void addGradGradStabilizationMatrix(typename LEAF::StackVariables const &stack, real64(&matrix)[LEAF::maxSupportPoints *NUMDOFSPERTRIALSUPPORTPOINT][LEAF::maxSupportPoints *NUMDOFSPERTRIALSUPPORTPOINT], real64 const scaleFactor=1.0) const
Add stabilization of grad-grad bilinear form to input matrix.
Definition: FiniteElementBase.hpp:365
geos::finiteElement::FiniteElementBase::plusNaFi
static GEOS_HOST_DEVICE void plusNaFi(real64 const (&N)[NUM_SUPPORT_POINTS], real64 const (&forcingTerm_detJxW)[3], real64(&R)[NUM_SUPPORT_POINTS][3])
Product of each shape function with a vector forcing term.
geos::finiteElement::FiniteElementBase::getFunctionSpace
constexpr static GEOS_HOST_DEVICE PDEUtilities::FunctionSpace getFunctionSpace()
Getter for the function space.
geos::finiteElement::FiniteElementBase::addEvaluatedGradGradStabilizationVector
GEOS_HOST_DEVICE GEOS_FORCE_INLINE void addEvaluatedGradGradStabilizationVector(typename LEAF::StackVariables const &stack, real64 const (&dofs)[LEAF::maxSupportPoints][NUMDOFSPERTRIALSUPPORTPOINT], real64(&targetVector)[LEAF::maxSupportPoints][NUMDOFSPERTRIALSUPPORTPOINT], real64 const scaleFactor=1.0) const
Add a grad-grad stabilization operator evaluated at a provided vector of dofs to input vector.
Definition: FiniteElementBase.hpp:422
geos::finiteElement::FiniteElementBase::getGradN
GEOS_HOST_DEVICE real64 getGradN(localIndex const k, localIndex const q, real64 const (&X)[LEAF::maxSupportPoints][3], real64(&gradN)[LEAF::maxSupportPoints][3]) const
Get the shape function gradients.
geos::finiteElement::FiniteElementBase::value
static GEOS_HOST_DEVICE void value(real64 const (&N)[NUM_SUPPORT_POINTS], real64 const (&var)[NUM_SUPPORT_POINTS], real64 &value)
Compute the interpolated value of a variable.
geos::finiteElement::FiniteElementBase::numSamplingPointsPerDirection
constexpr static int numSamplingPointsPerDirection
Number of sampling points.
Definition: FiniteElementBase.hpp:54
geos::finiteElement::FiniteElementBase::getGradNView
arrayView4d< real64 const > getGradNView() const
Getter for m_viewGradN.
Definition: FiniteElementBase.hpp:737
geos::finiteElement::FiniteElementBase::FiniteElementBase
FiniteElementBase(FiniteElementBase &&)=default
Default Move constructor.
geos::finiteElement::FiniteElementBase::gradient
static GEOS_HOST_DEVICE void gradient(GRADIENT_TYPE const &gradN, real64 const (&var)[NUM_SUPPORT_POINTS][3], real64(&gradVar)[3][3])
Calculate the gradient of a vector valued support field at a point using the input basis function gra...
geos::finiteElement::FiniteElementBase::getGradN
GEOS_HOST_DEVICE real64 getGradN(localIndex const k, localIndex const q, int const X, typename LEAF::StackVariables const &stack, real64(&gradN)[LEAF::maxSupportPoints][3]) const
Get the shape function gradients.
geos::finiteElement::FiniteElementBase::m_viewDetJ
arrayView2d< real64 const > m_viewDetJ
Definition: FiniteElementBase.hpp:758
geos::finiteElement::FiniteElementBase::setGradNView
void setGradNView(arrayView4d< real64 const > const &source)
Sets m_viewGradN equal to an input view.
Definition: FiniteElementBase.hpp:706
geos::finiteElement::FiniteElementBase::valueAndGradient
static GEOS_HOST_DEVICE void valueAndGradient(real64 const (&N)[NUM_SUPPORT_POINTS], GRADIENT_TYPE const &gradN, real64 const (&var)[NUM_SUPPORT_POINTS], real64 &value, real64(&gradVar)[3])
Calculate the value and gradient of a scalar valued support field at a point using the input basis fu...
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
geos::arrayView4d
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
Definition: DataTypes.hpp:228
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:196
geos::finiteElement::FiniteElementBase::FunctionSpaceHelper
An helper struct to determine the function space.
Definition: FiniteElementBase.hpp:219
geos::finiteElement::FiniteElementBase::MeshData
Variables used to initialize the class.
Definition: FiniteElementBase.hpp:111
Generated by doxygen 1.9.1