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 
44 template< int ORDER >
45 class BB_Tetrahedron final : public FiniteElementBase
46 {
47 public:
48 
50  static constexpr int order = ORDER;
51 
53  constexpr static localIndex numNodes = ( ORDER + 1 ) * ( ORDER + 2 ) * ( ORDER + 3 ) / 6;
54 
56  constexpr static localIndex numNodesPerFace = ( ORDER + 1 ) * ( ORDER + 2 ) / 2;
57 
59  constexpr static localIndex maxSupportPoints = numNodes;
60 
62  constexpr static localIndex numQuadraturePoints = numNodes;
63 
65  USING_FINITEELEMENTBASE
68  virtual ~BB_Tetrahedron() = default;
69 
70  GEOS_HOST_DEVICE
71  virtual localIndex getNumQuadraturePoints() const override
72  {
73  return numQuadraturePoints;
74  }
75 
81  GEOS_HOST_DEVICE
82  GEOS_FORCE_INLINE
83  static localIndex getNumQuadraturePoints( StackVariables const & stack )
84  {
85  GEOS_UNUSED_VAR( stack );
86  return numQuadraturePoints;
87  }
88 
89  GEOS_HOST_DEVICE
90  GEOS_FORCE_INLINE
91  virtual localIndex getNumSupportPoints() const override
92  {
93  return numNodes;
94  }
95 
96  GEOS_HOST_DEVICE
97  GEOS_FORCE_INLINE
98  virtual localIndex getMaxSupportPoints() const override
99  {
100  return maxSupportPoints;
101  }
102 
108  GEOS_HOST_DEVICE
109  GEOS_FORCE_INLINE
110  static localIndex getNumSupportPoints( StackVariables const & stack )
111  {
112  GEOS_UNUSED_VAR( stack );
113  return numNodes;
114  }
115 
116 
122  GEOS_HOST_DEVICE
123  GEOS_FORCE_INLINE
124  static real64 jacobianDeterminant( real64 const (&X)[4][3] )
125 
126  {
127  real64 m[3][3] = {};
128  for( int i = 0; i < 3; i++ )
129  {
130  for( int j = 0; j < 3; j++ )
131  {
132  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
133  }
134  }
135  return LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) );
136  }
137 
144  GEOS_HOST_DEVICE
145  GEOS_FORCE_INLINE
146  static real64 faceJacobianDeterminant( localIndex face,
147  real64 const (&X)[4][3] )
148  {
149  int i1 = ( face + 1 ) % 4;
150  int i2 = ( face + 2 ) % 4;
151  int i3 = ( face + 3 ) % 4;
152 
153  real64 ab[3] = { X[ i2 ][ 0 ] - X[ i1 ][ 0 ],
154  X[ i2 ][ 1 ] - X[ i1 ][ 1 ],
155  X[ i2 ][ 2 ] - X[ i1 ][ 2 ] };
156  real64 ac[3] = { X[ i3 ][ 0 ] - X[ i1 ][ 0 ],
157  X[ i3 ][ 1 ] - X[ i1 ][ 1 ],
158  X[ i3 ][ 2 ] - X[ i1 ][ 2 ] };
159  real64 term1 = ab[1] * ac[2] - ab[2] * ac[1];
160  real64 term2 = ab[2] * ac[0] - ab[0] * ac[2];
161  real64 term3 = ab[0] * ac[1] - ab[1] * ac[0];
162  return LvArray::math::sqrt( ( term1 * term1 + term2 * term2 + term3 * term3 ) );
163 
164  }
165 
172  GEOS_HOST_DEVICE
173  GEOS_FORCE_INLINE
174  static void calcN( localIndex const,
175  real64 (& N)[numNodes] )
176  {
177  for( int a=0; a < numNodes; ++a )
178  {
179  N[ a ] = 0;
180  }
181  GEOS_ERROR( "Bernstein-Bézier basis is modal, not nodal. No quadrature points are defined." );
182  }
183 
192  GEOS_HOST_DEVICE
193  GEOS_FORCE_INLINE
194  static void calcN( localIndex const q,
195  StackVariables const & stack,
196  real64 ( & N )[numNodes] )
197  {
198  GEOS_UNUSED_VAR( stack );
199  return calcN( q, N );
200  }
206  GEOS_HOST_DEVICE
207  GEOS_FORCE_INLINE
208  static void calcN( real64 const ( &lambda)[4],
209  real64 (& N)[numNodes] )
210  {
211  N[ 0 ] = 6.0;
212  int prev;
213  int c;
214  int limits[ 4 ] = { 1, 1, 1, 1 };
215  for( int np = 1; np <= ORDER; np++ )
216  {
217  prev = np * ( np + 1 ) * ( np + 2 ) / 6 - 1;
218  c = ( np + 1 ) * ( np + 2 ) * ( np + 3 ) / 6 - 1;
219  for( int i = 0; i < 4; i++ )
220  {
221  int denominator = i == 0 ? np : 1;
222  int offset = 0;
223  int c1 = np - 1;
224  int c2 = i + np - 2;
225  int repetitionCount = i == 0 ? 1 : limits[ i - 1 ];
226  for( int j = 0; j < limits[ i ]; j++ )
227  {
228  if( j - offset >= repetitionCount )
229  {
230  denominator++;
231  offset += repetitionCount;
232  repetitionCount = repetitionCount * c1 / c2;
233  c1--;
234  c2--;
235  }
236  N[ c-- ] = N[ prev - j ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
237  }
238  }
239  for( int i = 1; i < 4; i++ )
240  {
241  limits[ i ] += limits[ i - 1 ];
242  }
243  }
244  }
245 
254  GEOS_HOST_DEVICE
255  GEOS_FORCE_INLINE
256  static void calcN( real64 const (&X)[4][3],
257  real64 const (&coords)[3],
258  real64 (& N)[numNodes] )
259  {
260  real64 lambda[4] = {};
261  real64 m[3][3] = {};
262  for( int i = 0; i < 3; i++ )
263  {
264  for( int j = 0; j < 3; j++ )
265  {
266  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
267  }
268  }
269  real64 den = LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) );
270  for( int i = 0; i < 3; i++ )
271  {
272  for( int j = 0; j < 3; j++ )
273  {
274  m[ i ][ j ] = coords[ j ] - X[ 0 ][ j ];
275  }
276  lambda[ i + 1 ] = LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) ) / den;
277  for( int j = 0; j < 3; j++ )
278  {
279  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
280  }
281  }
282  lambda[ 0 ] = 1.0 - lambda[ 1 ] - lambda[ 2 ] - lambda[ 3 ];
283  return calcN( lambda, N );
284  }
285 
293  GEOS_HOST_DEVICE
294  GEOS_FORCE_INLINE
295  static void calcNandGradN( real64 const ( &lambda)[4],
296  real64 const ( &N)[numNodes],
297  real64 (& gradN)[numNodes][ 4 ] )
298  {
299  gradN[ 0 ][ 0 ] = 0.0;
300  gradN[ 0 ][ 1 ] = 0.0;
301  gradN[ 0 ][ 2 ] = 0.0;
302  gradN[ 0 ][ 3 ] = 0.0;
303  N[ 0 ] = 6.0;
304  int prev;
305  int c;
306  int limits[ 4 ] = { 1, 1, 1, 1 };
307  for( int np = 1; np <= ORDER; np++ )
308  {
309  prev = np * ( np + 1 ) * ( np + 2 ) / 6 - 1;
310  c = ( np + 1 ) * ( np + 2 ) * ( np + 3 ) / 6 - 1;
311  for( int i = 0; i < 4; i++ )
312  {
313  int denominator = i == 0 ? np : 1;
314  int offset = 0;
315  int c1 = np - 1;
316  int c2 = i + np - 2;
317  int repetitionCount = i == 0 ? 1 : limits[ i - 1 ];
318  for( int j = 0; j < limits[ i ]; j++ )
319  {
320  if( j - offset >= repetitionCount )
321  {
322  denominator++;
323  offset += repetitionCount;
324  repetitionCount = repetitionCount * c1 / c2;
325  c1--;
326  c2--;
327  }
328  gradN[ c ][ 0 ] = gradN[ prev - j ][ 0 ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
329  gradN[ c ][ 1 ] = gradN[ prev - j ][ 1 ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
330  gradN[ c ][ 2 ] = gradN[ prev - j ][ 2 ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
331  gradN[ c ][ 3 ] = gradN[ prev - j ][ 3 ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
332  gradN[ c ][ 3 - i ] += N[ prev - j ] * ( np + 3 ) / denominator;
333  N[ c-- ] = N[ prev - j ] * lambda[ 3 - i ] * ( np + 3 ) / denominator;
334  }
335  }
336  for( int i = 1; i < 4; i++ )
337  {
338  limits[ i ] += limits[ i - 1 ];
339  }
340  }
341  }
342 
343 
352  GEOS_HOST_DEVICE
353  GEOS_FORCE_INLINE
354  static void calcNandGradN( real64 const (&X)[4][3],
355  real64 const (&coords)[3],
356  real64 (& N)[numNodes],
357  real64 (& gradN)[numNodes][ 4 ] )
358  {
359  real64 lambda[4] = {};
360  real64 m[3][3] = {};
361  real64 dNdLambda[numNodes][4];
362  for( int i = 0; i < 3; i++ )
363  {
364  for( int j = 0; j < 3; j++ )
365  {
366  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
367  }
368  }
369  real64 den = LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) );
370  for( int i = 0; i < 3; i++ )
371  {
372  for( int j = 0; j < 3; j++ )
373  {
374  m[ i ][ j ] = coords[ j ] - X[ 0 ][ j ];
375  }
376  lambda[ i + 1 ] = LvArray::math::abs( LvArray::tensorOps::determinant< 3 >( m ) ) / den;
377  for( int j = 0; j < 3; j++ )
378  {
379  m[ i ][ j ] = X[ i + 1 ][ j ] - X[ 0 ][ j ];
380  }
381  }
382  lambda[ 0 ] = 1.0 - lambda[ 1 ] - lambda[ 2 ] - lambda[ 3 ];
383  calcNandGradN( lambda, N, dNdLambda );
384  for( int i = 0; i < numNodes; i++ )
385  {
386  for( int j = 0; j < 3; j++ )
387  {
388  gradN[ i ][ j ] =
389  ( ( ( 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 ] ) )*
390  ( dNdLambda[ i ][ 1 ] - dNdLambda[ i ][ 0 ] ) +
391  ( ( 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 ]) *
392  ( X[ 3 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) )* ( dNdLambda[ i ][ 2 ] - dNdLambda[ i ][ 0 ] ) +
393  ( ( 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 ]) *
394  ( X[ 1 ][ (j+2)%3 ] - X[ 0 ][ (j+2)%3 ] ) )* ( dNdLambda[ i ][ 3 ] - dNdLambda[ i ][ 0 ] )) / den;
395  }
396  }
397  }
398 
399 
409  GEOS_HOST_DEVICE
410  GEOS_FORCE_INLINE
411  static real64 calcGradN( localIndex const q,
412  real64 const (&X)[numNodes][3],
413  real64 ( & gradN )[numNodes][3] )
414  {
415  gradN[ q ][ 0 ] = 0.0;
416  gradN[ q ][ 1 ] = 0.0;
417  gradN[ q ][ 2 ] = 0.0;
418  GEOS_UNUSED_VAR( q, X, gradN );
419  GEOS_ERROR( "Bernstein-Bézier basis is modal, not nodal. No quadrature points are defined." );
420  return 0;
421  }
432  GEOS_HOST_DEVICE
433  GEOS_FORCE_INLINE
434  static real64 calcGradN( localIndex const q,
435  real64 const (&X)[numNodes][3],
436  StackVariables const & stack,
437  real64 ( & gradN )[numNodes][3] )
438  {
439  GEOS_UNUSED_VAR( stack );
440  return calcGradN( q, X, gradN );
441  }
442 
450  GEOS_HOST_DEVICE
451  GEOS_FORCE_INLINE
452  constexpr static real64 integralTerm( const int a, const int b, const int c )
453  {
454  real64 res = 1.0;
455  int num = a;
456  int den = c;
457  for( int i = b; i > 0; i-- )
458  {
459  res *= ( (real64) num ) / i;
460  num--;
461  }
462  for( int i = num; i > 0; i-- )
463  {
464  res *= ( (real64) i ) / den;
465  den--;
466  }
467  for( int i = den; i > 0; i-- )
468  {
469  res /= i;
470  }
471 
472  return res;
473  }
474 
489  GEOS_HOST_DEVICE
490  GEOS_FORCE_INLINE
491  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 )
492  {
493  return (i1+j1+k1+l1+dim)* (i2==0 ? 1 : i2) * (j2==0 ? 1 : j2)* (k2==0 ? 1 : k2)* (l2==0 ? 1 : l2);
494  }
495 
496 
510  GEOS_HOST_DEVICE
511  GEOS_FORCE_INLINE
512  static constexpr real64 computeSuperpositionIntegral( const int i1, const int j1, const int k1, const int l1,
513  const int i2, const int j2, const int k2, const int l2 )
514  {
515  return (integralTerm( i1+i2, i1, i2 )*
516  integralTerm( j1+j2, j1, j2 )*
517  integralTerm( k1+k2, k1, k2 )*
518  integralTerm( l1+l2, l1, l2 ))/
519  integralTerm( i1+j1+k1+l1+i2+j2+k2+l2+3,
520  i1+j1+k1+l1+3, i2+j2+k2+l2+3 );
521  }
522 
534  GEOS_HOST_DEVICE
535  GEOS_FORCE_INLINE
536  constexpr static real64 computeFaceSuperpositionIntegral( const int i1, const int j1, const int k1,
537  const int i2, const int j2, const int k2 )
538  {
539  return ((i1+k1+j1+3)*(i2+j2+k2+3)*
540  integralTerm( i1+i2, i1, i2 )*
541  integralTerm( j1+j2, j1, j2 )*
542  integralTerm( k1+k2, k1, k2 ))/
543  integralTerm( i1+j1+k1+i2+j2+k2+2,
544  i1+j1+k1+2, i2+j2+k2+2 );
545  }
546 
555  template< int I, int J, int K >
556  GEOS_HOST_DEVICE
557  GEOS_FORCE_INLINE
558  static
559  constexpr
560  int
561  dofIndex()
562  {
563  return ( ORDER - I ) * ( ORDER - I + 1 ) * ( ORDER - I + 2) / 6 +
564  ( ORDER - I - J ) * ( ORDER - I - J + 1 ) / 2 +
565  ( ORDER - I - J - K );
566  }
567 
577  GEOS_HOST_DEVICE
578  GEOS_FORCE_INLINE
579  static
580  constexpr
581  int
582  dofIndex( const int i, const int j, const int k )
583  {
584  return ( ORDER - i ) * ( ORDER - i + 1 ) * ( ORDER - i + 2) / 6 +
585  ( ORDER - i - j ) * ( ORDER - i - j + 1 ) / 2 +
586  ( ORDER - i - j - k );
587  }
588 
589 
590 
597  template< int C, int VTX >
598  GEOS_HOST_DEVICE
599  GEOS_FORCE_INLINE
600  static
601  constexpr
602  int
603  indexToIJKL()
604  {
605  static_assert( VTX >= 0 && VTX < 4 );
606  // compute the indices of c in the current element using tetrahedral and triangular roots
607  constexpr int cc1 = C + 1;
608  constexpr real64 tetr = cbrt( 3.0 * cc1 + sqrt( 9.0 * cc1 * cc1 - 1.0 / 27.0 ) )
609  + cbrt( 3.0 * cc1 - sqrt( 9.0 * cc1 * cc1 - 1.0 / 27.0 ) ) - 2;
610  constexpr int i = ORDER - round( tetr * 10.0 ) / 10;
611  constexpr int cc2 = C - ( ORDER - i ) * ( ORDER - i + 1 ) * ( ORDER - i + 2 ) / 6 + 1;
612  constexpr real64 trir = ( sqrt( 8.0 * cc2 + 1.0 ) - 1.0 ) / 2.0 - 1;
613  constexpr int j = ORDER - i - round( trir * 10.0 ) / 10;
614  constexpr int k = ORDER - i - j - ( C - (ORDER - i ) * ( ORDER - i + 1 ) * ( ORDER - i + 2 ) / 6
615  - ( ORDER - i - j ) * ( ORDER - i - j + 1 ) / 2 );
616  if constexpr ( VTX == 0 )
617  {
618  return i;
619  }
620  else if constexpr ( VTX == 1)
621  {
622  return j;
623  }
624  else if constexpr ( VTX == 2)
625  {
626  return k;
627  }
628  else if constexpr ( VTX == 3)
629  {
630  return ORDER - i - j - k;
631  }
632  return -1;
633  }
634 
642  template< typename FUNC, int... Is >
643  GEOS_HOST_DEVICE
644  GEOS_FORCE_INLINE
645  static constexpr void loop_impl( FUNC && func, std::integer_sequence< int, Is... > )
646  {
647  ( func( std::integral_constant< int, Is >{} ), ... );
648  }
649 
657  template< int N, typename FUNC >
658  GEOS_HOST_DEVICE
659  static constexpr void loop( FUNC const & func )
660  {
661  if constexpr (N > 0)
662  {
663  loop< N - 1 >( func );
664  func( std::integral_constant< int, N - 1 >{} );
665  }
666  }
667 
668 
669 
675  template< typename FUNC >
676  GEOS_HOST_DEVICE
677  GEOS_FORCE_INLINE
678  static constexpr void barycentricCoordinateLoop( FUNC && func )
679  {
680  loop_impl( std::forward< FUNC >( func ), std::make_integer_sequence< int, 4 >{} );
681  }
682 
688  template< typename FUNC >
689  GEOS_HOST_DEVICE
690  GEOS_FORCE_INLINE
691  static constexpr void basisLoop( FUNC && func )
692  {
693  loop< ORDER + 1 >( [&func] ( auto const iic )
694  {
695  constexpr int i = decltype(iic)::value;
696  constexpr int i1 = ORDER - i;
697  loop< ORDER + 1 >( [&func] ( auto const jjc )
698  {
699  constexpr int j = decltype(jjc)::value;
700  constexpr int j1 = ORDER - j;
701  if constexpr ( j1 <= ORDER - i1 )
702  {
703  loop< ORDER + 1 >( [&func] ( auto const kkc )
704  {
705  constexpr int k = decltype(kkc)::value;
706  constexpr int k1 = ORDER - k;
707  if constexpr ( k1 <= (ORDER - i1 - j1) )
708  {
709  constexpr int l1 = ORDER - i1 - j1 - k1;
710  constexpr int c1 = dofIndex< i1, j1, k1 >();
711  func( std::integral_constant< int, c1 >{},
712  std::integral_constant< int, i1 >{},
713  std::integral_constant< int, j1 >{},
714  std::integral_constant< int, k1 >{},
715  std::integral_constant< int, l1 >{} );
716  }
717  } );
718  }
719  } );
720  } );
721  }
722 
723 
737  template< int c1, int i1, int j1, int k1, int l1, typename F, int... Is >
738  GEOS_HOST_DEVICE
739  static constexpr void call_matching_cases( F && func, std::integer_sequence< int, Is... > )
740  {
741 
742  auto check_i1 = [&]( auto I ) {
743  if constexpr (i1 == decltype(I)::value) {
744  func( std::integral_constant< int, 0 >{},
745  std::integral_constant< int, i1 >{},
746  std::integral_constant< int, c1 >{},
747  std::integral_constant< int, j1 >{},
748  std::integral_constant< int, k1 >{},
749  std::integral_constant< int, l1 >{} );
750  }
751  };
752  (check_i1( std::integral_constant< int, Is >{} ), ...);
753 
754  auto check_j1 = [&]( auto I ) {
755  if constexpr (j1 == decltype(I)::value) {
756  func( std::integral_constant< int, 1 >{},
757  std::integral_constant< int, j1 >{},
758  std::integral_constant< int, c1 >{},
759  std::integral_constant< int, i1 >{},
760  std::integral_constant< int, k1 >{},
761  std::integral_constant< int, l1 >{} );
762  }
763  };
764  (check_j1( std::integral_constant< int, Is >{} ), ...);
765 
766  auto check_k1 = [&]( auto I ) {
767  if constexpr (k1 == decltype(I)::value) {
768  func( std::integral_constant< int, 2 >{},
769  std::integral_constant< int, k1 >{},
770  std::integral_constant< int, c1 >{},
771  std::integral_constant< int, i1 >{},
772  std::integral_constant< int, j1 >{},
773  std::integral_constant< int, l1 >{} );
774  }
775  };
776  (check_k1( std::integral_constant< int, Is >{} ), ...);
777 
778  auto check_l1 = [&]( auto I ) {
779  if constexpr (l1 == decltype(I)::value) {
780  func( std::integral_constant< int, 3 >{},
781  std::integral_constant< int, l1 >{},
782  std::integral_constant< int, c1 >{},
783  std::integral_constant< int, i1 >{},
784  std::integral_constant< int, j1 >{},
785  std::integral_constant< int, k1 >{} );
786  }
787  };
788  (check_l1( std::integral_constant< int, Is >{} ), ...);
789  }
790 
791 
799  template< int... Is, typename FUNC >
800  GEOS_HOST_DEVICE
801  static constexpr void conditionalBasisLoop( FUNC const & func )
802  {
803  loop< ORDER + 1 >( [&]( auto const i ) {
804  using i_t = decltype(i);
805  constexpr int iVal = i_t::value;
806  constexpr int i1 = ORDER - iVal;
807 
808  loop< ORDER + 1 >( [&]( auto const j ) {
809  using j_t = decltype(j);
810  constexpr int jVal = j_t::value;
811  constexpr int j1 = ORDER - jVal;
812 
813  constexpr bool valid_j1_i1 = (j1 <= iVal);
814  if constexpr (valid_j1_i1) {
815 
816  loop< ORDER + 1 >( [&]( auto const k ) {
817  using k_t = decltype(k);
818  constexpr int kVal = k_t::value;
819  constexpr int k1 = ORDER - kVal;
820 
821  constexpr bool valid_k1 = (k1 <= (ORDER - i1 - j1));
822  if constexpr (valid_k1) {
823  constexpr int l1 = ORDER - i1 - j1 - k1;
824  constexpr int c1 = dofIndex< i1, j1, k1 >();
825 
826  call_matching_cases< c1, i1, j1, k1, l1 >( func, std::integer_sequence< int, Is... >{} );
827 
828  // Pour chaque Is...
829  // (void)(((i1 == Is) &&
830  // (void(func(std::integral_constant<int, 0>{},
831  // std::integral_constant<int, i1>{},
832  // std::integral_constant<int, c1>{},
833  // std::integral_constant<int, j1>{},
834  // std::integral_constant<int, k1>{},
835  // std::integral_constant<int, l1>{})), 1)) || ...);
836 
837  // (void)(((j1 == Is) &&
838  // (void(func(std::integral_constant<int, 1>{},
839  // std::integral_constant<int, j1>{},
840  // std::integral_constant<int, c1>{},
841  // std::integral_constant<int, i1>{},
842  // std::integral_constant<int, k1>{},
843  // std::integral_constant<int, l1>{})), 1)) || ...);
844 
845  // (void)(((k1==Is) &&
846  // (void(func(std::integral_constant<int, 2>{},
847  // std::integral_constant<int, k1>{},
848  // std::integral_constant<int, c1>{},
849  // std::integral_constant<int, i1>{},
850  // std::integral_constant<int, j1>{},
851  // std::integral_constant<int, l1>{})), 1)) || ...);
852 
853  // (void)(((l1 ==Is) &&
854  // (void(func(std::integral_constant<int, 3>{},
855  // std::integral_constant<int, l1>{},
856  // std::integral_constant<int, c1>{},
857  // std::integral_constant<int, i1>{},
858  // std::integral_constant<int, j1>{},
859  // std::integral_constant<int, k1>{})), 1)) || ...);
860 
861 
862 
863  // (void)(((i1 ==Is) &&
864  // (void(func(
865  // HDInt<0>{},
866  // HDInt<i1>{},
867  // HDInt<c1>{},
868  // HDInt<j1>{},
869  // HDInt<k1>{},
870  // HDInt<l1>{})), 1)) || ...);
871 
872  // (void)(((j1==Is) &&
873  // (void(func(HDInt<1>{},
874  // HDInt<j1>{},
875  // HDInt<c1>{},
876  // HDInt<i1>{},
877  // HDInt<k1>{},
878  // HDInt<l1>{})), 1)) || ...);
879 
880  // (void)(((k1==Is) &&
881  // (void(func(
882  // HDInt<2>{},
883  // HDInt<k1>{},
884  // HDInt<c1>{},
885  // HDInt<i1>{},
886  // HDInt<j1>{},
887  // HDInt<l1>{})), 1)) || ...);
888 
889  // (void)(((l1==Is) &&
890  // (void(func(
891  // HDInt<3>{},
892  // HDInt<l1>{},
893  // HDInt<c1>{},
894  // HDInt<i1>{},
895  // HDInt<j1>{},
896  // HDInt<k1>{})), 1)) || ...);
897  // i1 matching
898  // callIfMatch<i1, Is...>([&] {
899  // func(HDInt<0>{}, HDInt<i1>{}, HDInt<c1>{}, HDInt<j1>{}, HDInt<k1>{}, HDInt<l1>{});
900  // });
901 
902  // callIfMatch<j1, Is...>([&] {
903  // func(HDInt<1>{}, HDInt<j1>{}, HDInt<c1>{}, HDInt<i1>{}, HDInt<k1>{}, HDInt<l1>{});
904  // });
905 
906  // callIfMatch<k1, Is...>([&] {
907  // func(HDInt<2>{}, HDInt<k1>{}, HDInt<c1>{}, HDInt<i1>{}, HDInt<j1>{}, HDInt<l1>{});
908  // });
909 
910  // callIfMatch<l1, Is...>([&] {
911  // func(HDInt<3>{}, HDInt<l1>{}, HDInt<c1>{}, HDInt<i1>{}, HDInt<j1>{}, HDInt<k1>{});
912  // });
913  // callIfMatch<i1, Is...>([&] {
914  // constexpr auto h0 = HDInt<0>{};
915  // constexpr auto hi1 = HDInt<i1>{};
916  // constexpr auto hc1 = HDInt<c1>{};
917  // constexpr auto hj1 = HDInt<j1>{};
918  // constexpr auto hk1 = HDInt<k1>{};
919  // constexpr auto hl1 = HDInt<l1>{};
920 
921  // func(h0, hi1, hc1, hj1, hk1, hl1);
922  // });
923 
924  // callIfMatch<j1, Is...>([&] {
925  // constexpr auto h1 = HDInt<1>{};
926  // constexpr auto hj1 = HDInt<j1>{};
927  // constexpr auto hc1 = HDInt<c1>{};
928  // constexpr auto hi1 = HDInt<i1>{};
929  // constexpr auto hk1 = HDInt<k1>{};
930  // constexpr auto hl1 = HDInt<l1>{};
931 
932  // func(h1, hj1, hc1, hi1, hk1, hl1);
933  // });
934 
935  // callIfMatch<k1, Is...>([&] {
936  // constexpr auto h2 = HDInt<2>{};
937  // constexpr auto hk1 = HDInt<k1>{};
938  // constexpr auto hc1 = HDInt<c1>{};
939  // constexpr auto hi1 = HDInt<i1>{};
940  // constexpr auto hj1 = HDInt<j1>{};
941  // constexpr auto hl1 = HDInt<l1>{};
942 
943  // func(h2, hk1, hc1, hi1, hj1, hl1);
944  // });
945 
946  // callIfMatch<l1, Is...>([&] {
947  // constexpr auto h3 = HDInt<3>{};
948  // constexpr auto hl1 = HDInt<l1>{};
949  // constexpr auto hc1 = HDInt<c1>{};
950  // constexpr auto hi1 = HDInt<i1>{};
951  // constexpr auto hj1 = HDInt<j1>{};
952  // constexpr auto hk1 = HDInt<k1>{};
953 
954  // func(h3, hl1, hc1, hi1, hj1, hk1);
955  // });
956 
957  }
958  } );
959  }
960  } );
961  } );
962  }
968  template< typename FUNC >
969  GEOS_HOST_DEVICE
970  GEOS_FORCE_INLINE
971  static constexpr void faceBarycentricCoordinateLoop( FUNC && func )
972  {
973  loop< 3 >( std::forward< FUNC >( func ) );
974  }
975 
983  GEOS_HOST_DEVICE
984  GEOS_FORCE_INLINE
985  static
986  void
987  computeReferenceMassMatrix( arraySlice2d< real64 > const & m )
988  {
989  basisLoop( [ &m ] ( auto const cc1, auto const ii1, auto const jj1, auto const kk1, auto const ll1 )
990  {
991  constexpr int c1 = decltype(cc1)::value;
992  constexpr int i1 = decltype(ii1)::value;
993  constexpr int j1 = decltype(jj1)::value;
994  constexpr int k1 = decltype(kk1)::value;
995  constexpr int l1 = decltype(ll1)::value;
996  // Needed for compilers that do not support constexpr lambdas
997  GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
998  basisLoop( [ &m ] ( auto const cc2, auto const ii2, auto const jj2, auto const kk2, auto const ll2 )
999  {
1000  constexpr int c2 = decltype(cc2)::value;
1001  constexpr int i2 = decltype(ii2)::value;
1002  constexpr int j2 = decltype(jj2)::value;
1003  constexpr int k2 = decltype(kk2)::value;
1004  constexpr int l2 = decltype(ll2)::value;
1005  constexpr real64 val = computeSuperpositionIntegral( i1, j1, k1, l1, i2, j2, k2, l2 );
1006  m[ c1 ][ c2 ] = val;
1007  } );
1008  } );
1009  }
1010 
1022  template< typename DAMPING >
1023  GEOS_HOST_DEVICE
1024  GEOS_FORCE_INLINE
1025  static
1026  void
1027  computeReferenceDampingMatrix( real64 (& d)[numNodes][numNodes], bool const face1Damped, bool const face2Damped, bool const face3Damped, bool const face4Damped )
1028  {
1029 
1030  for( int c1 = 0; c1 < numNodes; c1++ )
1031  {
1032  for( int c2 = 0; c2 < numNodes; c2++ )
1033  {
1034  d[ c1 ][ c2 ] = 0;
1035  }
1036  }
1037  conditionalBasisLoop< 0 >( [&] ( auto const f1, auto const, auto const c1, auto const i1, auto const j1, auto const k1 )
1038  {
1039  conditionalBasisLoop< 0 >( [&] ( auto const f2, auto const, auto const c2, auto const i2, auto const j2, auto const k2 )
1040  {
1041  if constexpr ( f1 == f2 )
1042  {
1043  constexpr real64 val = computeFaceSuperpositionIntegral( i1, j1, k1, i2, j2, k2 );
1044  if( ( f1 == 0 && face1Damped ) ||
1045  ( f1 == 1 && face2Damped ) ||
1046  ( f1 == 2 && face3Damped ) ||
1047  ( f1 == 3 && face4Damped ) )
1048  {
1049  d[ c1 ][ c2 ] += val;
1050  }
1051  }
1052  } );
1053  } );
1054  }
1055 
1062  template< typename FUNC >
1063  GEOS_HOST_DEVICE
1064  GEOS_FORCE_INLINE
1065  static
1066  constexpr
1067  void
1068  computeMassTerm( real64 const (&X)[4][3],
1069  FUNC && func )
1070  {
1071  // real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1072  // basisLoop( [&func, &detJ] ( auto const cc1, auto const ii1, auto const jj1, auto const kk1, auto const ll1 )
1073  // {
1074  // constexpr int c1 = cc1;
1075  // constexpr int i1 = ii1;
1076  // constexpr int j1 = jj1;
1077  // constexpr int k1 = kk1;
1078  // constexpr int l1 = ll1;
1079  // //Needed for compilors that do not support constexpr lambdas
1080  // GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
1081  // basisLoop( [&func, &detJ] ( auto const c2, auto const ii2, auto const jj2, auto const kk2, auto const ll2 )
1082  // {
1083  // constexpr int i2 = decltype(ii2)::value;
1084  // constexpr int j2 = decltype(jj2)::value;
1085  // constexpr int k2 = decltype(kk2)::value;
1086  // constexpr int l2 = decltype(ll2)::value;
1087  // constexpr real64 val = computeSuperpositionIntegral( i1, j1, k1, l1, i2, j2, k2, l2 );
1088  // func( c1, c2, val * detJ );
1089  // } );
1090  // } );
1091  real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1092  basisLoop( [&func, &detJ] ( auto cc1, auto ii1, auto jj1, auto kk1, auto ll1 )
1093  {
1094  constexpr int c1 = decltype(cc1)::value;
1095  constexpr int i1 = decltype(ii1)::value;
1096  constexpr int j1 = decltype(jj1)::value;
1097  constexpr int k1 = decltype(kk1)::value;
1098  constexpr int l1 = decltype(ll1)::value;
1099  GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
1100 
1101  basisLoop( [&func, &detJ, i1, j1, k1, l1] ( auto c2, auto ii2, auto jj2, auto kk2, auto ll2 )
1102  {
1103  constexpr int c2v = decltype(c2)::value;
1104  constexpr int i2 = decltype(ii2)::value;
1105  constexpr int j2 = decltype(jj2)::value;
1106  constexpr int k2 = decltype(kk2)::value;
1107  constexpr int l2 = decltype(ll2)::value;
1108 
1109  constexpr real64 val = computeSuperpositionIntegral( i1, j1, k1, l1, i2, j2, k2, l2 );
1110  func( c1, c2v, val * detJ );
1111  } );
1112  } );
1113  }
1114 
1124  GEOS_HOST_DEVICE
1125  GEOS_FORCE_INLINE
1126  static constexpr real64
1127  computeFluxDerivativeFactor( real64 const (&X)[4][3], int x1, int x2, int o1, int o2 ) // Order x1, x2, o1, o2
1128  {
1129 
1130  real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1131 
1132  // height with respect to o2
1133 
1134  real64 detJf1 = LvArray::math::abs( faceJacobianDeterminant( o1, X ));
1135 
1136  real64 scal = 1.0;
1137 
1138  if( o1 != o2 )
1139  {
1140 
1141  // scalar product
1142 
1143  real64 detJf2 = LvArray::math::abs( faceJacobianDeterminant( o2, X ));
1144 
1145  real64 el2[3][2] = { { edgeLength2( x1, x2, X ), edgeLength2( o1, o2, X ) },
1146 
1147  { edgeLength2( x1, o1, X ), edgeLength2( x2, o2, X ) },
1148 
1149  { edgeLength2( x1, o2, X ), edgeLength2( x2, o1, X ) } };
1150 
1151  real64 h2 = (el2[1][0]+el2[1][1])-(el2[2][0]+el2[2][1]);
1152 
1153  h2 = (4.0 * el2[0][0]*el2[0][1] - h2 * h2)/16.0;
1154 
1155  scal = (4.0*h2-detJf1 * detJf1 - detJf2 * detJf2 ) / (2.0 * detJf1 * detJf2);
1156 
1157  }
1158 
1159  return scal * detJf1 / detJ;
1160 
1161  }
1162 
1169  template< typename FUNC >
1170  GEOS_HOST_DEVICE
1171  GEOS_FORCE_INLINE
1172  static
1173  constexpr
1174  void
1175  computeStiffnessTerm( real64 const (&X)[4][3],
1176  FUNC && func )
1177  {
1178  real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1179  real64 dLambdadX[4][3] = {};
1180  for( int j = 0; j < 3; j++ )
1181  {
1182  dLambdadX[1][j] =
1183  ( ( 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;
1184  dLambdadX[2][j] =
1185  ( ( 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;
1186  dLambdadX[3][j] =
1187  ( ( 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;
1188  dLambdadX[0][j] = -dLambdadX[1][j] - dLambdadX[2][j] - dLambdadX[3][j];
1189  }
1190  basisLoop( [&func, &dLambdadX, &detJ] ( auto const cc1, auto const ci1, auto const cj1, auto const ck1, auto const cl1 )
1191  {
1192 
1193  constexpr int c1 = decltype(cc1)::value;
1194  constexpr int i1 = decltype(ci1)::value;
1195  constexpr int j1 = decltype(cj1)::value;
1196  constexpr int k1 = decltype(ck1)::value;
1197  constexpr int l1 = decltype(cl1)::value;
1198  //Not used in some combinations, but needed for constexpr
1199  GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
1200  basisLoop( [&func, &dLambdadX, &detJ] ( auto const cc2, auto const ci2, auto const cj2, auto const ck2, auto const cl2 )
1201  {
1202  constexpr int c2 = decltype(cc2)::value;
1203  constexpr int i2 = decltype(ci2)::value;
1204  constexpr int j2 = decltype(cj2)::value;
1205  constexpr int k2 = decltype(ck2)::value;
1206  constexpr int l2 = decltype(cl2)::value;
1207  //Not used in some combinations, but needed for constexpr
1208  GEOS_UNUSED_VAR( c2, i2, j2, k2, l2 );
1209  barycentricCoordinateLoop( [&func, &dLambdadX, &detJ] ( auto const cd1 )
1210  {
1211  constexpr int d1 = decltype(cd1)::value;
1212  barycentricCoordinateLoop( [&func, &dLambdadX, &detJ] ( auto const cd2 )
1213  {
1214  constexpr int d2 = decltype(cd2)::value;
1215  constexpr int ii1 = i1 + ( d1 == 0 ) * ( -1 );
1216  constexpr int ij1 = j1 + ( d1 == 1 ) * ( -1 );
1217  constexpr int ik1 = k1 + ( d1 == 2 ) * ( -1 );
1218  constexpr int il1 = l1 + ( d1 == 3 ) * ( -1 );
1219  constexpr int ii2 = i2 + ( d2 == 0 ) * ( -1 );
1220  constexpr int ij2 = j2 + ( d2 == 1 ) * ( -1 );
1221  constexpr int ik2 = k2 + ( d2 == 2 ) * ( -1 );
1222  constexpr int il2 = l2 + ( d2 == 3 ) * ( -1 );
1223  constexpr real64 factor1 = correctionFactorDerivative( i1, j1, k1, l1, ii1, ij1, ik1, il1, 3 );
1224  constexpr real64 factor2 = correctionFactorDerivative( i2, j2, k2, l2, ii2, ij2, ik2, il2, 3 );
1225  if constexpr (ii1 >= 0 && ij1 >= 0 && ik1 >= 0 && il1 >= 0 &&
1226  ii2 >= 0 && ij2 >= 0 && ik2 >= 0 && il2 >= 0)
1227  {
1228  constexpr real64 val = computeSuperpositionIntegral( ii1, ij1, ik1, il1, ii2, ij2, ik2, il2 ) * factor1 * factor2;
1229  func( c1, c2, val * detJ * ( dLambdadX[d1][0]*dLambdadX[d2][0] + dLambdadX[d1][1]*dLambdadX[d2][1] + dLambdadX[d1][2]*dLambdadX[d2][2] ) );
1230  }
1231  } );
1232  } );
1233  } );
1234  } );
1235  }
1236 
1237 
1245  GEOS_HOST_DEVICE
1246  GEOS_FORCE_INLINE
1247  static real64 edgeLength2( localIndex i1,
1248  localIndex i2,
1249  real64 const (&X)[4][3] )
1250  {
1251  real64 ab[3] = { X[ i2 ][ 0 ] - X[ i1 ][ 0 ],
1252  X[ i2 ][ 1 ] - X[ i1 ][ 1 ],
1253  X[ i2 ][ 2 ] - X[ i1 ][ 2 ] };
1254  return ab[0] * ab[0] + ab[1] * ab[1]+ ab[2]*ab[2];
1255  }
1256 
1271  template< typename FUNCP, typename FUNCF >
1272  GEOS_HOST_DEVICE
1273  GEOS_FORCE_INLINE
1274  static
1275  constexpr
1276  void
1277  computeSurfaceTerms( real64 const (&X)[4][3],
1278  FUNCP && funcP,
1279  FUNCF && funcF )
1280  {
1281  real64 detJf[4] = { faceJacobianDeterminant( 0, X ), faceJacobianDeterminant( 1, X ),
1282  faceJacobianDeterminant( 2, X ), faceJacobianDeterminant( 3, X ) };
1283  conditionalBasisLoop< 0, 1 >( [&] ( auto const cf1, auto const cd, auto const cc1, auto const ci1, auto const cj1, auto const ck1 )
1284  {
1285  conditionalBasisLoop< 0 >( [&] ( auto const cf2, auto const, auto const cc2, auto const ci2, auto const cj2, auto const ck2 )
1286  {
1287  constexpr int c1 = decltype(cc1)::value;
1288  constexpr int i1 = decltype(ci1)::value;
1289  constexpr int j1 = decltype(cj1)::value;
1290  constexpr int k1 = decltype(ck1)::value;
1291  // Not used in some combinations, but needed for constexpr
1292  GEOS_UNUSED_VAR( c1, i1, j1, k1 );
1293  constexpr int f2 = decltype(cf2)::value;
1294  constexpr int c2 = decltype(cc2)::value;
1295  constexpr int i2 = decltype(ci2)::value;
1296  constexpr int j2 = decltype(cj2)::value;
1297  constexpr int k2 = decltype(ck2)::value;
1298  // Not used in some combinations, but needed for constexpr
1299  GEOS_UNUSED_VAR( c2, i2, j2, k2, f2 );
1300 
1301  // The second function is nonzero on the face indexed by f2, so we integrate on this face.
1302  if constexpr ( std::decay_t< decltype(cf1) >::value == std::decay_t< decltype(cf2) >::value)
1303  {
1304  // compute penalty term iff the other function is also nonzero on the same face (i.e., d1==0)
1305  if constexpr ( std::decay_t< decltype(cd) >::value == 0 )
1306  {
1307  constexpr real64 val = computeFaceSuperpositionIntegral( i1, j1, k1, i2, j2, k2 );
1308  funcP( c1, c2, f2, i1, j1, k1, i2, j2, k2, val * detJf[ f2 ] );
1309  }
1310  // Compute flux term. This is nonzero in two cases.
1311  // first case: function has exponent 1 wrt to the same face. In this case, one can derive it once wrt to the
1312  // corresponding lambda and it will obtain a nonzero function on the face.
1313  if constexpr ( std::decay_t< decltype(cd) >::value == 1 )
1314  {
1315  constexpr real64 derFactor = ( i1 + j1 + k1 + 4 );
1316  constexpr real64 val = computeFaceSuperpositionIntegral( i1, j1, k1, i2, j2, k2 ) * derFactor;
1317  funcF( c1, c2, f2, -1, i1, j1, k1, i2, j2, k2, val * detJf[ f2 ] );
1318  }
1319  // second case: function has exponent zero wrt f2.
1320  // In this case, one can derive it wrt to any other face.
1321  else if constexpr ( std::decay_t< decltype(cd) >::value == 0 )
1322  {
1323  faceBarycentricCoordinateLoop( [ &funcF, &detJf ] ( auto const cl )
1324  {
1325  constexpr int l = decltype(cl)::value;
1326  constexpr int ii1 = i1 + ( l == 0 ) * ( -1 );
1327  constexpr int ij1 = j1 + ( l == 1 ) * ( -1 );
1328  constexpr int ik1 = k1 + ( l == 2 ) * ( -1 );
1329  if constexpr (ii1 >= 0 && ij1 >= 0 && ik1 >= 0)
1330  {
1331  constexpr real64 derFactor = ( ii1 + ij1 + ik1 + 4 );
1332  constexpr real64 val = computeFaceSuperpositionIntegral( ii1, ij1, ik1, i2, j2, k2 ) * derFactor;
1333  constexpr int f = l >= f2 ? l + 1 : l;
1334  funcF( c1, c2, f2, f, i1, j1, k1, i2, j2, k2, val * detJf[f2] );
1335  }
1336  } );
1337  }
1338  }
1339  } );
1340  } );
1341  }
1342 
1343 };
1344 
1348 using BB1_Tetrahedron = BB_Tetrahedron< 1 >;
1352 using BB2_Tetrahedron = BB_Tetrahedron< 2 >;
1356 using BB3_Tetrahedron = BB_Tetrahedron< 3 >;
1357 //using BB4_Tetrahedron = BB_Tetrahedron< 4 >;
1358 //using BB5_Tetrahedron = BB_Tetrahedron< 5 >;
1359 
1360 }
1361 }
1362 
1363 #endif // GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_BBTETRAHEDRON_HPP_
USING_FINITEELEMENTBASE
#define USING_FINITEELEMENTBASE
Macro to simplify name resolution in derived classes.
Definition: FiniteElementBase.hpp:1063
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(msg)
Raise a hard error and terminate the program.
Definition: Logger.hpp:157
geos::finiteElement::BB_Tetrahedron
Definition: BB_Tetrahedron.hpp:46
geos::finiteElement::BB_Tetrahedron::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:512
geos::finiteElement::BB_Tetrahedron::getMaxSupportPoints
GEOS_HOST_DEVICE virtual GEOS_FORCE_INLINE localIndex getMaxSupportPoints() const override
Get the maximum number of support points for this element.
Definition: BB_Tetrahedron.hpp:98
geos::finiteElement::BB_Tetrahedron::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:1277
geos::finiteElement::BB_Tetrahedron::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:208
geos::finiteElement::BB_Tetrahedron::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:434
geos::finiteElement::BB_Tetrahedron::edgeLength2
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE real64 edgeLength2(localIndex i1, localIndex i2, real64 const (&X)[4][3])
Definition: BB_Tetrahedron.hpp:1247
geos::finiteElement::BB_Tetrahedron::getNumQuadraturePoints
virtual GEOS_HOST_DEVICE localIndex getNumQuadraturePoints() const override
Virtual getter for the number of quadrature points per element.
Definition: BB_Tetrahedron.hpp:71
geos::finiteElement::BB_Tetrahedron::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:1027
geos::finiteElement::BB_Tetrahedron::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:801
geos::finiteElement::BB_Tetrahedron::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:452
geos::finiteElement::BB_Tetrahedron::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:603
geos::finiteElement::BB_Tetrahedron::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:971
geos::finiteElement::BB_Tetrahedron::getNumSupportPoints
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE localIndex getNumSupportPoints(StackVariables const &stack)
Get the number of support points.
Definition: BB_Tetrahedron.hpp:110
geos::finiteElement::BB_Tetrahedron::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:1127
geos::finiteElement::BB_Tetrahedron::numNodesPerFace
constexpr static localIndex numNodesPerFace
The number of shape functions per face.
Definition: BB_Tetrahedron.hpp:56
geos::finiteElement::BB_Tetrahedron::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:256
geos::finiteElement::BB_Tetrahedron::numQuadraturePoints
constexpr static localIndex numQuadraturePoints
The number of quadrature points per element.
Definition: BB_Tetrahedron.hpp:62
geos::finiteElement::BB_Tetrahedron::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:411
geos::finiteElement::BB_Tetrahedron::maxSupportPoints
constexpr static localIndex maxSupportPoints
The maximum number of support points per element.
Definition: BB_Tetrahedron.hpp:59
geos::finiteElement::BB_Tetrahedron::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:354
geos::finiteElement::BB_Tetrahedron::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:194
geos::finiteElement::BB_Tetrahedron::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:174
geos::finiteElement::BB_Tetrahedron::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:739
geos::finiteElement::BB_Tetrahedron::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:124
geos::finiteElement::BB_Tetrahedron::order
static constexpr int order
The order of the finite element.
Definition: BB_Tetrahedron.hpp:50
geos::finiteElement::BB_Tetrahedron::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:491
geos::finiteElement::BB_Tetrahedron::barycentricCoordinateLoop
GEOS_HOST_DEVICE static constexpr GEOS_FORCE_INLINE void barycentricCoordinateLoop(FUNC &&func)
Helper function for loop over barycentric coordinates.
Definition: BB_Tetrahedron.hpp:678
geos::finiteElement::BB_Tetrahedron::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:645
geos::finiteElement::BB_Tetrahedron::~BB_Tetrahedron
virtual ~BB_Tetrahedron()=default
geos::finiteElement::BB_Tetrahedron::getNumQuadraturePoints
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE localIndex getNumQuadraturePoints(StackVariables const &stack)
Get the number of quadrature points.
Definition: BB_Tetrahedron.hpp:83
geos::finiteElement::BB_Tetrahedron::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:295
geos::finiteElement::BB_Tetrahedron::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:146
geos::finiteElement::BB_Tetrahedron::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:536
geos::finiteElement::BB_Tetrahedron::getNumSupportPoints
GEOS_HOST_DEVICE virtual GEOS_FORCE_INLINE localIndex getNumSupportPoints() const override
Virtual getter for the number of support points per element.
Definition: BB_Tetrahedron.hpp:91
geos::finiteElement::BB_Tetrahedron::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:987
geos::finiteElement::BB_Tetrahedron::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:582
geos::finiteElement::BB_Tetrahedron::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:1175
geos::finiteElement::BB_Tetrahedron::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:561
geos::finiteElement::BB_Tetrahedron::loop
static constexpr GEOS_HOST_DEVICE void loop(FUNC const &func)
Helper function for static for loop.
Definition: BB_Tetrahedron.hpp:659
geos::finiteElement::BB_Tetrahedron::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:691
geos::finiteElement::BB_Tetrahedron::numNodes
constexpr static localIndex numNodes
The number of shape functions per element.
Definition: BB_Tetrahedron.hpp:53
geos::finiteElement::BB_Tetrahedron::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:1068
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], real64 &value)
Compute the interpolated value of a variable.
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
Generated by doxygen 1.9.1