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 
674  template< typename FUNC >
675  GEOS_HOST_DEVICE
676  GEOS_FORCE_INLINE
677  static constexpr void barycentricCoordinateLoop( FUNC && func )
678  {
679  loop_impl( std::forward< FUNC >( func ), std::make_integer_sequence< int, 4 >{} );
680  }
681 
687  template< typename FUNC >
688  GEOS_HOST_DEVICE
689  GEOS_FORCE_INLINE
690  static constexpr void basisLoop( FUNC && func )
691  {
692  loop< ORDER + 1 >( [&func] ( auto const iic )
693  {
694  constexpr int i = decltype(iic)::value;
695  constexpr int i1 = ORDER - i;
696  loop< ORDER + 1 >( [&func] ( auto const jjc )
697  {
698  constexpr int j = decltype(jjc)::value;
699  constexpr int j1 = ORDER - j;
700  if constexpr ( j1 <= ORDER - i1 )
701  {
702  loop< ORDER + 1 >( [&func] ( auto const kkc )
703  {
704  constexpr int k = decltype(kkc)::value;
705  constexpr int k1 = ORDER - k;
706  if constexpr ( k1 <= (ORDER - i1 - j1) )
707  {
708  constexpr int l1 = ORDER - i1 - j1 - k1;
709  constexpr int c1 = dofIndex< i1, j1, k1 >();
710  func( std::integral_constant< int, c1 >{},
711  std::integral_constant< int, i1 >{},
712  std::integral_constant< int, j1 >{},
713  std::integral_constant< int, k1 >{},
714  std::integral_constant< int, l1 >{} );
715  }
716  } );
717  }
718  } );
719  } );
720  }
721 
735  template< int c1, int i1, int j1, int k1, int l1, typename F, int... Is >
736  GEOS_HOST_DEVICE
737  static constexpr void call_matching_cases( F && func, std::integer_sequence< int, Is... > )
738  {
739 
740  auto check_i1 = [&]( auto I ) {
741  if constexpr (i1 == decltype(I)::value) {
742  func( std::integral_constant< int, 0 >{},
743  std::integral_constant< int, i1 >{},
744  std::integral_constant< int, c1 >{},
745  std::integral_constant< int, j1 >{},
746  std::integral_constant< int, k1 >{},
747  std::integral_constant< int, l1 >{} );
748  }
749  };
750  (check_i1( std::integral_constant< int, Is >{} ), ...);
751 
752  auto check_j1 = [&]( auto I ) {
753  if constexpr (j1 == decltype(I)::value) {
754  func( std::integral_constant< int, 1 >{},
755  std::integral_constant< int, j1 >{},
756  std::integral_constant< int, c1 >{},
757  std::integral_constant< int, i1 >{},
758  std::integral_constant< int, k1 >{},
759  std::integral_constant< int, l1 >{} );
760  }
761  };
762  (check_j1( std::integral_constant< int, Is >{} ), ...);
763 
764  auto check_k1 = [&]( auto I ) {
765  if constexpr (k1 == decltype(I)::value) {
766  func( std::integral_constant< int, 2 >{},
767  std::integral_constant< int, k1 >{},
768  std::integral_constant< int, c1 >{},
769  std::integral_constant< int, i1 >{},
770  std::integral_constant< int, j1 >{},
771  std::integral_constant< int, l1 >{} );
772  }
773  };
774  (check_k1( std::integral_constant< int, Is >{} ), ...);
775 
776  auto check_l1 = [&]( auto I ) {
777  if constexpr (l1 == decltype(I)::value) {
778  func( std::integral_constant< int, 3 >{},
779  std::integral_constant< int, l1 >{},
780  std::integral_constant< int, c1 >{},
781  std::integral_constant< int, i1 >{},
782  std::integral_constant< int, j1 >{},
783  std::integral_constant< int, k1 >{} );
784  }
785  };
786  (check_l1( std::integral_constant< int, Is >{} ), ...);
787  }
788 
789 
797  template< int... Is, typename FUNC >
798  GEOS_HOST_DEVICE
799  static constexpr void conditionalBasisLoop( FUNC const & func )
800  {
801  loop< ORDER + 1 >( [&]( auto const i ) {
802  using i_t = decltype(i);
803  constexpr int iVal = i_t::value;
804  constexpr int i1 = ORDER - iVal;
805 
806  loop< ORDER + 1 >( [&]( auto const j ) {
807  using j_t = decltype(j);
808  constexpr int jVal = j_t::value;
809  constexpr int j1 = ORDER - jVal;
810 
811  constexpr bool valid_j1_i1 = (j1 <= iVal);
812  if constexpr (valid_j1_i1) {
813 
814  loop< ORDER + 1 >( [&]( auto const k ) {
815  using k_t = decltype(k);
816  constexpr int kVal = k_t::value;
817  constexpr int k1 = ORDER - kVal;
818 
819  constexpr bool valid_k1 = (k1 <= (ORDER - i1 - j1));
820  if constexpr (valid_k1) {
821  constexpr int l1 = ORDER - i1 - j1 - k1;
822  constexpr int c1 = dofIndex< i1, j1, k1 >();
823 
824  call_matching_cases< c1, i1, j1, k1, l1 >( func, std::integer_sequence< int, Is... >{} );
825 
826  // Pour chaque Is...
827  // (void)(((i1 == Is) &&
828  // (void(func(std::integral_constant<int, 0>{},
829  // std::integral_constant<int, i1>{},
830  // std::integral_constant<int, c1>{},
831  // std::integral_constant<int, j1>{},
832  // std::integral_constant<int, k1>{},
833  // std::integral_constant<int, l1>{})), 1)) || ...);
834 
835  // (void)(((j1 == Is) &&
836  // (void(func(std::integral_constant<int, 1>{},
837  // std::integral_constant<int, j1>{},
838  // std::integral_constant<int, c1>{},
839  // std::integral_constant<int, i1>{},
840  // std::integral_constant<int, k1>{},
841  // std::integral_constant<int, l1>{})), 1)) || ...);
842 
843  // (void)(((k1==Is) &&
844  // (void(func(std::integral_constant<int, 2>{},
845  // std::integral_constant<int, k1>{},
846  // std::integral_constant<int, c1>{},
847  // std::integral_constant<int, i1>{},
848  // std::integral_constant<int, j1>{},
849  // std::integral_constant<int, l1>{})), 1)) || ...);
850 
851  // (void)(((l1 ==Is) &&
852  // (void(func(std::integral_constant<int, 3>{},
853  // std::integral_constant<int, l1>{},
854  // std::integral_constant<int, c1>{},
855  // std::integral_constant<int, i1>{},
856  // std::integral_constant<int, j1>{},
857  // std::integral_constant<int, k1>{})), 1)) || ...);
858 
859 
860 
861  // (void)(((i1 ==Is) &&
862  // (void(func(
863  // HDInt<0>{},
864  // HDInt<i1>{},
865  // HDInt<c1>{},
866  // HDInt<j1>{},
867  // HDInt<k1>{},
868  // HDInt<l1>{})), 1)) || ...);
869 
870  // (void)(((j1==Is) &&
871  // (void(func(HDInt<1>{},
872  // HDInt<j1>{},
873  // HDInt<c1>{},
874  // HDInt<i1>{},
875  // HDInt<k1>{},
876  // HDInt<l1>{})), 1)) || ...);
877 
878  // (void)(((k1==Is) &&
879  // (void(func(
880  // HDInt<2>{},
881  // HDInt<k1>{},
882  // HDInt<c1>{},
883  // HDInt<i1>{},
884  // HDInt<j1>{},
885  // HDInt<l1>{})), 1)) || ...);
886 
887  // (void)(((l1==Is) &&
888  // (void(func(
889  // HDInt<3>{},
890  // HDInt<l1>{},
891  // HDInt<c1>{},
892  // HDInt<i1>{},
893  // HDInt<j1>{},
894  // HDInt<k1>{})), 1)) || ...);
895  // i1 matching
896  // callIfMatch<i1, Is...>([&] {
897  // func(HDInt<0>{}, HDInt<i1>{}, HDInt<c1>{}, HDInt<j1>{}, HDInt<k1>{}, HDInt<l1>{});
898  // });
899 
900  // callIfMatch<j1, Is...>([&] {
901  // func(HDInt<1>{}, HDInt<j1>{}, HDInt<c1>{}, HDInt<i1>{}, HDInt<k1>{}, HDInt<l1>{});
902  // });
903 
904  // callIfMatch<k1, Is...>([&] {
905  // func(HDInt<2>{}, HDInt<k1>{}, HDInt<c1>{}, HDInt<i1>{}, HDInt<j1>{}, HDInt<l1>{});
906  // });
907 
908  // callIfMatch<l1, Is...>([&] {
909  // func(HDInt<3>{}, HDInt<l1>{}, HDInt<c1>{}, HDInt<i1>{}, HDInt<j1>{}, HDInt<k1>{});
910  // });
911  // callIfMatch<i1, Is...>([&] {
912  // constexpr auto h0 = HDInt<0>{};
913  // constexpr auto hi1 = HDInt<i1>{};
914  // constexpr auto hc1 = HDInt<c1>{};
915  // constexpr auto hj1 = HDInt<j1>{};
916  // constexpr auto hk1 = HDInt<k1>{};
917  // constexpr auto hl1 = HDInt<l1>{};
918 
919  // func(h0, hi1, hc1, hj1, hk1, hl1);
920  // });
921 
922  // callIfMatch<j1, Is...>([&] {
923  // constexpr auto h1 = HDInt<1>{};
924  // constexpr auto hj1 = HDInt<j1>{};
925  // constexpr auto hc1 = HDInt<c1>{};
926  // constexpr auto hi1 = HDInt<i1>{};
927  // constexpr auto hk1 = HDInt<k1>{};
928  // constexpr auto hl1 = HDInt<l1>{};
929 
930  // func(h1, hj1, hc1, hi1, hk1, hl1);
931  // });
932 
933  // callIfMatch<k1, Is...>([&] {
934  // constexpr auto h2 = HDInt<2>{};
935  // constexpr auto hk1 = HDInt<k1>{};
936  // constexpr auto hc1 = HDInt<c1>{};
937  // constexpr auto hi1 = HDInt<i1>{};
938  // constexpr auto hj1 = HDInt<j1>{};
939  // constexpr auto hl1 = HDInt<l1>{};
940 
941  // func(h2, hk1, hc1, hi1, hj1, hl1);
942  // });
943 
944  // callIfMatch<l1, Is...>([&] {
945  // constexpr auto h3 = HDInt<3>{};
946  // constexpr auto hl1 = HDInt<l1>{};
947  // constexpr auto hc1 = HDInt<c1>{};
948  // constexpr auto hi1 = HDInt<i1>{};
949  // constexpr auto hj1 = HDInt<j1>{};
950  // constexpr auto hk1 = HDInt<k1>{};
951 
952  // func(h3, hl1, hc1, hi1, hj1, hk1);
953  // });
954 
955  }
956  } );
957  }
958  } );
959  } );
960  }
966  template< typename FUNC >
967  GEOS_HOST_DEVICE
968  GEOS_FORCE_INLINE
969  static constexpr void faceBarycentricCoordinateLoop( FUNC && func )
970  {
971  loop< 3 >( std::forward< FUNC >( func ) );
972  }
973 
981  GEOS_HOST_DEVICE
982  GEOS_FORCE_INLINE
983  static
984  void
985  computeReferenceMassMatrix( arraySlice2d< real64 > const & m )
986  {
987  basisLoop( [ &m ] ( auto const cc1, auto const ii1, auto const jj1, auto const kk1, auto const ll1 )
988  {
989  constexpr int c1 = decltype(cc1)::value;
990  constexpr int i1 = decltype(ii1)::value;
991  constexpr int j1 = decltype(jj1)::value;
992  constexpr int k1 = decltype(kk1)::value;
993  constexpr int l1 = decltype(ll1)::value;
994  // Needed for compilers that do not support constexpr lambdas
995  GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
996  basisLoop( [ &m ] ( auto const cc2, auto const ii2, auto const jj2, auto const kk2, auto const ll2 )
997  {
998  constexpr int c2 = decltype(cc2)::value;
999  constexpr int i2 = decltype(ii2)::value;
1000  constexpr int j2 = decltype(jj2)::value;
1001  constexpr int k2 = decltype(kk2)::value;
1002  constexpr int l2 = decltype(ll2)::value;
1003  constexpr real64 val = computeSuperpositionIntegral( i1, j1, k1, l1, i2, j2, k2, l2 );
1004  m[ c1 ][ c2 ] = val;
1005  } );
1006  } );
1007  }
1008 
1020  template< typename DAMPING >
1021  GEOS_HOST_DEVICE
1022  GEOS_FORCE_INLINE
1023  static
1024  void
1025  computeReferenceDampingMatrix( real64 (& d)[numNodes][numNodes], bool const face1Damped, bool const face2Damped, bool const face3Damped, bool const face4Damped )
1026  {
1027 
1028  for( int c1 = 0; c1 < numNodes; c1++ )
1029  {
1030  for( int c2 = 0; c2 < numNodes; c2++ )
1031  {
1032  d[ c1 ][ c2 ] = 0;
1033  }
1034  }
1035  conditionalBasisLoop< 0 >( [&] ( auto const f1, auto const, auto const c1, auto const i1, auto const j1, auto const k1 )
1036  {
1037  conditionalBasisLoop< 0 >( [&] ( auto const f2, auto const, auto const c2, auto const i2, auto const j2, auto const k2 )
1038  {
1039  if constexpr ( f1 == f2 )
1040  {
1041  constexpr real64 val = computeFaceSuperpositionIntegral( i1, j1, k1, i2, j2, k2 );
1042  if( ( f1 == 0 && face1Damped ) ||
1043  ( f1 == 1 && face2Damped ) ||
1044  ( f1 == 2 && face3Damped ) ||
1045  ( f1 == 3 && face4Damped ) )
1046  {
1047  d[ c1 ][ c2 ] += val;
1048  }
1049  }
1050  } );
1051  } );
1052  }
1053 
1060  template< typename FUNC >
1061  GEOS_HOST_DEVICE
1062  GEOS_FORCE_INLINE
1063  static
1064  constexpr
1065  void
1066  computeMassTerm( real64 const (&X)[4][3],
1067  FUNC && func )
1068  {
1069  real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1070  basisLoop( [&func, &detJ] ( auto cc1, auto ii1, auto jj1, auto kk1, auto ll1 )
1071  {
1072  constexpr int c1 = decltype(cc1)::value;
1073  constexpr int i1 = decltype(ii1)::value;
1074  constexpr int j1 = decltype(jj1)::value;
1075  constexpr int k1 = decltype(kk1)::value;
1076  constexpr int l1 = decltype(ll1)::value;
1077  GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
1078 
1079  basisLoop( [&func, &detJ, i1, j1, k1, l1] ( auto c2, auto ii2, auto jj2, auto kk2, auto ll2 )
1080  {
1081  constexpr int c2v = decltype(c2)::value;
1082  constexpr int i2 = decltype(ii2)::value;
1083  constexpr int j2 = decltype(jj2)::value;
1084  constexpr int k2 = decltype(kk2)::value;
1085  constexpr int l2 = decltype(ll2)::value;
1086 
1087  constexpr real64 val = computeSuperpositionIntegral( i1, j1, k1, l1, i2, j2, k2, l2 );
1088  func( c1, c2v, val * detJ );
1089  } );
1090  } );
1091  }
1092 
1102  GEOS_HOST_DEVICE
1103  GEOS_FORCE_INLINE
1104  static constexpr real64
1105  computeFluxDerivativeFactor( real64 const (&X)[4][3], int x1, int x2, int o1, int o2 ) // Order x1, x2, o1, o2
1106  {
1107 
1108  real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1109 
1110  // height with respect to o2
1111 
1112  real64 detJf1 = LvArray::math::abs( faceJacobianDeterminant( o1, X ));
1113 
1114  real64 scal = 1.0;
1115 
1116  if( o1 != o2 )
1117  {
1118 
1119  // scalar product
1120 
1121  real64 detJf2 = LvArray::math::abs( faceJacobianDeterminant( o2, X ));
1122 
1123  real64 el2[3][2] = { { edgeLength2( x1, x2, X ), edgeLength2( o1, o2, X ) },
1124 
1125  { edgeLength2( x1, o1, X ), edgeLength2( x2, o2, X ) },
1126 
1127  { edgeLength2( x1, o2, X ), edgeLength2( x2, o1, X ) } };
1128 
1129  real64 h2 = (el2[1][0]+el2[1][1])-(el2[2][0]+el2[2][1]);
1130 
1131  h2 = (4.0 * el2[0][0]*el2[0][1] - h2 * h2)/16.0;
1132 
1133  scal = (4.0*h2-detJf1 * detJf1 - detJf2 * detJf2 ) / (2.0 * detJf1 * detJf2);
1134 
1135  }
1136 
1137  return scal * detJf1 / detJ;
1138 
1139  }
1140 
1147  template< typename FUNC >
1148  GEOS_HOST_DEVICE
1149  GEOS_FORCE_INLINE
1150  static
1151  constexpr
1152  void
1153  computeStiffnessTerm( real64 const (&X)[4][3],
1154  FUNC && func )
1155  {
1156  real64 detJ = LvArray::math::abs( jacobianDeterminant( X ));
1157  real64 dLambdadX[4][3] = {};
1158  for( int j = 0; j < 3; j++ )
1159  {
1160  dLambdadX[1][j] =
1161  ( ( 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;
1162  dLambdadX[2][j] =
1163  ( ( 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;
1164  dLambdadX[3][j] =
1165  ( ( 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;
1166  dLambdadX[0][j] = -dLambdadX[1][j] - dLambdadX[2][j] - dLambdadX[3][j];
1167  }
1168  basisLoop( [&func, &dLambdadX, &detJ] ( auto const cc1, auto const ci1, auto const cj1, auto const ck1, auto const cl1 )
1169  {
1170 
1171  constexpr int c1 = decltype(cc1)::value;
1172  constexpr int i1 = decltype(ci1)::value;
1173  constexpr int j1 = decltype(cj1)::value;
1174  constexpr int k1 = decltype(ck1)::value;
1175  constexpr int l1 = decltype(cl1)::value;
1176  //Not used in some combinations, but needed for constexpr
1177  GEOS_UNUSED_VAR( c1, i1, j1, k1, l1 );
1178  basisLoop( [&func, &dLambdadX, &detJ] ( auto const cc2, auto const ci2, auto const cj2, auto const ck2, auto const cl2 )
1179  {
1180  constexpr int c2 = decltype(cc2)::value;
1181  constexpr int i2 = decltype(ci2)::value;
1182  constexpr int j2 = decltype(cj2)::value;
1183  constexpr int k2 = decltype(ck2)::value;
1184  constexpr int l2 = decltype(cl2)::value;
1185  //Not used in some combinations, but needed for constexpr
1186  GEOS_UNUSED_VAR( c2, i2, j2, k2, l2 );
1187  barycentricCoordinateLoop( [&func, &dLambdadX, &detJ] ( auto const cd1 )
1188  {
1189  constexpr int d1 = decltype(cd1)::value;
1190  barycentricCoordinateLoop( [&func, &dLambdadX, &detJ] ( auto const cd2 )
1191  {
1192  constexpr int d2 = decltype(cd2)::value;
1193  constexpr int ii1 = i1 + ( d1 == 0 ) * ( -1 );
1194  constexpr int ij1 = j1 + ( d1 == 1 ) * ( -1 );
1195  constexpr int ik1 = k1 + ( d1 == 2 ) * ( -1 );
1196  constexpr int il1 = l1 + ( d1 == 3 ) * ( -1 );
1197  constexpr int ii2 = i2 + ( d2 == 0 ) * ( -1 );
1198  constexpr int ij2 = j2 + ( d2 == 1 ) * ( -1 );
1199  constexpr int ik2 = k2 + ( d2 == 2 ) * ( -1 );
1200  constexpr int il2 = l2 + ( d2 == 3 ) * ( -1 );
1201  constexpr real64 factor1 = correctionFactorDerivative( i1, j1, k1, l1, ii1, ij1, ik1, il1, 3 );
1202  constexpr real64 factor2 = correctionFactorDerivative( i2, j2, k2, l2, ii2, ij2, ik2, il2, 3 );
1203  if constexpr (ii1 >= 0 && ij1 >= 0 && ik1 >= 0 && il1 >= 0 &&
1204  ii2 >= 0 && ij2 >= 0 && ik2 >= 0 && il2 >= 0)
1205  {
1206  constexpr real64 val = computeSuperpositionIntegral( ii1, ij1, ik1, il1, ii2, ij2, ik2, il2 ) * factor1 * factor2;
1207  func( c1, c2, val * detJ * ( dLambdadX[d1][0]*dLambdadX[d2][0] + dLambdadX[d1][1]*dLambdadX[d2][1] + dLambdadX[d1][2]*dLambdadX[d2][2] ) );
1208  }
1209  } );
1210  } );
1211  } );
1212  } );
1213  }
1214 
1215 
1223  GEOS_HOST_DEVICE
1224  GEOS_FORCE_INLINE
1225  static real64 edgeLength2( localIndex i1,
1226  localIndex i2,
1227  real64 const (&X)[4][3] )
1228  {
1229  real64 ab[3] = { X[ i2 ][ 0 ] - X[ i1 ][ 0 ],
1230  X[ i2 ][ 1 ] - X[ i1 ][ 1 ],
1231  X[ i2 ][ 2 ] - X[ i1 ][ 2 ] };
1232  return ab[0] * ab[0] + ab[1] * ab[1]+ ab[2]*ab[2];
1233  }
1234 
1249  template< typename FUNCP, typename FUNCF >
1250  GEOS_HOST_DEVICE
1251  GEOS_FORCE_INLINE
1252  static
1253  constexpr
1254  void
1255  computeSurfaceTerms( real64 const (&X)[4][3],
1256  FUNCP && funcP,
1257  FUNCF && funcF )
1258  {
1259  real64 detJf[4] = { faceJacobianDeterminant( 0, X ), faceJacobianDeterminant( 1, X ),
1260  faceJacobianDeterminant( 2, X ), faceJacobianDeterminant( 3, X ) };
1261  conditionalBasisLoop< 0, 1 >( [&] ( auto const cf1, auto const cd, auto const cc1, auto const ci1, auto const cj1, auto const ck1 )
1262  {
1263  conditionalBasisLoop< 0 >( [&] ( auto const cf2, auto const, auto const cc2, auto const ci2, auto const cj2, auto const ck2 )
1264  {
1265  constexpr int c1 = decltype(cc1)::value;
1266  constexpr int i1 = decltype(ci1)::value;
1267  constexpr int j1 = decltype(cj1)::value;
1268  constexpr int k1 = decltype(ck1)::value;
1269  // Not used in some combinations, but needed for constexpr
1270  GEOS_UNUSED_VAR( c1, i1, j1, k1 );
1271  constexpr int f2 = decltype(cf2)::value;
1272  constexpr int c2 = decltype(cc2)::value;
1273  constexpr int i2 = decltype(ci2)::value;
1274  constexpr int j2 = decltype(cj2)::value;
1275  constexpr int k2 = decltype(ck2)::value;
1276  // Not used in some combinations, but needed for constexpr
1277  GEOS_UNUSED_VAR( c2, i2, j2, k2, f2 );
1278 
1279  // The second function is nonzero on the face indexed by f2, so we integrate on this face.
1280  if constexpr ( std::decay_t< decltype(cf1) >::value == std::decay_t< decltype(cf2) >::value)
1281  {
1282  // compute penalty term iff the other function is also nonzero on the same face (i.e., d1==0)
1283  if constexpr ( std::decay_t< decltype(cd) >::value == 0 )
1284  {
1285  constexpr real64 val = computeFaceSuperpositionIntegral( i1, j1, k1, i2, j2, k2 );
1286  funcP( c1, c2, f2, i1, j1, k1, i2, j2, k2, val * detJf[ f2 ] );
1287  }
1288  // Compute flux term. This is nonzero in two cases.
1289  // first case: function has exponent 1 wrt to the same face. In this case, one can derive it once wrt to the
1290  // corresponding lambda and it will obtain a nonzero function on the face.
1291  if constexpr ( std::decay_t< decltype(cd) >::value == 1 )
1292  {
1293  constexpr real64 derFactor = ( i1 + j1 + k1 + 4 );
1294  constexpr real64 val = computeFaceSuperpositionIntegral( i1, j1, k1, i2, j2, k2 ) * derFactor;
1295  funcF( c1, c2, f2, -1, i1, j1, k1, i2, j2, k2, val * detJf[ f2 ] );
1296  }
1297  // second case: function has exponent zero wrt f2.
1298  // In this case, one can derive it wrt to any other face.
1299  else if constexpr ( std::decay_t< decltype(cd) >::value == 0 )
1300  {
1301  faceBarycentricCoordinateLoop( [ &funcF, &detJf ] ( auto const cl )
1302  {
1303  constexpr int l = decltype(cl)::value;
1304  constexpr int ii1 = i1 + ( l == 0 ) * ( -1 );
1305  constexpr int ij1 = j1 + ( l == 1 ) * ( -1 );
1306  constexpr int ik1 = k1 + ( l == 2 ) * ( -1 );
1307  if constexpr (ii1 >= 0 && ij1 >= 0 && ik1 >= 0)
1308  {
1309  constexpr real64 derFactor = ( ii1 + ij1 + ik1 + 4 );
1310  constexpr real64 val = computeFaceSuperpositionIntegral( ii1, ij1, ik1, i2, j2, k2 ) * derFactor;
1311  constexpr int f = l >= f2 ? l + 1 : l;
1312  funcF( c1, c2, f2, f, i1, j1, k1, i2, j2, k2, val * detJf[f2] );
1313  }
1314  } );
1315  }
1316  }
1317  } );
1318  } );
1319  }
1320 
1321 };
1322 
1326 using BB1_Tetrahedron = BB_Tetrahedron< 1 >;
1330 using BB2_Tetrahedron = BB_Tetrahedron< 2 >;
1334 using BB3_Tetrahedron = BB_Tetrahedron< 3 >;
1335 //using BB4_Tetrahedron = BB_Tetrahedron< 4 >;
1336 //using BB5_Tetrahedron = BB_Tetrahedron< 5 >;
1337 
1338 }
1339 }
1340 
1341 #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:1255
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:1225
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:1025
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:799
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:969
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:1105
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:737
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:677
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:985
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:1153
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:690
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:1066
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