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 Total, S.A
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  {
107  GEOS_HOST_DEVICE
108  StackVariables()
109  {}
110  };
111 
116  template< typename SUBREGION_TYPE >
117  struct MeshData
118  {
122  MeshData()
123  {}
124  };
125 
126 
135  template< typename SUBREGION_TYPE >
136  static void fillMeshData( NodeManager const & nodeManager,
137  EdgeManager const & edgeManager,
138  FaceManager const & faceManager,
139  SUBREGION_TYPE const & cellSubRegion,
140  MeshData< SUBREGION_TYPE > & meshData )
141  {
142  GEOS_UNUSED_VAR( nodeManager,
143  edgeManager,
144  faceManager,
145  cellSubRegion,
146  meshData );
147  }
148 
159  template< typename LEAF, typename SUBREGION_TYPE >
160  static void initialize( NodeManager const & nodeManager,
161  EdgeManager const & edgeManager,
162  FaceManager const & faceManager,
163  SUBREGION_TYPE const & cellSubRegion,
164  typename LEAF::template MeshData< SUBREGION_TYPE > & meshData
165  )
166  {
167  LEAF::template fillMeshData< SUBREGION_TYPE >( nodeManager,
168  edgeManager,
169  faceManager,
170  cellSubRegion,
171  meshData );
172  }
173 
174 
181  template< typename SUBREGION_TYPE >
182  GEOS_HOST_DEVICE
183  GEOS_FORCE_INLINE
184  static void setupStack( localIndex const & cellIndex,
185  MeshData< SUBREGION_TYPE > const & meshData,
186  StackVariables & stack )
187  {
188  GEOS_UNUSED_VAR( cellIndex,
189  meshData,
190  stack );
191 
192  }
193 
202  template< typename LEAF, typename SUBREGION_TYPE >
203  GEOS_HOST_DEVICE
204  void setup( localIndex const & cellIndex,
205  typename LEAF::template MeshData< SUBREGION_TYPE > const & meshData,
206  typename LEAF::StackVariables & stack ) const
207  {
208  LEAF::setupStack( cellIndex, meshData, stack );
209  }
210 
215  GEOS_HOST_DEVICE
216  virtual localIndex getNumQuadraturePoints() const = 0;
217 
222  GEOS_HOST_DEVICE
223  virtual localIndex getNumSupportPoints() const = 0;
224 
230  template< int N >
231  struct FunctionSpaceHelper
232  {};
233 
240  template< int N >
241  GEOS_HOST_DEVICE
242  constexpr static PDEUtilities::FunctionSpace getFunctionSpace();
243 
250  template< typename LEAF >
251  GEOS_HOST_DEVICE
252  localIndex numSupportPoints( typename LEAF::StackVariables const & stack ) const
253  {
254  return LEAF::getNumSupportPoints( stack );
255  }
256 
263  GEOS_HOST_DEVICE
264  virtual localIndex getMaxSupportPoints() const = 0;
265 
277  template< typename LEAF >
278  GEOS_HOST_DEVICE
279  real64 getGradN( localIndex const k,
280  localIndex const q,
281  real64 const (&X)[LEAF::maxSupportPoints][3],
282  real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
283 
296  template< typename LEAF >
297  GEOS_HOST_DEVICE
298  real64 getGradN( localIndex const k,
299  localIndex const q,
300  real64 const (&X)[LEAF::maxSupportPoints][3],
301  typename LEAF::StackVariables const & stack,
302  real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
303 
315  template< typename LEAF >
316  GEOS_HOST_DEVICE
317  real64 getGradN( localIndex const k,
318  localIndex const q,
319  int const X,
320  real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
333  template< typename LEAF >
334  GEOS_HOST_DEVICE
335  real64 getGradN( localIndex const k,
336  localIndex const q,
337  int const X,
338  typename LEAF::StackVariables const & stack,
339  real64 ( &gradN )[LEAF::maxSupportPoints][3] ) const;
340 
341 
352  template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, localIndex MAXSUPPORTPOINTS, bool UPPER >
353  GEOS_HOST_DEVICE
354  GEOS_FORCE_INLINE
355  static void addGradGradStabilization( StackVariables const & stack,
356  real64 ( & matrix )
357  [MAXSUPPORTPOINTS * NUMDOFSPERTRIALSUPPORTPOINT]
358  [MAXSUPPORTPOINTS * NUMDOFSPERTRIALSUPPORTPOINT],
359  real64 const & scaleFactor )
360  {
361  GEOS_UNUSED_VAR( stack,
362  matrix,
363  scaleFactor );
364  }
365 
366 
376  template< typename LEAF, localIndex NUMDOFSPERTRIALSUPPORTPOINT, bool UPPER = false >
377  GEOS_HOST_DEVICE
378  void addGradGradStabilizationMatrix( typename LEAF::StackVariables const & stack,
379  real64 ( & matrix )
380  [LEAF::maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT]
381  [LEAF::maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT],
382  real64 const scaleFactor = 1.0 ) const
383  {
384  LEAF::template addGradGradStabilization< NUMDOFSPERTRIALSUPPORTPOINT,
385  LEAF::maxSupportPoints,
386  UPPER >( stack,
387  matrix,
388  scaleFactor );
389  }
390 
404  template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, localIndex MAXSUPPORTPOINTS >
405  GEOS_HOST_DEVICE
406  GEOS_FORCE_INLINE
407  static void addEvaluatedGradGradStabilization( StackVariables const & stack,
408  real64 const ( &dofs )[MAXSUPPORTPOINTS][NUMDOFSPERTRIALSUPPORTPOINT],
409  real64 ( & targetVector )[MAXSUPPORTPOINTS][NUMDOFSPERTRIALSUPPORTPOINT],
410  real64 const scaleFactor )
411  {
412  GEOS_UNUSED_VAR( stack );
413  GEOS_UNUSED_VAR( dofs );
414  GEOS_UNUSED_VAR( targetVector );
415  GEOS_UNUSED_VAR( scaleFactor );
416  }
417 
431  template< typename LEAF, localIndex NUMDOFSPERTRIALSUPPORTPOINT >
432  GEOS_HOST_DEVICE
433  GEOS_FORCE_INLINE
434  void
435  addEvaluatedGradGradStabilizationVector( typename LEAF::StackVariables const & stack,
436  real64 const ( &dofs )[LEAF::maxSupportPoints]
437  [NUMDOFSPERTRIALSUPPORTPOINT],
438  real64 ( & targetVector )[LEAF::maxSupportPoints]
439  [NUMDOFSPERTRIALSUPPORTPOINT],
440  real64 const scaleFactor = 1.0 ) const
441  {
442  LEAF::template
443  addEvaluatedGradGradStabilization< NUMDOFSPERTRIALSUPPORTPOINT, LEAF::maxSupportPoints >( stack,
444  dofs,
445  targetVector,
446  scaleFactor );
447  }
448 
453 
473  template< int NUM_SUPPORT_POINTS >
474  GEOS_HOST_DEVICE
475  static
476  void value( real64 const (&N)[NUM_SUPPORT_POINTS],
477  real64 const (&var)[NUM_SUPPORT_POINTS],
478  real64 & value );
479 
485  template< int NUM_SUPPORT_POINTS,
486  int NUM_COMPONENTS >
487  GEOS_HOST_DEVICE
488  static
489  void value( real64 const (&N)[NUM_SUPPORT_POINTS],
490  real64 const (&var)[NUM_SUPPORT_POINTS][NUM_COMPONENTS],
491  real64 ( &value )[NUM_COMPONENTS] );
492 
494 
499 
518  template< int NUM_SUPPORT_POINTS,
519  typename GRADIENT_TYPE >
520  GEOS_HOST_DEVICE
521  static void symmetricGradient( GRADIENT_TYPE const & gradN,
522  real64 const (&var)[NUM_SUPPORT_POINTS][3],
523  real64 ( &gradVar )[6] );
524 
537  template< int NUM_SUPPORT_POINTS,
538  typename GRADIENT_TYPE >
539  GEOS_HOST_DEVICE
540  static real64 symmetricGradientTrace( GRADIENT_TYPE const & gradN,
541  real64 const (&var)[NUM_SUPPORT_POINTS][3] );
542 
559  template< int NUM_SUPPORT_POINTS,
560  typename GRADIENT_TYPE >
561  GEOS_HOST_DEVICE
562  static void gradient( GRADIENT_TYPE const & gradN,
563  real64 const (&var)[NUM_SUPPORT_POINTS],
564  real64 ( &gradVar )[3] );
565 
576  template< int NUM_SUPPORT_POINTS,
577  typename GRADIENT_TYPE >
578  GEOS_HOST_DEVICE
579  static void gradient( GRADIENT_TYPE const & gradN,
580  real64 const (&var)[NUM_SUPPORT_POINTS][3],
581  real64 ( &gradVar )[3][3] );
583 
588 
602  template< int NUM_SUPPORT_POINTS,
603  typename GRADIENT_TYPE >
604  GEOS_HOST_DEVICE
605  static void valueAndGradient( real64 const (&N)[NUM_SUPPORT_POINTS],
606  GRADIENT_TYPE const & gradN,
607  real64 const (&var)[NUM_SUPPORT_POINTS],
608  real64 & value,
609  real64 ( &gradVar )[3] );
610 
612 
620 
641  template< int NUM_SUPPORT_POINTS,
642  typename GRADIENT_TYPE >
643  GEOS_HOST_DEVICE
644  static void plusGradNajAij( GRADIENT_TYPE const & gradN,
645  real64 const (&var_detJxW)[6],
646  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
647 
653  template< int NUM_SUPPORT_POINTS,
654  typename GRADIENT_TYPE >
655  GEOS_HOST_DEVICE
656  static void plusGradNajAij( GRADIENT_TYPE const & gradN,
657  real64 const (&var_detJxW)[3][3],
658  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
659 
668  template< int NUM_SUPPORT_POINTS >
669  GEOS_HOST_DEVICE
670  static void plusNaFi( real64 const (&N)[NUM_SUPPORT_POINTS],
671  real64 const (&forcingTerm_detJxW)[3],
672  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
673 
691  template< int NUM_SUPPORT_POINTS,
692  typename GRADIENT_TYPE >
693  GEOS_HOST_DEVICE
694  static void plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
695  real64 const (&var_detJxW)[3][3],
696  real64 const (&N)[NUM_SUPPORT_POINTS],
697  real64 const (&forcingTerm_detJxW)[3],
698  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
699 
705  template< int NUM_SUPPORT_POINTS,
706  typename GRADIENT_TYPE >
707  GEOS_HOST_DEVICE
708  static void plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
709  real64 const (&var_detJxW)[6],
710  real64 const (&N)[NUM_SUPPORT_POINTS],
711  real64 const (&forcingTerm_detJxW)[3],
712  real64 ( &R )[NUM_SUPPORT_POINTS][3] );
713 
714 
719  void setGradNView( arrayView4d< real64 const > const & source )
720  {
721  GEOS_ERROR_IF_NE_MSG( source.size( 1 ),
722  getNumQuadraturePoints(),
723  "2nd-dimension of gradN array does not match number of quadrature points" );
724  GEOS_ERROR_IF_NE_MSG( source.size( 2 ),
725  getMaxSupportPoints(),
726  "3rd-dimension of gradN array does not match number of support points" );
727  GEOS_ERROR_IF_NE_MSG( source.size( 3 ),
728  3,
729  "4th-dimension of gradN array does not match 3" );
730 
731  m_viewGradN = source;
732  }
733 
738  void setDetJView( arrayView2d< real64 const > const & source )
739  {
740  GEOS_ERROR_IF_NE_MSG( source.size( 1 ),
741  getNumQuadraturePoints(),
742  "2nd-dimension of gradN array does not match number of quadrature points" );
743  m_viewDetJ = source;
744  }
745 
750  arrayView4d< real64 const > getGradNView() const
751  {
752  return m_viewGradN;
753  }
754 
759  arrayView2d< real64 const > getDetJView() const
760  {
761  return m_viewDetJ;
762  }
763 
764 
765 protected:
767  arrayView4d< real64 const > m_viewGradN;
768 
771  arrayView2d< real64 const > m_viewDetJ;
772 };
773 
775 
776 //*************************************************************************************************
777 //***** Definitions *******************************************************************************
778 //*************************************************************************************************
779 
780 template<>
781 struct FiniteElementBase::FunctionSpaceHelper< 1 >
782 {
783  GEOS_HOST_DEVICE
784  constexpr static PDEUtilities::FunctionSpace getFunctionSpace()
785  {
786  return PDEUtilities::FunctionSpace::H1;
787  }
788 };
789 
790 template<>
791 struct FiniteElementBase::FunctionSpaceHelper< 3 >
792 {
793  GEOS_HOST_DEVICE
794  constexpr static PDEUtilities::FunctionSpace getFunctionSpace()
795  {
796  return PDEUtilities::FunctionSpace::H1vector;
797  }
798 };
799 
800 template< int N >
801 GEOS_HOST_DEVICE
802 constexpr PDEUtilities::FunctionSpace FiniteElementBase::getFunctionSpace()
803 {
804  return FunctionSpaceHelper< N >::getFunctionSpace();
805 }
806 
807 template< typename LEAF >
808 GEOS_HOST_DEVICE
809 GEOS_FORCE_INLINE
810 real64 FiniteElementBase::getGradN( localIndex const k,
811  localIndex const q,
812  real64 const (&X)[LEAF::maxSupportPoints][3],
813  real64 (& gradN)[LEAF::maxSupportPoints][3] ) const
814 {
815  GEOS_UNUSED_VAR( k );
816  return LEAF::calcGradN( q, X, 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  real64 const (&X)[LEAF::maxSupportPoints][3],
825  typename LEAF::StackVariables const & stack,
826  real64 ( & gradN )[LEAF::maxSupportPoints][3] ) const
827 {
828  GEOS_UNUSED_VAR( k );
829  return LEAF::calcGradN( q, X, stack, gradN );
830 }
831 
832 template< typename LEAF >
833 GEOS_HOST_DEVICE
834 GEOS_FORCE_INLINE
835 real64 FiniteElementBase::getGradN( localIndex const k,
836  localIndex const q,
837  int const X,
838  real64 (& gradN)[LEAF::maxSupportPoints][3] ) const
839 {
840  GEOS_UNUSED_VAR( X );
841 
842  LvArray::tensorOps::copy< LEAF::maxSupportPoints, 3 >( gradN, m_viewGradN[ k ][ q ] );
843 
844  return m_viewDetJ( k, q );
845 }
846 
847 template< typename LEAF >
848 GEOS_HOST_DEVICE
849 GEOS_FORCE_INLINE
850 real64 FiniteElementBase::getGradN( localIndex const k,
851  localIndex const q,
852  int const X,
853  typename LEAF::StackVariables const & stack,
854  real64 (& gradN)[LEAF::maxSupportPoints][3] ) const
855 {
856  GEOS_UNUSED_VAR( X );
857  GEOS_UNUSED_VAR( stack );
858 
859  LvArray::tensorOps::copy< LEAF::maxSupportPoints, 3 >( gradN, m_viewGradN[ k ][ q ] );
860 
861  return m_viewDetJ( k, q );
862 }
863 
864 //*************************************************************************************************
865 //***** Interpolated Value Functions **************************************************************
866 //*************************************************************************************************
867 
868 template< int NUM_SUPPORT_POINTS >
869 GEOS_HOST_DEVICE
870 GEOS_FORCE_INLINE
871 void FiniteElementBase::value( real64 const (&N)[NUM_SUPPORT_POINTS],
872  real64 const (&var)[NUM_SUPPORT_POINTS],
873  real64 & value )
874 {
875  value = LvArray::tensorOps::AiBi< NUM_SUPPORT_POINTS >( N, var );
876 }
877 
878 template< int NUM_SUPPORT_POINTS,
879  int NUM_COMPONENTS >
880 GEOS_HOST_DEVICE
881 GEOS_FORCE_INLINE
882 void FiniteElementBase::value( real64 const (&N)[NUM_SUPPORT_POINTS],
883  real64 const (&var)[NUM_SUPPORT_POINTS][NUM_COMPONENTS],
884  real64 (& value)[NUM_COMPONENTS] )
885 {
886 
887  LvArray::tensorOps::Ri_eq_AjiBj< 3, NUM_SUPPORT_POINTS >( value, var, N );
888 }
889 
890 
891 //*************************************************************************************************
892 //***** Variable Gradient Functions ***************************************************************
893 //*************************************************************************************************
894 
895 template< int NUM_SUPPORT_POINTS,
896  typename GRADIENT_TYPE >
897 GEOS_HOST_DEVICE
898 GEOS_FORCE_INLINE
899 void FiniteElementBase::symmetricGradient( GRADIENT_TYPE const & gradN,
900  real64 const (&var)[NUM_SUPPORT_POINTS][3],
901  real64 (& gradVar)[6] )
902 {
903  gradVar[0] = gradN[0][0] * var[0][0];
904  gradVar[1] = gradN[0][1] * var[0][1];
905  gradVar[2] = gradN[0][2] * var[0][2];
906  gradVar[3] = gradN[0][2] * var[0][1] + gradN[0][1] * var[0][2];
907  gradVar[4] = gradN[0][2] * var[0][0] + gradN[0][0] * var[0][2];
908  gradVar[5] = gradN[0][1] * var[0][0] + gradN[0][0] * var[0][1];
909 
910  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
911  {
912  gradVar[0] = gradVar[0] + gradN[a][0] * var[ a ][0];
913  gradVar[1] = gradVar[1] + gradN[a][1] * var[ a ][1];
914  gradVar[2] = gradVar[2] + gradN[a][2] * var[ a ][2];
915  gradVar[3] = gradVar[3] + gradN[a][2] * var[ a ][1] + gradN[a][1] * var[ a ][2];
916  gradVar[4] = gradVar[4] + gradN[a][2] * var[ a ][0] + gradN[a][0] * var[ a ][2];
917  gradVar[5] = gradVar[5] + gradN[a][1] * var[ a ][0] + gradN[a][0] * var[ a ][1];
918  }
919 }
920 
921 template< int NUM_SUPPORT_POINTS,
922  typename GRADIENT_TYPE >
923 GEOS_HOST_DEVICE
924 GEOS_FORCE_INLINE
925 real64 FiniteElementBase::symmetricGradientTrace( GRADIENT_TYPE const & gradN,
926  real64 const (&var)[NUM_SUPPORT_POINTS][3] )
927 {
928  real64 result = gradN[0][0] * var[0][0] + gradN[0][1] * var[0][1] + gradN[0][2] * var[0][2];
929 
930  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
931  {
932  result = result + gradN[a][0] * var[a][0] + gradN[a][1] * var[a][1] + gradN[a][2] * var[a][2];
933  }
934  return result;
935 }
936 
937 template< int NUM_SUPPORT_POINTS,
938  typename GRADIENT_TYPE >
939 GEOS_HOST_DEVICE
940 GEOS_FORCE_INLINE
941 void FiniteElementBase::gradient( GRADIENT_TYPE const & gradN,
942  real64 const (&var)[NUM_SUPPORT_POINTS],
943  real64 (& gradVar)[3] )
944 {
945  LvArray::tensorOps::Ri_eq_AjiBj< 3, NUM_SUPPORT_POINTS >( gradVar, gradN, var );
946 }
947 
948 template< int NUM_SUPPORT_POINTS,
949  typename GRADIENT_TYPE >
950 GEOS_HOST_DEVICE
951 GEOS_FORCE_INLINE
952 void FiniteElementBase::gradient( GRADIENT_TYPE const & gradN,
953  real64 const (&var)[NUM_SUPPORT_POINTS][3],
954  real64 (& gradVar)[3][3] )
955 {
956  LvArray::tensorOps::Rij_eq_AkiBkj< 3, 3, NUM_SUPPORT_POINTS >( gradVar, var, gradN );
957 }
958 
959 
960 
961 template< int NUM_SUPPORT_POINTS,
962  typename GRADIENT_TYPE >
963 GEOS_HOST_DEVICE
964 GEOS_FORCE_INLINE
965 void FiniteElementBase::valueAndGradient( real64 const (&N)[NUM_SUPPORT_POINTS],
966  GRADIENT_TYPE const & gradN,
967  real64 const (&var)[NUM_SUPPORT_POINTS],
968  real64 & value,
969  real64 (& gradVar)[3] )
970 {
971  value = N[0] * var[0];
972  for( int i = 0; i < 3; ++i )
973  {
974  gradVar[i] = var[0] * gradN[0][i];
975  }
976 
977  for( int a=1; a<NUM_SUPPORT_POINTS; ++a )
978  {
979  value = value + N[a] * var[a];
980  for( int i = 0; i < 3; ++i )
981  {
982  gradVar[i] = gradVar[i] + var[ a ] * gradN[a][i];
983  }
984  }
985 }
986 
987 
988 
989 template< int NUM_SUPPORT_POINTS,
990  typename GRADIENT_TYPE >
991 GEOS_HOST_DEVICE
992 GEOS_FORCE_INLINE
993 void FiniteElementBase::plusGradNajAij( GRADIENT_TYPE const & gradN,
994  real64 const (&var_detJxW)[6],
995  real64 (& R)[NUM_SUPPORT_POINTS][3] )
996 {
997  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
998  {
999  R[a][0] = R[a][0] + var_detJxW[0] * gradN[a][0] + var_detJxW[5] * gradN[a][1] + var_detJxW[4] * gradN[a][2];
1000  R[a][1] = R[a][1] + var_detJxW[5] * gradN[a][0] + var_detJxW[1] * gradN[a][1] + var_detJxW[3] * gradN[a][2];
1001  R[a][2] = R[a][2] + var_detJxW[4] * gradN[a][0] + var_detJxW[3] * gradN[a][1] + var_detJxW[2] * gradN[a][2];
1002  }
1003 }
1004 
1005 
1006 template< int NUM_SUPPORT_POINTS,
1007  typename GRADIENT_TYPE >
1008 GEOS_HOST_DEVICE
1009 GEOS_FORCE_INLINE
1010 void FiniteElementBase::plusGradNajAij( GRADIENT_TYPE const & gradN,
1011  real64 const (&var_detJxW)[3][3],
1012  real64 (& R)[NUM_SUPPORT_POINTS][3] )
1013 {
1014  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
1015  {
1016  LvArray::tensorOps::Ri_add_AijBj< 3, 3 >( R[a], var_detJxW, gradN[a] );
1017  }
1018 }
1019 
1020 template< int NUM_SUPPORT_POINTS >
1021 GEOS_HOST_DEVICE
1022 GEOS_FORCE_INLINE
1023 void FiniteElementBase::plusNaFi( real64 const (&N)[NUM_SUPPORT_POINTS],
1024  real64 const (&var_detJxW)[3],
1025  real64 ( & R )[NUM_SUPPORT_POINTS][3] )
1026 {
1027  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
1028  {
1029  LvArray::tensorOps::scaledAdd< 3 >( R[a], var_detJxW, N[a] );
1030  }
1031 }
1032 
1033 
1034 template< int NUM_SUPPORT_POINTS,
1035  typename GRADIENT_TYPE >
1036 GEOS_HOST_DEVICE
1037 GEOS_FORCE_INLINE
1038 void FiniteElementBase::plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
1039  real64 const (&var_detJxW)[6],
1040  real64 const (&N)[NUM_SUPPORT_POINTS],
1041  real64 const (&forcingTerm_detJxW)[3],
1042  real64 (& R)[NUM_SUPPORT_POINTS][3] )
1043 {
1044  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
1045  {
1046  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];
1047  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];
1048  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];
1049  }
1050 }
1051 
1052 template< int NUM_SUPPORT_POINTS,
1053  typename GRADIENT_TYPE >
1054 GEOS_HOST_DEVICE
1055 GEOS_FORCE_INLINE
1056 void FiniteElementBase::plusGradNajAijPlusNaFi( GRADIENT_TYPE const & gradN,
1057  real64 const (&var_detJxW)[3][3],
1058  real64 const (&N)[NUM_SUPPORT_POINTS],
1059  real64 const (&forcingTerm_detJxW)[3],
1060  real64 (& R)[NUM_SUPPORT_POINTS][3] )
1061 {
1062  for( int a=0; a<NUM_SUPPORT_POINTS; ++a )
1063  {
1064  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];
1065  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];
1066  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];
1067  }
1068 }
1070 
1071 }
1072 }
1073 
1074 
1076 #define USING_FINITEELEMENTBASE \
1077  using FiniteElementBase::value; \
1078  using FiniteElementBase::symmetricGradient; \
1079  using FiniteElementBase::gradient; \
1080  using FiniteElementBase::valueAndGradient; \
1081  using FiniteElementBase::plusGradNajAij; \
1082  using FiniteElementBase::plusNaFi; \
1083  using FiniteElementBase::plusGradNajAijPlusNaFi;
1084 
1085 #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:136
geos::finiteElement::FiniteElementBase::getDetJView
arrayView2d< real64 const > getDetJView() const
Getter for m_viewDetJ.
Definition: FiniteElementBase.hpp:759
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:184
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:204
geos::finiteElement::FiniteElementBase::m_viewGradN
arrayView4d< real64 const > m_viewGradN
View to potentially hold pre-calculated shape function gradients.
Definition: FiniteElementBase.hpp:767
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:355
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:738
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:252
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:160
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:407
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:378
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:435
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:750
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:771
geos::finiteElement::FiniteElementBase::setGradNView
void setGradNView(arrayView4d< real64 const > const &source)
Sets m_viewGradN equal to an input view.
Definition: FiniteElementBase.hpp:719
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:232
geos::finiteElement::FiniteElementBase::MeshData
Variables used to initialize the class.
Definition: FiniteElementBase.hpp:118
Generated by doxygen 1.9.1