GEOS
BB_Tetrahedron.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) 2018-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_BBTETRAHEDRON_HPP_
21 #define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_BBTETRAHEDRON_HPP_
22 
23 #include "FiniteElementBase.hpp"
24 #include <utility>
25 
26 
27 
28 namespace geos
29 {
30 namespace finiteElement
31 {
32 
33 namespace
34 {
35 template< int ORDER >
36 constexpr int BB_Tetrahedron_NumNodes = ( ORDER + 1 ) * ( ORDER + 2 ) * ( ORDER + 3 ) / 6;
37 }
38 
50 template< int ORDER >
51 class BB_Tetrahedron_impl : public FiniteElementBase_impl< BB_Tetrahedron_NumNodes< ORDER >,
52  4,
53  BB_Tetrahedron_NumNodes< ORDER > >
54 {
55 public:
57  using Base = FiniteElementBase_impl< BB_Tetrahedron_NumNodes< ORDER >,
58  4,
59  BB_Tetrahedron_NumNodes< ORDER > >;
60 
62  using StackVariables = typename Base::StackVariables;
63 
65  template< typename SubregionType >
66  using MeshData = typename Base::template MeshData< SubregionType >;
67 
69  using Base::numNodes;
70 
72  using Base::numQuadraturePoints;
73 
75  using Base::maxSupportPoints;
76 
78  static constexpr int order = ORDER;
79 
81  constexpr static localIndex numNodesPerFace = ( ORDER + 1 ) * ( ORDER + 2 ) / 2;
82 
84  GEOS_HOST_DEVICE
85  static localIndex getNumQuadraturePoints()
86  {
87  return numQuadraturePoints;
88  }
89 
95  GEOS_HOST_DEVICE
96  GEOS_FORCE_INLINE
97  static localIndex getNumQuadraturePoints( StackVariables const & stack )
98  {
99  GEOS_UNUSED_VAR( stack );
100  return numQuadraturePoints;
101  }
102 
104  GEOS_HOST_DEVICE
105  GEOS_FORCE_INLINE
106  static localIndex getNumSupportPoints()
107  {
108  return numNodes;
109  }
110 
112  GEOS_HOST_DEVICE
113  GEOS_FORCE_INLINE
114  static localIndex getMaxSupportPoints()
115  {
116  return maxSupportPoints;
117  }
118 
124  GEOS_HOST_DEVICE
125  GEOS_FORCE_INLINE
126  static localIndex getNumSupportPoints( StackVariables const & stack )
127  {
128  GEOS_UNUSED_VAR( stack );
129  return numNodes;
130  }
131 
132 
138  GEOS_HOST_DEVICE
139  GEOS_FORCE_INLINE
140  static real64 jacobianDeterminant( real64 const (&X)[4][3] )
141  {
142  real64 m[3][3] = {};
143  for( int i = 0; i < 3; i++ )
144  {
145  for( int j = 0; j < 3; j++ )
146  {
147  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
148  }
149  }
150  return LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) );
151  }
152 
159  GEOS_HOST_DEVICE
160  GEOS_FORCE_INLINE
161  static real64 faceJacobianDeterminant( localIndex face,
162  real64 const (&X)[4][3] )
163  {
164  int i1 = ( face + 1 ) % 4;
165  int i2 = ( face + 2 ) % 4;
166  int i3 = ( face + 3 ) % 4;
167 
168  real64 ab[3] = { X[ i2 ][ 0 ] - X[ i1 ][ 0 ],
169  X[ i2 ][ 1 ] - X[ i1 ][ 1 ],
170  X[ i2 ][ 2 ] - X[ i1 ][ 2 ] };
171  real64 ac[3] = { X[ i3 ][ 0 ] - X[ i1 ][ 0 ],
172  X[ i3 ][ 1 ] - X[ i1 ][ 1 ],
173  X[ i3 ][ 2 ] - X[ i1 ][ 2 ] };
174  real64 term1 = ab[1] * ac[2] - ab[2] * ac[1];
175  real64 term2 = ab[2] * ac[0] - ab[0] * ac[2];
176  real64 term3 = ab[0] * ac[1] - ab[1] * ac[0];
177  return LvArray::math::sqrt( ( term1 * term1 + term2 * term2 + term3 * term3 ) );
178 
179  }
180 
187  GEOS_HOST_DEVICE
188  GEOS_FORCE_INLINE
189  static void calcN( localIndex const,
190  real64 (& N)[numNodes] )
191  {
192  for( int a=0; a < numNodes; ++a )
193  {
194  N[ a ] = 0;
195  }
196  GEOS_ERROR( "Bernstein-Bézier basis is modal, not nodal. No quadrature points are defined." );
197  }
198 
207  GEOS_HOST_DEVICE
208  GEOS_FORCE_INLINE
209  static void calcN( localIndex const q,
210  StackVariables const & stack,
211  real64 ( & N )[numNodes] )
212  {
213  GEOS_UNUSED_VAR( stack );
214  return calcN( q, N );
215  }
221  GEOS_HOST_DEVICE
222  GEOS_FORCE_INLINE
223  static void calcN( real64 const ( &lambda)[4],
224  real64 (& N)[numNodes] )
225  {
226  N[ 0 ] = 6.0;
227  int prev;
228  int c;
229  int limits[ 4 ] = { 1, 1, 1, 1 };
230  for( int np = 1; np <= ORDER; np++ )
231  {
232  prev = np * ( np + 1 ) * ( np + 2 ) / 6 - 1;
233  c = ( np + 1 ) * ( np + 2 ) * ( np + 3 ) / 6 - 1;
234  for( int i = 0; i < 4; i++ )
235  {
236  int denominator = i == 0 ? np : 1;
237  int offset = 0;
238  int c1 = np - 1;
239  int c2 = i + np - 2;
240  int repetitionCount = i == 0 ? 1 : limits[ i - 1 ];
241  for( int j = 0; j < limits[ i ]; j++ )
242  {
243  if( j - offset >= repetitionCount )
244  {
245  denominator++;
246  offset += repetitionCount;
247  repetitionCount = repetitionCount * c1 / c2;
248  c1--;
249  c2--;
250  }
251  N[ c-- ] = N[ prev - j ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
252  }
253  }
254  for( int i = 1; i < 4; i++ )
255  {
256  limits[ i ] += limits[ i - 1 ];
257  }
258  }
259  }
260 
269  GEOS_HOST_DEVICE
270  GEOS_FORCE_INLINE
271  static void calcN( real64 const (&X)[4][3],
272  real64 const (&coords)[3],
273  real64 (& N)[numNodes] )
274  {
275  real64 lambda[4] = {};
276  real64 m[3][3] = {};
277  for( int i = 0; i < 3; i++ )
278  {
279  for( int j = 0; j < 3; j++ )
280  {
281  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
282  }
283  }
284  real64 den = LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) );
285  for( int i = 0; i < 3; i++ )
286  {
287  for( int j = 0; j < 3; j++ )
288  {
289  m[ i ][ j ] = coords[ j ] - X[ 0 ][ j ];
290  }
291  lambda[ i + 1 ] = LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) ) / den;
292  for( int j = 0; j < 3; j++ )
293  {
294  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
295  }
296  }
297  lambda[ 0 ] = 1.0 - lambda[ 1 ] - lambda[ 2 ] - lambda[ 3 ];
298  return calcN( lambda, N );
299  }
300 
308  GEOS_HOST_DEVICE
309  GEOS_FORCE_INLINE
310  static void calcNandGradN( real64 const ( &lambda)[4],
311  real64 const ( &N)[numNodes],
312  real64 (& gradN)[numNodes][ 4 ] )
313  {
314  gradN[ 0 ][ 0 ] = 0.0;
315  gradN[ 0 ][ 1 ] = 0.0;
316  gradN[ 0 ][ 2 ] = 0.0;
317  gradN[ 0 ][ 3 ] = 0.0;
318  N[ 0 ] = 6.0;
319  int prev;
320  int c;
321  int limits[ 4 ] = { 1, 1, 1, 1 };
322  for( int np = 1; np <= ORDER; np++ )
323  {
324  prev = np * ( np + 1 ) * ( np + 2 ) / 6 - 1;
325  c = ( np + 1 ) * ( np + 2 ) * ( np + 3 ) / 6 - 1;
326  for( int i = 0; i < 4; i++ )
327  {
328  int denominator = i == 0 ? np : 1;
329  int offset = 0;
330  int c1 = np - 1;
331  int c2 = i + np - 2;
332  int repetitionCount = i == 0 ? 1 : limits[ i - 1 ];
333  for( int j = 0; j < limits[ i ]; j++ )
334  {
335  if( j - offset >= repetitionCount )
336  {
337  denominator++;
338  offset += repetitionCount;
339  repetitionCount = repetitionCount * c1 / c2;
340  c1--;
341  c2--;
342  }
343  gradN[ c ][ 0 ] = gradN[ prev - j ][ 0 ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
344  gradN[ c ][ 1 ] = gradN[ prev - j ][ 1 ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
345  gradN[ c ][ 2 ] = gradN[ prev - j ][ 2 ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
346  gradN[ c ][ 3 ] = gradN[ prev - j ][ 3 ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
347  gradN[ c ][ 3 - i ] += N[ prev - j ] * ( np + 3 ) / denominator;
348  N[ c-- ] = N[ prev - j ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
349  }
350  }
351  for( int i = 1; i < 4; i++ )
352  {
353  limits[ i ] += limits[ i - 1 ];
354  }
355  }
356  }
357 
358 
367  GEOS_HOST_DEVICE
368  GEOS_FORCE_INLINE
369  static void calcNandGradN( real64 const (&X)[4][3],
370  real64 const (&coords)[3],
371  real64 (& N)[numNodes],
372  real64 (& gradN)[numNodes][ 4 ] )
373  {
374  real64 lambda[4] = {};
375  real64 m[3][3] = {};
376  real64 dNdLambda[numNodes][4];
377  for( int i = 0; i < 3; i++ )
378  {
379  for( int j = 0; j < 3; j++ )
380  {
381  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
382  }
383  }
384  real64 den = LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) );
385  for( int i = 0; i < 3; i++ )
386  {
387  for( int j = 0; j < 3; j++ )
388  {
389  m[ i ][ j ] = coords[ j ] - X[ 0 ][ j ];
390  }
391  lambda[ i + 1 ] = LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) ) / den;
392  for( int j = 0; j < 3; j++ )
393  {
394  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
395  }
396  }
397  lambda[ 0 ] = 1.0 - lambda[ 1 ] - lambda[ 2 ] - lambda[ 3 ];
398  calcNandGradN( lambda, N, dNdLambda );
399  for( int i = 0; i < numNodes; i++ )
400  {
401  for( int j = 0; j < 3; j++ )
402  {
403  gradN[ i ][ j ] =
404  ( ( ( X[ 2 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 3 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) - ( X[ 3 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 2 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) )*
405  ( dNdLambda[ i ][ 1 ] - dNdLambda[ i ][ 0 ] ) +
406  ( ( X[ 3 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 1 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) - ( X[ 1 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) *
407  ( X[ 3 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) )* ( dNdLambda[ i ][ 2 ] - dNdLambda[ i ][ 0 ] ) +
408  ( ( X[ 1 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 2 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) - ( X[ 2 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) *
409  ( X[ 1 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) )* ( dNdLambda[ i ][ 3 ] - dNdLambda[ i ][ 0 ] )) / den;
410  }
411  }
412  }
413 
421  GEOS_HOST_DEVICE
422  GEOS_FORCE_INLINE
423  static real64 calcJacobian( localIndex const q,
424  real64 const (&X)[numNodes][3],
425  real64 (& J)[3][3] )
426  {
427  GEOS_UNUSED_VAR( q );
428  for( int i = 0; i < 3; ++i )
429  {
430  for( int j = 0; j < 3; ++j )
431  {
432  J[i][j] = 0;
433  }
434  }
435  real64 m[3][3] = {};
436  for( int i = 0; i < 3; i++ )
437  {
438  for( int j = 0; j < 3; j++ )
439  {
440  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
441  }
442  }
443  real64 const detJ = LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) );
444  return detJ;
445  }
446 
455  GEOS_HOST_DEVICE
456  GEOS_FORCE_INLINE
457  static real64 calcJacobian( localIndex const q,
458  real64 const (&X)[numNodes][3],
459  StackVariables const & stack,
460  real64 (& J)[3][3] )
461  {
462  GEOS_UNUSED_VAR( stack );
463  return calcJacobian( q, X, J );
464  }
465 
475  GEOS_HOST_DEVICE
476  GEOS_FORCE_INLINE
477  static real64 calcGradN( localIndex const q,
478  real64 const (&X)[numNodes][3],
479  real64 ( & gradN )[numNodes][3] )
480  {
481  gradN[ q ][ 0 ] = 0.0;
482  gradN[ q ][ 1 ] = 0.0;
483  gradN[ q ][ 2 ] = 0.0;
484  GEOS_UNUSED_VAR( q, X, gradN );
485  GEOS_ERROR( "Bernstein-Bézier basis is modal, not nodal. No quadrature points are defined." );
486  return 0;
487  }
498  GEOS_HOST_DEVICE
499  GEOS_FORCE_INLINE
500  static real64 calcGradN( localIndex const q,
501  real64 const (&X)[numNodes][3],
502  StackVariables const & stack,
503  real64 ( & gradN )[numNodes][3] )
504  {
505  GEOS_UNUSED_VAR( stack );
506  return calcGradN( q, X, gradN );
507  }
508 
516  GEOS_HOST_DEVICE
517  GEOS_FORCE_INLINE
518  constexpr static real64 integralTerm( const int a, const int b, const int c )
519  {
520  real64 res = 1.0;
521  int num = a;
522  int den = c;
523  for( int i = b; i > 0; i-- )
524  {
525  res *= ( (real64) num ) / i;
526  num--;
527  }
528  for( int i = num; i > 0; i-- )
529  {
530  res *= ( (real64) i ) / den;
531  den--;
532  }
533  for( int i = den; i > 0; i-- )
534  {
535  res /= i;
536  }
537 
538  return res;
539  }
540 
555  GEOS_HOST_DEVICE
556  GEOS_FORCE_INLINE
557  static constexpr real64 correctionFactorDerivative( int const i1, int const j1, int const k1, int const l1, int const i2, int const j2, int const k2, int const l2, int const dim )
558  {
559  return (i1+j1+k1+l1+dim)* (i2==0 ? 1 : i2) * (j2==0 ? 1 : j2)* (k2==0 ? 1 : k2)* (l2==0 ? 1 : l2);
560  }
561 
562 
576  GEOS_HOST_DEVICE
577  GEOS_FORCE_INLINE
578  static constexpr real64 computeSuperpositionIntegral( const int i1, const int j1, const int k1, const int l1,
579  const int i2, const int j2, const int k2, const int l2 )
580  {
581  return (integralTerm( i1+i2, i1, i2 )*
582  integralTerm( j1+j2, j1, j2 )*
583  integralTerm( k1+k2, k1, k2 )*
584  integralTerm( l1+l2, l1, l2 ))/
585  integralTerm( i1+j1+k1+l1+i2+j2+k2+l2+3,
586  i1+j1+k1+l1+3, i2+j2+k2+l2+3 );
587  }
588 
600  GEOS_HOST_DEVICE
601  GEOS_FORCE_INLINE
602  constexpr static real64 computeFaceSuperpositionIntegral( const int i1, const int j1, const int k1,
603  const int i2, const int j2, const int k2 )
604  {
605  return ((i1+k1+j1+3)*(i2+j2+k2+3)*
606  integralTerm( i1+i2, i1, i2 )*
607  integralTerm( j1+j2, j1, j2 )*
608  integralTerm( k1+k2, k1, k2 ))/
609  integralTerm( i1+j1+k1+i2+j2+k2+2,
610  i1+j1+k1+2, i2+j2+k2+2 );
611  }
612 
621  template< int I, int J, int K >
622  GEOS_HOST_DEVICE
623  GEOS_FORCE_INLINE
624  static
625  constexpr
626  int
627  dofIndex()
628  {
629  return ( ORDER - I ) * ( ORDER - I + 1 ) * ( ORDER - I + 2) / 6 +
630  ( ORDER - I - J ) * ( ORDER - I - J + 1 ) / 2 +
631  ( ORDER - I - J - K );
632  }
633 
643  GEOS_HOST_DEVICE
644  GEOS_FORCE_INLINE
645  static
646  constexpr
647  int
648  dofIndex( const int i, const int j, const int k )
649  {
650  return ( ORDER - i ) * ( ORDER - i + 1 ) * ( ORDER - i + 2) / 6 +
651  ( ORDER - i - j ) * ( ORDER - i - j + 1 ) / 2 +
652  ( ORDER - i - j - k );
653  }
654 
655 
656 
663  template< int C, int VTX >
664  GEOS_HOST_DEVICE
665  GEOS_FORCE_INLINE
666  static
667  constexpr
668  int
669  indexToIJKL()
670  {
671  static_assert( VTX >= 0 && VTX < 4 );
672  // compute the indices of c in the current element using tetrahedral and triangular roots
673  constexpr int cc1 = C + 1;
674  constexpr real64 tetr = cbrt( 3.0 * cc1 + sqrt( 9.0 * cc1 * cc1 - 1.0 / 27.0 ) )
675  + cbrt( 3.0 * cc1 - sqrt( 9.0 * cc1 * cc1 - 1.0 / 27.0 ) ) - 2;
676  constexpr int i = ORDER - round( tetr * 10.0 ) / 10;
677  constexpr int cc2 = C - ( ORDER - i ) * ( ORDER - i + 1 ) * ( ORDER - i + 2 ) / 6 + 1;
678  constexpr real64 trir = ( sqrt( 8.0 * cc2 + 1.0 ) - 1.0 ) / 2.0 - 1;
679  constexpr int j = ORDER - i - round( trir * 10.0 ) / 10;
680  constexpr int k = ORDER - i - j - ( C - (ORDER - i ) * ( ORDER - i + 1 ) * ( ORDER - i + 2 ) / 6
681  - ( ORDER - i - j ) * ( ORDER - i - j + 1 ) / 2 );
682  if constexpr ( VTX == 0 )
683  {
684  return i;
685  }
686  else if constexpr ( VTX == 1)
687  {
688  return j;
689  }
690  else if constexpr ( VTX == 2)
691  {
692  return k;
693  }
694  else if constexpr ( VTX == 3)
695  {
696  return ORDER - i - j - k;
697  }
698  return -1;
699  }
700 
708  template< typename FUNC, int... Is >
709  GEOS_HOST_DEVICE
710  GEOS_FORCE_INLINE
711  static constexpr void loop_impl( FUNC && func, std::integer_sequence< int, Is... > )
712  {
713  ( func( std::integral_constant< int, Is >{} ), ... );
714  }
715 
723  template< int N, typename FUNC >
724  GEOS_HOST_DEVICE
725  static constexpr void loop( FUNC const & func )
726  {
727  if constexpr (N > 0)
728  {
729  loop< N - 1 >( func );
730  func( std::integral_constant< int, N - 1 >{} );
731  }
732  }
733 
734 
740  template< typename FUNC >
741  GEOS_HOST_DEVICE
742  GEOS_FORCE_INLINE
743  static constexpr void barycentricCoordinateLoop( FUNC && func )
744  {
745  loop_impl( std::forward< FUNC >( func ), std::make_integer_sequence< int, 4 >{} );
746  }
747 
753  template< typename FUNC >
754  GEOS_HOST_DEVICE
755  GEOS_FORCE_INLINE
756  static constexpr void basisLoop( FUNC && func )
757  {
758  loop< ORDER + 1 >( [&func] ( auto const iic )
759  {
760  constexpr int i = decltype(iic)::value;
761  constexpr int i1 = ORDER - i;
762  loop< ORDER + 1 >( [&func] ( auto const jjc )
763  {
764  constexpr int j = decltype(jjc)::value;
765  constexpr int j1 = ORDER - j;
766  if constexpr ( j1 <= ORDER - i1 )
767  {
768  loop< ORDER + 1 >( [&func] ( auto const kkc )
769  {
770  constexpr int k = decltype(kkc)::value;
771  constexpr int k1 = ORDER - k;
772  if constexpr ( k1 <= (ORDER - i1 - j1) )
773  {
774  constexpr int l1 = ORDER - i1 - j1 - k1;
775  constexpr int c1 = dofIndex< i1, j1, k1 >();
776  func( std::integral_constant< int, c1 >{},
777  std::integral_constant< int, i1 >{},
778  std::integral_constant< int, j1 >{},
779  std::integral_constant< int, k1 >{},
780  std::integral_constant< int, l1 >{} );
781  }
782  } );
783  }
784  } );
785  } );
786  }
787 
801  template< int c1, int i1, int j1, int k1, int l1, typename F, int... Is >
802  GEOS_HOST_DEVICE
803  static constexpr void call_matching_cases( F && func, std::integer_sequence< int, Is... > )
804  {
805 
806  auto check_i1 = [&]( auto I ) {
807  if constexpr (i1 == decltype(I)::value) {
808  func( std::integral_constant< int, 0 >{},
809  std::integral_constant< int, i1 >{},
810  std::integral_constant< int, c1 >{},
811  std::integral_constant< int, j1 >{},
812  std::integral_constant< int, k1 >{},
813  std::integral_constant< int, l1 >{} );
814  }
815  };
816  (check_i1( std::integral_constant< int, Is >{} ), ...);
817 
818  auto check_j1 = [&]( auto I ) {
819  if constexpr (j1 == decltype(I)::value) {
820  func( std::integral_constant< int, 1 >{},
821  std::integral_constant< int, j1 >{},
822  std::integral_constant< int, c1 >{},
823  std::integral_constant< int, i1 >{},
824  std::integral_constant< int, k1 >{},
825  std::integral_constant< int, l1 >{} );
826  }
827  };
828  (check_j1( std::integral_constant< int, Is >{} ), ...);
829 
830  auto check_k1 = [&]( auto I ) {
831  if constexpr (k1 == decltype(I)::value) {
832  func( std::integral_constant< int, 2 >{},
833  std::integral_constant< int, k1 >{},
834  std::integral_constant< int, c1 >{},
835  std::integral_constant< int, i1 >{},
836  std::integral_constant< int, j1 >{},
837  std::integral_constant< int, l1 >{} );
838  }
839  };
840  (check_k1( std::integral_constant< int, Is >{} ), ...);
841 
842  auto check_l1 = [&]( auto I ) {
843  if constexpr (l1 == decltype(I)::value) {
844  func( std::integral_constant< int, 3 >{},
845  std::integral_constant< int, l1 >{},
846  std::integral_constant< int, c1 >{},
847  std::integral_constant< int, i1 >{},
848  std::integral_constant< int, j1 >{},
849  std::integral_constant< int, k1 >{} );
850  }
851  };
852  (check_l1( std::integral_constant< int, Is >{} ), ...);
853  }
854 
855 
863  template< int... Is, typename FUNC >
864  GEOS_HOST_DEVICE
865  static constexpr void conditionalBasisLoop( FUNC const & func )
866  {
867  loop< ORDER + 1 >( [&]( auto const i ) {
868  using i_t = decltype(i);
869  constexpr int iVal = i_t::value;
870  constexpr int i1 = ORDER - iVal;
871 
872  loop< ORDER + 1 >( [&]( auto const j ) {
873  using j_t = decltype(j);
874  constexpr int jVal = j_t::value;
875  constexpr int j1 = ORDER - jVal;
876 
877  constexpr bool valid_j1_i1 = (j1 <= iVal);
878  if constexpr (valid_j1_i1) {
879 
880  loop< ORDER + 1 >( [&]( auto const k ) {
881  using k_t = decltype(k);
882  constexpr int kVal = k_t::value;
883  constexpr int k1 = ORDER - kVal;
884 
885  constexpr bool valid_k1 = (k1 <= (ORDER - i1 - j1));
886  if constexpr (valid_k1) {
887  constexpr int l1 = ORDER - i1 - j1 - k1;
888  constexpr int c1 = dofIndex< i1, j1, k1 >();
889 
890  call_matching_cases< c1, i1, j1, k1, l1 >( func, std::integer_sequence< int, Is... >{} );
891 
892  // Pour chaque Is...
893  // (void)(((i1 == Is) &&
894  // (void(func(std::integral_constant<int, 0>{},
895  // std::integral_constant<int, i1>{},
896  // std::integral_constant<int, c1>{},
897  // std::integral_constant<int, j1>{},
898  // std::integral_constant<int, k1>{},
899  // std::integral_constant<int, l1>{})), 1)) || ...);
900 
901  // (void)(((j1 == Is) &&
902  // (void(func(std::integral_constant<int, 1>{},
903  // std::integral_constant<int, j1>{},
904  // std::integral_constant<int, c1>{},
905  // std::integral_constant<int, i1>{},
906  // std::integral_constant<int, k1>{},
907  // std::integral_constant<int, l1>{})), 1)) || ...);
908 
909  // (void)(((k1==Is) &&
910  // (void(func(std::integral_constant<int, 2>{},
911  // std::integral_constant<int, k1>{},
912  // std::integral_constant<int, c1>{},
913  // std::integral_constant<int, i1>{},
914  // std::integral_constant<int, j1>{},
915  // std::integral_constant<int, l1>{})), 1)) || ...);
916 
917  // (void)(((l1 ==Is) &&
918  // (void(func(std::integral_constant<int, 3>{},
919  // std::integral_constant<int, l1>{},
920  // std::integral_constant<int, c1>{},
921  // std::integral_constant<int, i1>{},
922  // std::integral_constant<int, j1>{},
923  // std::integral_constant<int, k1>{})), 1)) || ...);
924 
925 
926 
927  // (void)(((i1 ==Is) &&
928  // (void(func(
929  // HDInt<0>{},
930  // HDInt<i1>{},
931  // HDInt<c1>{},
932  // HDInt<j1>{},
933  // HDInt<k1>{},
934  // HDInt<l1>{})), 1)) || ...);
935 
936  // (void)(((j1==Is) &&
937  // (void(func(HDInt<1>{},
938  // HDInt<j1>{},
939  // HDInt<c1>{},
940  // HDInt<i1>{},
941  // HDInt<k1>{},
942  // HDInt<l1>{})), 1)) || ...);
943 
944  // (void)(((k1==Is) &&
945  // (void(func(
946  // HDInt<2>{},
947  // HDInt<k1>{},
948  // HDInt<c1>{},
949  // HDInt<i1>{},
950  // HDInt<j1>{},
951  // HDInt<l1>{})), 1)) || ...);
952 
953  // (void)(((l1==Is) &&
954  // (void(func(
955  // HDInt<3>{},
956  // HDInt<l1>{},
957  // HDInt<c1>{},
958  // HDInt<i1>{},
959  // HDInt<j1>{},
960  // HDInt<k1>{})), 1)) || ...);
961  // i1 matching
962  // callIfMatch<i1, Is...>([&] {
963  // func(HDInt<0>{}, HDInt<i1>{}, HDInt<c1>{}, HDInt<j1>{}, HDInt<k1>{}, HDInt<l1>{});
964  // });
965 
966  // callIfMatch<j1, Is...>([&] {
967  // func(HDInt<1>{}, HDInt<j1>{}, HDInt<c1>{}, HDInt<i1>{}, HDInt<k1>{}, HDInt<l1>{});
968  // });
969 
970  // callIfMatch<k1, Is...>([&] {
971  // func(HDInt<2>{}, HDInt<k1>{}, HDInt<c1>{}, HDInt<i1>{}, HDInt<j1>{}, HDInt<l1>{});
972  // });
973 
974  // callIfMatch<l1, Is...>([&] {
975  // func(HDInt<3>{}, HDInt<l1>{}, HDInt<c1>{}, HDInt<i1>{}, HDInt<j1>{}, HDInt<k1>{});
976  // });
977  // callIfMatch<i1, Is...>([&] {
978  // constexpr auto h0 = HDInt<0>{};
979  // constexpr auto hi1 = HDInt<i1>{};
980  // constexpr auto hc1 = HDInt<c1>{};
981  // constexpr auto hj1 = HDInt<j1>{};
982  // constexpr auto hk1 = HDInt<k1>{};
983  // constexpr auto hl1 = HDInt<l1>{};
984 
985  // func(h0, hi1, hc1, hj1, hk1, hl1);
986  // });
987 
988  // callIfMatch<j1, Is...>([&] {
989  // constexpr auto h1 = HDInt<1>{};
990  // constexpr auto hj1 = HDInt<j1>{};
991  // constexpr auto hc1 = HDInt<c1>{};
992  // constexpr auto hi1 = HDInt<i1>{};
993  // constexpr auto hk1 = HDInt<k1>{};
994  // constexpr auto hl1 = HDInt<l1>{};
995 
996  // func(h1, hj1, hc1, hi1, hk1, hl1);
997  // });
998 
999  // callIfMatch<k1, Is...>([&] {
1000  // constexpr auto h2 = HDInt<2>{};
1001  // constexpr auto hk1 = HDInt<k1>{};
1002  // constexpr auto hc1 = HDInt<c1>{};
1003  // constexpr auto hi1 = HDInt<i1>{};
1004  // constexpr auto hj1 = HDInt<j1>{};
1005  // constexpr auto hl1 = HDInt<l1>{};
1006 
1007  // func(h2, hk1, hc1, hi1, hj1, hl1);
1008  // });
1009 
1010  // callIfMatch<l1, Is...>([&] {
1011  // constexpr auto h3 = HDInt<3>{};
1012  // constexpr auto hl1 = HDInt<l1>{};
1013  // constexpr auto hc1 = HDInt<c1>{};
1014  // constexpr auto hi1 = HDInt<i1>{};
1015  // constexpr auto hj1 = HDInt<j1>{};
1016  // constexpr auto hk1 = HDInt<k1>{};
1017 
1018  // func(h3, hl1, hc1, hi1, hj1, hk1);
1019  // });
1020 
1021  }
1022  } );
1023  }
1024  } );
1025  } );
1026  }
1032  template< typename FUNC >
1033  GEOS_HOST_DEVICE
1034  GEOS_FORCE_INLINE
1035  static constexpr void faceBarycentricCoordinateLoop( FUNC && func )
1036  {
1037  loop< 3 >( std::forward< FUNC >( func ) );
1038  }
1039 
1047  GEOS_HOST_DEVICE
1048  GEOS_FORCE_INLINE
1049  static
1050  void
1051  computeReferenceMassMatrix( arraySlice2d< real64 > const & m )
1052  {
1053  basisLoop( [ &m ] ( auto const cc1, auto const ii1, auto const jj1, auto const kk1, auto const ll1 )
1054  {
1055  constexpr int c1 = decltype(cc1)::value;
1056  constexpr int i1 = decltype(ii1)::value;
1057  constexpr int j1 = decltype(jj1)::value;
1058  constexpr int k1 = decltype(kk1)::value;
1059  constexpr int l1 = decltype(ll1)::value;
1060  // Needed for compilers that do not support constexpr lambdas
1061  GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
1062  basisLoop( [ &m ] ( auto const cc2, auto const ii2, auto const jj2, auto const kk2, auto const ll2 )
1063  {
1064  constexpr int c2 = decltype(cc2)::value;
1065  constexpr int i2 = decltype(ii2)::value;
1066  constexpr int j2 = decltype(jj2)::value;
1067  constexpr int k2 = decltype(kk2)::value;
1068  constexpr int l2 = decltype(ll2)::value;
1069  constexpr real64 val = computeSuperpositionIntegral( i1, j1, k1, l1, i2, j2, k2, l2 );
1070  m[ c1 ][ c2 ] = val;
1071  } );
1072  } );
1073  }
1074 
1086  template< typename DAMPING >
1087  GEOS_HOST_DEVICE
1088  GEOS_FORCE_INLINE
1089  static
1090  void
1091  computeReferenceDampingMatrix( real64 (& d)[numNodes][numNodes], bool const face1Damped, bool const face2Damped, bool const face3Damped, bool const face4Damped )
1092  {
1093 
1094  for( int c1 = 0; c1 < numNodes; c1++ )
1095  {
1096  for( int c2 = 0; c2 < numNodes; c2++ )
1097  {
1098  d[ c1 ][ c2 ] = 0;
1099  }
1100  }
1101  conditionalBasisLoop< 0 >( [&] ( auto const f1, auto const, auto const c1, auto const i1, auto const j1, auto const k1 )
1102  {
1103  conditionalBasisLoop< 0 >( [&] ( auto const f2, auto const, auto const c2, auto const i2, auto const j2, auto const k2 )
1104  {
1105  if constexpr ( f1 == f2 )
1106  {
1107  constexpr real64 val = computeFaceSuperpositionIntegral( i1, j1, k1, i2, j2, k2 );
1108  if( ( f1 == 0 && face1Damped ) ||
1109  ( f1 == 1 && face2Damped ) ||
1110  ( f1 == 2 && face3Damped ) ||
1111  ( f1 == 3 && face4Damped ) )
1112  {
1113  d[ c1 ][ c2 ] += val;
1114  }
1115  }
1116  } );
1117  } );
1118  }
1119 
1126  template< typename FUNC >
1127  GEOS_HOST_DEVICE
1128  GEOS_FORCE_INLINE
1129  static
1130  constexpr
1131  void
1132  computeMassTerm( real64 const (&X)[4][3],
1133  FUNC && func )
1134  {
1135  real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1136  basisLoop( [&func, &detJ] ( auto cc1, auto ii1, auto jj1, auto kk1, auto ll1 )
1137  {
1138  constexpr int c1 = decltype(cc1)::value;
1139  constexpr int i1 = decltype(ii1)::value;
1140  constexpr int j1 = decltype(jj1)::value;
1141  constexpr int k1 = decltype(kk1)::value;
1142  constexpr int l1 = decltype(ll1)::value;
1143  GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
1144 
1145  basisLoop( [&func, &detJ, i1, j1, k1, l1] ( auto c2, auto ii2, auto jj2, auto kk2, auto ll2 )
1146  {
1147  constexpr int c2v = decltype(c2)::value;
1148  constexpr int i2 = decltype(ii2)::value;
1149  constexpr int j2 = decltype(jj2)::value;
1150  constexpr int k2 = decltype(kk2)::value;
1151  constexpr int l2 = decltype(ll2)::value;
1152 
1153  constexpr real64 val = computeSuperpositionIntegral( i1, j1, k1, l1, i2, j2, k2, l2 );
1154  func( c1, c2v, val * detJ );
1155  } );
1156  } );
1157  }
1158 
1168  GEOS_HOST_DEVICE
1169  GEOS_FORCE_INLINE
1170  static constexpr real64
1171  computeFluxDerivativeFactor( real64 const (&X)[4][3], int x1, int x2, int o1, int o2 ) // Order x1, x2, o1, o2
1172  {
1173 
1174  real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1175 
1176  // height with respect to o2
1177 
1178  real64 detJf1 = LvArray::math::abs( faceJacobianDeterminant( o1, X ));
1179 
1180  real64 scal = 1.0;
1181 
1182  if( o1 != o2 )
1183  {
1184 
1185  // scalar product
1186 
1187  real64 detJf2 = LvArray::math::abs( faceJacobianDeterminant( o2, X ));
1188 
1189  real64 el2[3][2] = { { edgeLength2( x1, x2, X ), edgeLength2( o1, o2, X ) },
1190 
1191  { edgeLength2( x1, o1, X ), edgeLength2( x2, o2, X ) },
1192 
1193  { edgeLength2( x1, o2, X ), edgeLength2( x2, o1, X ) } };
1194 
1195  real64 h2 = (el2[1][0]+el2[1][1])-(el2[2][0]+el2[2][1]);
1196 
1197  h2 = (4.0 * el2[0][0]*el2[0][1] - h2 * h2)/16.0;
1198 
1199  scal = (4.0*h2-detJf1 * detJf1 - detJf2 * detJf2 ) / (2.0 * detJf1 * detJf2);
1200 
1201  }
1202 
1203  return scal * detJf1 / detJ;
1204 
1205  }
1206 
1213  template< typename FUNC >
1214  GEOS_HOST_DEVICE
1215  GEOS_FORCE_INLINE
1216  static
1217  constexpr
1218  void
1219  computeStiffnessTerm( real64 const (&X)[4][3],
1220  FUNC && func )
1221  {
1222  real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1223  real64 dLambdadX[4][3] = {};
1224  for( int j = 0; j < 3; j++ )
1225  {
1226  dLambdadX[1][j] =
1227  ( ( X[ 2 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 3 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) - ( X[ 3 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 2 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) ) / detJ;
1228  dLambdadX[2][j] =
1229  ( ( X[ 3 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 1 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) - ( X[ 1 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 3 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) ) / detJ;
1230  dLambdadX[3][j] =
1231  ( ( X[ 1 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 2 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) - ( X[ 2 ][ (j+1)%3 ] - X[ 0 ][ (j+1)%3 ]) * ( X[ 1 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) ) / detJ;
1232  dLambdadX[0][j] = -dLambdadX[1][j] - dLambdadX[2][j] - dLambdadX[3][j];
1233  }
1234  basisLoop( [&func, &dLambdadX, &detJ] ( auto const cc1, auto const ci1, auto const cj1, auto const ck1, auto const cl1 )
1235  {
1236 
1237  constexpr int c1 = decltype(cc1)::value;
1238  constexpr int i1 = decltype(ci1)::value;
1239  constexpr int j1 = decltype(cj1)::value;
1240  constexpr int k1 = decltype(ck1)::value;
1241  constexpr int l1 = decltype(cl1)::value;
1242  //Not used in some combinations, but needed for constexpr
1243  GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
1244  basisLoop( [&func, &dLambdadX, &detJ] ( auto const cc2, auto const ci2, auto const cj2, auto const ck2, auto const cl2 )
1245  {
1246  constexpr int c2 = decltype(cc2)::value;
1247  constexpr int i2 = decltype(ci2)::value;
1248  constexpr int j2 = decltype(cj2)::value;
1249  constexpr int k2 = decltype(ck2)::value;
1250  constexpr int l2 = decltype(cl2)::value;
1251  //Not used in some combinations, but needed for constexpr
1252  GEOS_UNUSED_VAR( c2, i2, j2, k2, l2 );
1253  barycentricCoordinateLoop( [&func, &dLambdadX, &detJ] ( auto const cd1 )
1254  {
1255  constexpr int d1 = decltype(cd1)::value;
1256  barycentricCoordinateLoop( [&func, &dLambdadX, &detJ] ( auto const cd2 )
1257  {
1258  constexpr int d2 = decltype(cd2)::value;
1259  constexpr int ii1 = i1 + ( d1 == 0 ) * ( -1 );
1260  constexpr int ij1 = j1 + ( d1 == 1 ) * ( -1 );
1261  constexpr int ik1 = k1 + ( d1 == 2 ) * ( -1 );
1262  constexpr int il1 = l1 + ( d1 == 3 ) * ( -1 );
1263  constexpr int ii2 = i2 + ( d2 == 0 ) * ( -1 );
1264  constexpr int ij2 = j2 + ( d2 == 1 ) * ( -1 );
1265  constexpr int ik2 = k2 + ( d2 == 2 ) * ( -1 );
1266  constexpr int il2 = l2 + ( d2 == 3 ) * ( -1 );
1267  constexpr real64 factor1 = correctionFactorDerivative( i1, j1, k1, l1, ii1, ij1, ik1, il1, 3 );
1268  constexpr real64 factor2 = correctionFactorDerivative( i2, j2, k2, l2, ii2, ij2, ik2, il2, 3 );
1269  if constexpr (ii1 >= 0 && ij1 >= 0 && ik1 >= 0 && il1 >= 0 &&
1270  ii2 >= 0 && ij2 >= 0 && ik2 >= 0 && il2 >= 0)
1271  {
1272  constexpr real64 val = computeSuperpositionIntegral( ii1, ij1, ik1, il1, ii2, ij2, ik2, il2 ) * factor1 * factor2;
1273  func( c1, c2, val * detJ * ( dLambdadX[d1][0]*dLambdadX[d2][0] + dLambdadX[d1][1]*dLambdadX[d2][1] + dLambdadX[d1][2]*dLambdadX[d2][2] ) );
1274  }
1275  } );
1276  } );
1277  } );
1278  } );
1279  }
1280 
1281 
1289  GEOS_HOST_DEVICE
1290  GEOS_FORCE_INLINE
1291  static real64 edgeLength2( localIndex i1,
1292  localIndex i2,
1293  real64 const (&X)[4][3] )
1294  {
1295  real64 ab[3] = { X[ i2 ][ 0 ] - X[ i1 ][ 0 ],
1296  X[ i2 ][ 1 ] - X[ i1 ][ 1 ],
1297  X[ i2 ][ 2 ] - X[ i1 ][ 2 ] };
1298  return ab[0] * ab[0] + ab[1] * ab[1]+ ab[2]*ab[2];
1299  }
1300 
1315  template< typename FUNCP, typename FUNCF >
1316  GEOS_HOST_DEVICE
1317  GEOS_FORCE_INLINE
1318  static
1319  constexpr
1320  void
1321  computeSurfaceTerms( real64 const (&X)[4][3],
1322  FUNCP && funcP,
1323  FUNCF && funcF )
1324  {
1325  real64 detJf[4] = { faceJacobianDeterminant( 0, X ), faceJacobianDeterminant( 1, X ),
1326  faceJacobianDeterminant( 2, X ), faceJacobianDeterminant( 3, X ) };
1327  conditionalBasisLoop< 0, 1 >( [&] ( auto const cf1, auto const cd, auto const cc1, auto const ci1, auto const cj1, auto const ck1 )
1328  {
1329  conditionalBasisLoop< 0 >( [&] ( auto const cf2, auto const, auto const cc2, auto const ci2, auto const cj2, auto const ck2 )
1330  {
1331  constexpr int c1 = decltype(cc1)::value;
1332  constexpr int i1 = decltype(ci1)::value;
1333  constexpr int j1 = decltype(cj1)::value;
1334  constexpr int k1 = decltype(ck1)::value;
1335  // Not used in some combinations, but needed for constexpr
1336  GEOS_UNUSED_VAR( c1, i1, j1, k1 );
1337  constexpr int f2 = decltype(cf2)::value;
1338  constexpr int c2 = decltype(cc2)::value;
1339  constexpr int i2 = decltype(ci2)::value;
1340  constexpr int j2 = decltype(cj2)::value;
1341  constexpr int k2 = decltype(ck2)::value;
1342  // Not used in some combinations, but needed for constexpr
1343  GEOS_UNUSED_VAR( c2, i2, j2, k2, f2 );
1344 
1345  // The second function is nonzero on the face indexed by f2, so we integrate on this face.
1346  if constexpr ( std::decay_t< decltype(cf1) >::value == std::decay_t< decltype(cf2) >::value)
1347  {
1348  // compute penalty term iff the other function is also nonzero on the same face (i.e., d1==0)
1349  if constexpr ( std::decay_t< decltype(cd) >::value == 0 )
1350  {
1351  constexpr real64 val = computeFaceSuperpositionIntegral( i1, j1, k1, i2, j2, k2 );
1352  funcP( c1, c2, f2, i1, j1, k1, i2, j2, k2, val * detJf[ f2 ] );
1353  }
1354  // Compute flux term. This is nonzero in two cases.
1355  // first case: function has exponent 1 wrt to the same face. In this case, one can derive it once wrt to the
1356  // corresponding lambda and it will obtain a nonzero function on the face.
1357  if constexpr ( std::decay_t< decltype(cd) >::value == 1 )
1358  {
1359  constexpr real64 derFactor = ( i1 + j1 + k1 + 4 );
1360  constexpr real64 val = computeFaceSuperpositionIntegral( i1, j1, k1, i2, j2, k2 ) * derFactor;
1361  funcF( c1, c2, f2, -1, i1, j1, k1, i2, j2, k2, val * detJf[ f2 ] );
1362  }
1363  // second case: function has exponent zero wrt f2.
1364  // In this case, one can derive it wrt to any other face.
1365  else if constexpr ( std::decay_t< decltype(cd) >::value == 0 )
1366  {
1367  faceBarycentricCoordinateLoop( [ &funcF, &detJf ] ( auto const cl )
1368  {
1369  constexpr int l = decltype(cl)::value;
1370  constexpr int ii1 = i1 + ( l == 0 ) * ( -1 );
1371  constexpr int ij1 = j1 + ( l == 1 ) * ( -1 );
1372  constexpr int ik1 = k1 + ( l == 2 ) * ( -1 );
1373  if constexpr (ii1 >= 0 && ij1 >= 0 && ik1 >= 0)
1374  {
1375  constexpr real64 derFactor = ( ii1 + ij1 + ik1 + 4 );
1376  constexpr real64 val = computeFaceSuperpositionIntegral( ii1, ij1, ik1, i2, j2, k2 ) * derFactor;
1377  constexpr int f = l >= f2 ? l + 1 : l;
1378  funcF( c1, c2, f2, f, i1, j1, k1, i2, j2, k2, val * detJf[f2] );
1379  }
1380  } );
1381  }
1382  }
1383  } );
1384  } );
1385  }
1386 
1387 };
1388 
1390 template< int ORDER >
1391 class BB_Tetrahedron final : public BB_Tetrahedron_impl< ORDER >, public FiniteElementBase
1392 {
1393 public:
1394 
1396  using ImplType = BB_Tetrahedron_impl< ORDER >;
1397 
1398  BB_Tetrahedron():
1399  FiniteElementBase( ImplType::numNodes,
1400  ImplType::maxSupportPoints,
1401  ImplType::numQuadraturePoints )
1402  { }
1403 
1404 #ifdef __CUDACC__
1405  #pragma diag_push
1406  #pragma nv_diag_suppress 20012
1407 #endif
1408  GEOS_HOST_DEVICE
1409  virtual ~BB_Tetrahedron() override final = default;
1410 #ifdef __CUDACC__
1411  #pragma diag_pop
1412 #endif
1413 
1418  ImplType * getImpl()
1419  {
1420  return static_cast< ImplType * >(this);
1421  }
1422 
1427  ImplType const * getImpl() const
1428  {
1429  return static_cast< ImplType const * >(this);
1430  }
1431 
1432 };
1433 
1435 using BB1_Tetrahedron = BB_Tetrahedron< 1 >;
1437 using BB2_Tetrahedron = BB_Tetrahedron< 2 >;
1439 using BB3_Tetrahedron = BB_Tetrahedron< 3 >;
1440 
1441 //using BB4_Tetrahedron = BB_Tetrahedron< 4 >;
1442 //using BB5_Tetrahedron = BB_Tetrahedron< 5 >;
1443 
1444 
1446 using BB1_Tetrahedron_impl = BB_Tetrahedron_impl< 1 >;
1448 using BB2_Tetrahedron_impl = BB_Tetrahedron_impl< 2 >;
1450 using BB3_Tetrahedron_impl = BB_Tetrahedron_impl< 3 >;
1451 
1452 // using BB2_Tetrahedron_impl = BB_Tetrahedron_impl< 4 >;
1453 // using BB3_Tetrahedron_impl = BB_Tetrahedron_impl< 5 >;
1454 
1455 }
1456 }
1457 
1458 #endif // GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_BBTETRAHEDRON_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
#define GEOS_ERROR(...)
Raise a hard error and terminate the program.
Definition: Logger.hpp:204
geos::finiteElement::BB_Tetrahedron_impl
Definition: BB_Tetrahedron.hpp:54
geos::finiteElement::BB_Tetrahedron_impl::computeFluxDerivativeFactor
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE real64 computeFluxDerivativeFactor(real64 const (&X)[4][3], int x1, int x2, int o1, int o2)
Function to compute the factor for the flux derivative term.
Definition: BB_Tetrahedron.hpp:1171
geos::finiteElement::BB_Tetrahedron_impl::getMaxSupportPoints
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE localIndex getMaxSupportPoints()
Get the maximum number of support points.
Definition: BB_Tetrahedron.hpp:114
geos::finiteElement::BB_Tetrahedron_impl::calcNandGradN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcNandGradN(real64 const (&lambda)[4], real64 const (&N)[numNodes], real64(&gradN)[numNodes][4])
Calculate the values and derivatives of shape functions with respect to barycentric coordinates at a ...
Definition: BB_Tetrahedron.hpp:310
geos::finiteElement::BB_Tetrahedron_impl::integralTerm
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 integralTerm(const int a, const int b, const int c)
Computes a! / ( b! * c! ) with b + c >= a >= b, c.
Definition: BB_Tetrahedron.hpp:518
geos::finiteElement::BB_Tetrahedron_impl::getNumQuadraturePoints
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE localIndex getNumQuadraturePoints(StackVariables const &stack)
Get the number of quadrature points.
Definition: BB_Tetrahedron.hpp:97
geos::finiteElement::BB_Tetrahedron_impl::jacobianDeterminant
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 jacobianDeterminant(real64 const (&X)[4][3])
Returns the determinant of the Jacobian of the element.
Definition: BB_Tetrahedron.hpp:140
geos::finiteElement::BB_Tetrahedron_impl::indexToIJKL
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE int indexToIJKL()
Computes the local degree of freedom index given the shape function indices for each vertex.
Definition: BB_Tetrahedron.hpp:669
geos::finiteElement::BB_Tetrahedron_impl::computeStiffnessTerm
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE void computeStiffnessTerm(real64 const (&X)[4][3], FUNC &&func)
Computes the non-zero contributions inside the element of the stiffness matrix R, i....
Definition: BB_Tetrahedron.hpp:1219
geos::finiteElement::BB_Tetrahedron_impl::order
static constexpr int order
The order of the finite element.
Definition: BB_Tetrahedron.hpp:78
geos::finiteElement::BB_Tetrahedron_impl::computeReferenceMassMatrix
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeReferenceMassMatrix(arraySlice2d< real64 > const &m)
Computes the reference mass matrix, i.e., the superposition matrix of the shape functions in barycent...
Definition: BB_Tetrahedron.hpp:1051
geos::finiteElement::BB_Tetrahedron_impl::calcN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcN(localIndex const, real64(&N)[numNodes])
Calculate shape functions values for each support point at a quadrature point.
Definition: BB_Tetrahedron.hpp:189
geos::finiteElement::BB_Tetrahedron_impl::MeshData
typename Base::template MeshData< SubregionType > MeshData
Mesh data structure for the element.
Definition: BB_Tetrahedron.hpp:66
geos::finiteElement::BB_Tetrahedron_impl::calcN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcN(localIndex const q, StackVariables const &stack, real64(&N)[numNodes])
Calculate shape functions values for each support point at a quadrature point.
Definition: BB_Tetrahedron.hpp:209
geos::finiteElement::BB_Tetrahedron_impl::edgeLength2
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 edgeLength2(localIndex i1, localIndex i2, real64 const (&X)[4][3])
Definition: BB_Tetrahedron.hpp:1291
geos::finiteElement::BB_Tetrahedron_impl::computeSurfaceTerms
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE void computeSurfaceTerms(real64 const (&X)[4][3], FUNCP &&funcP, FUNCF &&funcF)
Computes the non-zero contributions inside the element of the surface terms, including the value of t...
Definition: BB_Tetrahedron.hpp:1321
geos::finiteElement::BB_Tetrahedron_impl::conditionalBasisLoop
static constexpr GEOS_HOST_DEVICE void conditionalBasisLoop(FUNC const &func)
Helper function for loop over tet basis functions that have one index in a given set of indices....
Definition: BB_Tetrahedron.hpp:865
geos::finiteElement::BB_Tetrahedron_impl::loop
static constexpr GEOS_HOST_DEVICE void loop(FUNC const &func)
Helper function for static for loop.
Definition: BB_Tetrahedron.hpp:725
geos::finiteElement::BB_Tetrahedron_impl::computeMassTerm
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE void computeMassTerm(real64 const (&X)[4][3], FUNC &&func)
Computes the non-zero contributions inside the element of the mass matrix M, i.e.,...
Definition: BB_Tetrahedron.hpp:1132
geos::finiteElement::BB_Tetrahedron_impl::calcGradN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcGradN(localIndex const q, real64 const (&X)[numNodes][3], StackVariables const &stack, real64(&gradN)[numNodes][3])
Calculate the shape functions derivatives wrt the physical coordinates.
Definition: BB_Tetrahedron.hpp:500
geos::finiteElement::BB_Tetrahedron_impl::basisLoop
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE void basisLoop(FUNC &&func)
Helper function for loop over tet basis functions.
Definition: BB_Tetrahedron.hpp:756
geos::finiteElement::BB_Tetrahedron_impl::getNumQuadraturePoints
static GEOS_HOST_DEVICE localIndex getNumQuadraturePoints()
Get the number of quadrature points.
Definition: BB_Tetrahedron.hpp:85
geos::finiteElement::BB_Tetrahedron_impl::call_matching_cases
static constexpr GEOS_HOST_DEVICE void call_matching_cases(F &&func, std::integer_sequence< int, Is... >)
Helper function for loop over tet basis functions that have one index in a given set of indices....
Definition: BB_Tetrahedron.hpp:803
geos::finiteElement::BB_Tetrahedron_impl::barycentricCoordinateLoop
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE void barycentricCoordinateLoop(FUNC &&func)
Helper function for loop over barycentric coordinates.
Definition: BB_Tetrahedron.hpp:743
geos::finiteElement::BB_Tetrahedron_impl::computeFaceSuperpositionIntegral
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 computeFaceSuperpositionIntegral(const int i1, const int j1, const int k1, const int i2, const int j2, const int k2)
Computes the superposition integral over a face between Bernstein-Bézier functions whose indices are ...
Definition: BB_Tetrahedron.hpp:602
geos::finiteElement::BB_Tetrahedron_impl::getNumSupportPoints
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE localIndex getNumSupportPoints()
Get the number of support points.
Definition: BB_Tetrahedron.hpp:106
geos::finiteElement::BB_Tetrahedron_impl::calcN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcN(real64 const (&X)[4][3], real64 const (&coords)[3], real64(&N)[numNodes])
Calculate shape functions values at a single point, given the coordinates of the tetrahedron vertices...
Definition: BB_Tetrahedron.hpp:271
geos::finiteElement::BB_Tetrahedron_impl::calcNandGradN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcNandGradN(real64 const (&X)[4][3], real64 const (&coords)[3], real64(&N)[numNodes], real64(&gradN)[numNodes][4])
Calculate the shape functions values and derivatives at a single point, given the coorginates of the ...
Definition: BB_Tetrahedron.hpp:369
geos::finiteElement::BB_Tetrahedron_impl::faceJacobianDeterminant
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 faceJacobianDeterminant(localIndex face, real64 const (&X)[4][3])
Calculate the determinant of the jacobian on the face opposite to the given vertex.
Definition: BB_Tetrahedron.hpp:161
geos::finiteElement::BB_Tetrahedron_impl::calcJacobian
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcJacobian(localIndex const q, real64 const (&X)[numNodes][3], StackVariables const &stack, real64(&J)[3][3])
Calculate the Jacobian matrix at a quadrature point.
Definition: BB_Tetrahedron.hpp:457
geos::finiteElement::BB_Tetrahedron_impl::dofIndex
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE int dofIndex()
Computes the local degree of freedom index given the shape function indices (i, j,...
Definition: BB_Tetrahedron.hpp:627
geos::finiteElement::BB_Tetrahedron_impl::getNumSupportPoints
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE localIndex getNumSupportPoints(StackVariables const &stack)
Get the number of support points.
Definition: BB_Tetrahedron.hpp:126
geos::finiteElement::BB_Tetrahedron_impl::calcN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcN(real64 const (&lambda)[4], real64(&N)[numNodes])
Calculate shape functions values at a single point using De Casteljau's algorithm.
Definition: BB_Tetrahedron.hpp:223
geos::finiteElement::BB_Tetrahedron_impl::loop_impl
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE void loop_impl(FUNC &&func, std::integer_sequence< int, Is... >)
Helper function for static for loop.
Definition: BB_Tetrahedron.hpp:711
geos::finiteElement::BB_Tetrahedron_impl::computeReferenceDampingMatrix
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void computeReferenceDampingMatrix(real64(&d)[numNodes][numNodes], bool const face1Damped, bool const face2Damped, bool const face3Damped, bool const face4Damped)
Computes the reference damping matrix, i.e., the superposition matrix of the shape functions in baryc...
Definition: BB_Tetrahedron.hpp:1091
geos::finiteElement::BB_Tetrahedron_impl::faceBarycentricCoordinateLoop
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE void faceBarycentricCoordinateLoop(FUNC &&func)
Helper function for loop over barycentric coordinates of a face.
Definition: BB_Tetrahedron.hpp:1035
geos::finiteElement::BB_Tetrahedron_impl::calcJacobian
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcJacobian(localIndex const q, real64 const (&X)[numNodes][3], real64(&J)[3][3])
Calculate the Jacobian matrix at a quadrature point.
Definition: BB_Tetrahedron.hpp:423
geos::finiteElement::BB_Tetrahedron_impl::dofIndex
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE int dofIndex(const int i, const int j, const int k)
Computes the local degree of freedom index given the shape function indices (i, j,...
Definition: BB_Tetrahedron.hpp:648
geos::finiteElement::BB_Tetrahedron_impl::calcGradN
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 calcGradN(localIndex const q, real64 const (&X)[numNodes][3], real64(&gradN)[numNodes][3])
Calculate the shape functions derivatives wrt the physical coordinates.
Definition: BB_Tetrahedron.hpp:477
geos::finiteElement::BB_Tetrahedron_impl::correctionFactorDerivative
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE real64 correctionFactorDerivative(int const i1, int const j1, int const k1, int const l1, int const i2, int const j2, int const k2, int const l2, int const dim)
Computes the correction factor for the superposition integral in case a function has been derived onc...
Definition: BB_Tetrahedron.hpp:557
geos::finiteElement::BB_Tetrahedron_impl::numNodesPerFace
constexpr static localIndex numNodesPerFace
The number of shape functions per face.
Definition: BB_Tetrahedron.hpp:81
geos::finiteElement::BB_Tetrahedron_impl::computeSuperpositionIntegral
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE real64 computeSuperpositionIntegral(const int i1, const int j1, const int k1, const int l1, const int i2, const int j2, const int k2, const int l2)
Computes the superposition integral between Bernstein-Bézier functions indexed by (i1,...
Definition: BB_Tetrahedron.hpp:578
geos::finiteElement::BB_Tetrahedron
Definition: BB_Tetrahedron.hpp:1392
geos::finiteElement::BB_Tetrahedron::getImpl
ImplType * getImpl()
Get the device-compatible implementation type.
Definition: BB_Tetrahedron.hpp:1418
geos::finiteElement::BB_Tetrahedron::getImpl
ImplType const * getImpl() const
Get the device-compatible implementation type.
Definition: BB_Tetrahedron.hpp:1427
geos::finiteElement::FiniteElementBase_impl
Device-compatible (non virtual) Base class for all finite element formulations.
Definition: FiniteElementBase.hpp:42
geos::finiteElement::FiniteElementBase_impl< BB_Tetrahedron_NumNodes< ORDER >, 4, BB_Tetrahedron_NumNodes< ORDER > >::numQuadraturePoints
constexpr static localIndex numQuadraturePoints
The number of quadrature points per element.
Definition: FiniteElementBase.hpp:58
geos::finiteElement::FiniteElementBase_impl< BB_Tetrahedron_NumNodes< ORDER >, 4, BB_Tetrahedron_NumNodes< ORDER > >::maxSupportPoints
constexpr static localIndex maxSupportPoints
The maximum number of support points per element.
Definition: FiniteElementBase.hpp:52
geos::finiteElement::FiniteElementBase_impl< BB_Tetrahedron_NumNodes< ORDER >, 4, BB_Tetrahedron_NumNodes< ORDER > >::numNodes
constexpr static localIndex numNodes
The number of nodes per element.
Definition: FiniteElementBase.hpp:46
geos::finiteElement::FiniteElementBase
Base class for FEM element implementations.
Definition: FiniteElementBase.hpp:417
geos::arraySlice2d
ArraySlice< T, 2, USD > arraySlice2d
Alias for 2D array slice.
Definition: DataTypes.hpp:199
geos::real64
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
geos::localIndex
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
StackVariables
Kernel variables (dof numbers, jacobian and residual) located on the stack.
Generated by doxygen 1.9.1