GEOS
LagrangeBasis5GL.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 Total, S.A
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
16 #ifndef GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_ELEMENTFORMULATIONS_LAGRANGEBASIS5GL_HPP_
17 #define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_ELEMENTFORMULATIONS_LAGRANGEBASIS5GL_HPP_
18 
23 #include "common/DataTypes.hpp"
24 
25 namespace geos
26 {
27 namespace finiteElement
28 {
29 
40 {
41 public:
42 
44  constexpr static localIndex numSupportPoints = 6;
45 
47  static constexpr real64 sqrt_7_ = 2.64575131106459059;
48 
50  static constexpr real64 sqrt__7_plus_2sqrt7__ = 3.50592393273573196;
51 
53  static constexpr real64 sqrt__7_mins_2sqrt7__ = 1.30709501485960033;
54 
56  static constexpr real64 sqrt__7_plus_sqrt7_div2__ = 2.884939454396278;
57 
59  static constexpr real64 sqrt__7_mins_sqrt7_div2__ = 2.382671682055189;
60 
62  static constexpr real64 sqrt_inv21 = 0.218217890235992381;
63 
71  constexpr static real64 weight( const int q )
72  {
73  switch( q )
74  {
75  case 1:
76  case 4:
77  return (1.0/30.0)*(14.0-sqrt_7_);
78  case 2:
79  case 3:
80  return (1.0/30.0)*(14.0+sqrt_7_);
81  default:
82  return 1.0/15.0;
83  }
84  }
85 
94  constexpr static real64 parentSupportCoord( const localIndex supportPointIndex )
95  {
96 
97  real64 result=0.0;
98 
99  switch( supportPointIndex )
100  {
101  case 0:
102  result = -1.0;
103  break;
104 
105  case 1:
107  break;
108 
109  case 2:
111  break;
112 
113  case 3:
115  break;
116 
117  case 4:
119  break;
120 
121  case 5:
122  result = 1.0;
123  break;
124 
125  default:
126  break;
127  }
128 
129  return result;
130  }
131 
141  constexpr static real64 value( const int index,
142  const real64 xi )
143  {
144 
145  real64 result=0.0;
146 
147  switch( index )
148  {
149  case 0:
150  return result = LagrangeBasis5GL::value0( xi );
151  break;
152 
153  case 1:
154  return result = LagrangeBasis5GL::value1( xi );
155  break;
156 
157  case 2:
158  return result = LagrangeBasis5GL::value2( xi );
159  break;
160 
161  case 3:
162  return result = LagrangeBasis5GL::value3( xi );
163  break;
164 
165  case 4:
166  return result = LagrangeBasis5GL::value4( xi );
167  break;
168 
169  case 5:
170  return result = LagrangeBasis5GL::value5( xi );
171  break;
172 
173  default:
174  break;
175  }
176 
177  return result;
178 
179  }
180 
181 
189  constexpr static real64 value0( const real64 xi )
190  {
191  /* Define the two GL points needed to compute the basis function at point index 0. Here we need the points
192  at index 3,4 called lambda3, lambda4. */
193 
196 
197  return (-21.0/16.0)*(xi*xi*xi*xi*xi-xi*xi*xi*xi-(lambda3*lambda3+lambda4*lambda4)*xi*xi*xi+(lambda3*lambda3+lambda4*lambda4)*xi*xi+
198  lambda3*lambda3*lambda4*lambda4*xi-lambda3*lambda3*lambda4*lambda4);
199  }
200 
201 
209  constexpr static real64 value1( const real64 xi )
210  {
211  /* Define the two GL points needed to compute the basis function at point index 1. Here we need the points
212  at index 3 and 4 called lambda3, lambda4. */
213 
216 
217  return ((21.0/16.0)*sqrt__7_mins_sqrt7_div2__)*(xi*xi*xi*xi*xi-lambda4*xi*xi*xi*xi-(lambda3*lambda3+1)*xi*xi*xi+lambda4*(lambda3*lambda3+1)*xi*xi+
218  lambda3*lambda3*xi-lambda4*lambda3*lambda3);
219  }
220 
228  constexpr static real64 value2( const real64 xi )
229  {
230  /* Define the two GL points needed to compute the basis function at point index 2. Here we need the points
231  at index 3 and 4 called lambda3, lambda4. */
232 
235 
236  return ((-21.0/16.0)*sqrt__7_plus_sqrt7_div2__)*(xi*xi*xi*xi*xi-lambda3*xi*xi*xi*xi-(lambda4*lambda4+1)*xi*xi*xi+lambda3*(lambda4*lambda4+1)*xi*xi+
237  lambda4*lambda4*xi-lambda3*lambda4*lambda4);
238  }
239 
247  constexpr static real64 value3( const real64 xi )
248  {
249  /* Define the two GL points needed to compute the basis function at point index 3. Here we need the points
250  at index 3 and 4 called lambda1, lambda2. */
251 
254 
255  return ((21.0/16.0)*sqrt__7_plus_sqrt7_div2__)*(xi*xi*xi*xi*xi+lambda3*xi*xi*xi*xi-(lambda4*lambda4+1)*xi*xi*xi-lambda3*(lambda4*lambda4+1)*xi*xi+
256  lambda4*lambda4*xi+lambda3*lambda4*lambda4);
257  }
258 
266  constexpr static real64 value4( const real64 xi )
267  {
268  /* Define the two GL points needed to compute the basis function at point index 4. Here we need the points
269  at index 3 and 4 called lambda3, lambda4. */
270 
273 
274  return ((-21.0/16.0)*sqrt__7_mins_sqrt7_div2__)*(xi*xi*xi*xi*xi+lambda4*xi*xi*xi*xi-(lambda3*lambda3+1)*xi*xi*xi-lambda4*(lambda3*lambda3+1)*xi*xi+
275  lambda3*lambda3*xi+lambda4*lambda3*lambda3);
276  }
277 
285  constexpr static real64 value5( const real64 xi )
286  {
287  /* Define the two GL points needed to compute the basis function at point index 5. Here we need the points
288  at index 3 and 4 called lambda3, lambda4. */
289 
292 
293  return (21.0/16.0)*(xi*xi*xi*xi*xi+xi*xi*xi*xi-(lambda4*lambda4+lambda3*lambda3)*xi*xi*xi-(lambda3*lambda3+lambda4*lambda4)*xi*xi+
294  lambda3*lambda3*lambda4*lambda4*xi+lambda4*lambda4*lambda3*lambda3);
295  }
296 
307  constexpr static real64 gradient( const int index,
308  const real64 xi )
309  {
310 
311  real64 result=0.0;
312 
313  switch( index )
314  {
315  case 0:
316  result = LagrangeBasis5GL::gradient0( xi );
317  break;
318 
319  case 1:
320  result = LagrangeBasis5GL::gradient1( xi );
321  break;
322 
323  case 2:
324  result = LagrangeBasis5GL::gradient2( xi );
325  break;
326 
327  case 3:
328  result = LagrangeBasis5GL::gradient3( xi );
329  break;
330 
331  case 4:
332  result = LagrangeBasis5GL::gradient4( xi );
333  break;
334 
335  case 5:
336  result = LagrangeBasis5GL::gradient5( xi );
337  break;
338 
339  default:
340  break;
341  }
342 
343  return result;
344 
345  }
346 
355  constexpr static real64 gradient0( const real64 xi )
356  {
357 
360 
361  return (-21.0/16.0)*(5.0*xi*xi*xi*xi-4.0*xi*xi*xi-3.0*(lambda3*lambda3+lambda4*lambda4)*xi*xi+2.0*(lambda3*lambda3+lambda4*lambda4)*xi+lambda3*lambda3*lambda4*lambda4);
362 
363  }
364 
373  constexpr static real64 gradient1( const real64 xi )
374  {
375 
378 
379  return (21.0/16.0)*sqrt__7_mins_sqrt7_div2__*(5.0*xi*xi*xi*xi-4.0*lambda4*xi*xi*xi-3.0*(lambda3*lambda3+1.0)*xi*xi+2.0*lambda4*(lambda3*lambda3+1.0)*xi+lambda3*lambda3);
380 
381  }
382 
391  constexpr static real64 gradient2( const real64 xi )
392  {
393 
396 
397  return (-21.0/16.0)*sqrt__7_plus_sqrt7_div2__*(5.0*xi*xi*xi*xi-4.0*lambda3*xi*xi*xi-3.0*(lambda4*lambda4+1.0)*xi*xi+2.0*lambda3*(lambda4*lambda4+1.0)*xi+lambda4*lambda4);
398 
399  }
400 
409  constexpr static real64 gradient3( const real64 xi )
410  {
411 
414 
415  return (21.0/16.0)*sqrt__7_plus_sqrt7_div2__*(5.0*xi*xi*xi*xi+4.0*lambda3*xi*xi*xi-3.0*(lambda4*lambda4+1.0)*xi*xi-2*lambda3*(lambda4*lambda4+1.0)*xi+lambda4*lambda4);
416 
417  }
418 
427  constexpr static real64 gradient4( const real64 xi )
428  {
429 
432 
433  return (-21.0/16.0)*sqrt__7_mins_sqrt7_div2__*(5.0*xi*xi*xi*xi+4.0*lambda4*xi*xi*xi-3.0*(lambda3*lambda3+1.0)*xi*xi-2.0*lambda4*(lambda3*lambda3+1.0)*xi+lambda3*lambda3);
434  }
435 
444  constexpr static real64 gradient5( const real64 xi )
445  {
446 
449 
450  return (21.0/16.0)*(5.0*xi*xi*xi*xi+4.0*xi*xi*xi-3.0*(lambda3*lambda3+lambda4*lambda4)*xi*xi-2.0*(lambda3*lambda3+lambda4*lambda4)*xi+lambda3*lambda3*lambda4*lambda4);
451 
452  }
453 
463  constexpr static real64 gradientAt( const int q,
464  const int p )
465  {
466  switch( q )
467  {
468  case 0:
469  switch( p )
470  {
471  case 0: return -7.5000000000000000000;
472  case 1: return -1.7863649483390948939;
473  case 2: return 0.48495104785356916930;
474  }
475  break;
476  case 1:
477  switch( p )
478  {
479  case 0: return 10.14141593631966928023;
480  case 1: return 0.0;
481  case 2: return -1.72125695283023338321;
482  }
483  break;
484  case 2:
485  switch( p )
486  {
487  case 0: return -4.03618727030534800527;
488  case 1: return 2.5234267774294554319088;
489  case 2: return 0.0;
490  }
491  break;
492  case 3:
493  switch( p )
494  {
495  case 0: return 2.2446846481761668242712;
496  case 1: return -1.1528281585359293413318;
497  case 2: return 1.7529619663678659788775;
498  }
499  break;
500  case 4:
501  switch( p )
502  {
503  case 0: return -1.3499133141904880992312;
504  case 1: return 0.6535475074298001672007;
505  case 2: return -0.7863566722232407374395;
506  }
507  break;
508  case 5:
509  switch( p )
510  {
511  case 0: return 0.500000000000000000000;
512  case 1: return -0.2377811779842313638052;
513  case 2: return 0.2697006108320389724720;
514  }
515  break;
516  }
517  return 0;
518  }
519 
520  /* UNCRUSTIFY-OFF */
521 
554  /* UNCRUSTIFY-ON */
555 
557  {
559  constexpr static localIndex numSupportPoints = 36;
560 
570  constexpr static int linearIndex( const int i,
571  const int j )
572  {
573  return i + 6 * j;
574  }
575 
585  constexpr static void multiIndex( const int linearIndex,
586  int & i0,
587  int & i1 )
588  {
589 
590  i1 = linearIndex/6;
591 
592  i0 = linearIndex%6;
593 
594  }
595 
605  static void value( const real64 (& coords)[2],
606  real64 (& N)[numSupportPoints] )
607  {
608  for( int a=0; a<6; ++a )
609  {
610  for( int b=0; b<6; ++b )
611  {
612  const int lindex = LagrangeBasis5GL::TensorProduct2D::linearIndex( a, b );
613  N[ lindex ] = LagrangeBasis5GL::value( a, coords[0] ) *
614  LagrangeBasis5GL::value( b, coords[1] );
615  }
616  }
617  }
618 
619  };
620 
668  /* UNCRUSTIFY-ON */
669 
671  {
673  constexpr static localIndex numSupportPoints = 216;
674 
685  constexpr static int linearIndex( const int i,
686  const int j,
687  const int k )
688  {
689  return i + 6 * j + 36 * k;
690  }
691 
702  constexpr static void multiIndex( const int linearIndex,
703  int & i0,
704  int & i1,
705  int & i2 )
706  {
707 
708  i2 = linearIndex/36;
709 
710  i1 = (linearIndex%36)/6;
711 
712  i0 = (linearIndex%36)%6;
713 
714  }
715 
725  static void value( const real64 (& coords)[3],
726  real64 (& N)[numSupportPoints] )
727  {
728  for( int a=0; a<6; ++a )
729  {
730  for( int b=0; b<6; ++b )
731  {
732  for( int c=0; c<6; ++c )
733  {
734  const int lindex = LagrangeBasis5GL::TensorProduct3D::linearIndex( a, b, c );
735  N[ lindex ] = LagrangeBasis5GL::value( a, coords[0] ) *
736  LagrangeBasis5GL::value( b, coords[1] ) *
737  LagrangeBasis5GL::value( c, coords[2] );
738  }
739  }
740  }
741  }
742 
743  };
744 
745 };
746 
747 }
748 }
749 
750 
751 #endif /* GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_ELEMENTFORMULATIONS_LAGRANGEBASIS5GL_HPP_ */
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
#define GEOS_FORCE_INLINE
Marks a function or lambda for inlining.
Definition: GeosxMacros.hpp:51
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 value4(const real64 xi)
The value of the basis function for the 4 support point.
static constexpr real64 sqrt__7_mins_2sqrt7__
sqrt( 7 - 2 * sqrt(7) )
static constexpr real64 sqrt__7_plus_2sqrt7__
sqrt( 7 + 2 * sqrt(7) )
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 value(const int index, const real64 xi)
The value of the basis function for a support point evaluated at a point along the axes.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 gradientAt(const int q, const int p)
The gradient of the basis function for a support point evaluated at a given support point....
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 gradient5(const real64 xi)
The gradient of the basis function for support point 5 evaluated at a point along the axes.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 gradient1(const real64 xi)
The gradient of the basis function for support point 1 evaluated at a point along the axes.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 weight(const int q)
The value of the weight for the given support point.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 gradient0(const real64 xi)
The gradient of the basis function for support point 0 evaluated at a point along the axes.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 value3(const real64 xi)
The value of the basis function for the 3 support point.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 gradient3(const real64 xi)
The gradient of the basis function for support point 3 evaluated at a point along the axes.
static constexpr real64 sqrt__7_mins_sqrt7_div2__
sqrt( 7 - 2 * sqrt(7) )
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 parentSupportCoord(const localIndex supportPointIndex)
Calculate the parent coordinates for the xi0 direction, given the linear index of a support point....
static constexpr real64 sqrt__7_plus_sqrt7_div2__
sqrt( 7 + 2 * sqrt(7) )
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 gradient(const int index, const real64 xi)
The gradient of the basis function for a support point evaluated at a point along the axes.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 value2(const real64 xi)
The value of the basis function for the 2 support point.
static constexpr real64 sqrt_inv21
sqrt(1/21)
constexpr static localIndex numSupportPoints
The number of support points for the basis.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 gradient4(const real64 xi)
The gradient of the basis function for support point 4 evaluated at a point along the axes.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 value5(const real64 xi)
The value of the basis function for the 5 support point.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 gradient2(const real64 xi)
The gradient of the basis function for support point 2 evaluated at a point along the axes.
static constexpr real64 sqrt_7_
sqrt(7)
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 value1(const real64 xi)
The value of the basis function for the 1 support point.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE real64 value0(const real64 xi)
The value of the basis function for the 0 support point.
double real64
64-bit floating point type.
Definition: DataTypes.hpp:99
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:85
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE int linearIndex(const int i, const int j)
Calculates the linear index for support/quadrature points from ij coordinates.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE void multiIndex(const int linearIndex, int &i0, int &i1)
Calculate the Cartesian/TensorProduct index given the linear index of a support point.
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void value(const real64(&coords)[2], real64(&N)[numSupportPoints])
The value of the basis function for a support point evaluated at a point along the axes.
constexpr static localIndex numSupportPoints
The number of support points in the basis.
constexpr static localIndex numSupportPoints
The number of support points in the basis.
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void value(const real64(&coords)[3], real64(&N)[numSupportPoints])
The value of the basis function for a support point evaluated at a point along the axes.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE int linearIndex(const int i, const int j, const int k)
Calculates the linear index for support/quadrature points from ijk coordinates.
GEOS_HOST_DEVICE constexpr static GEOS_FORCE_INLINE void multiIndex(const int linearIndex, int &i0, int &i1, int &i2)
Calculate the Cartesian/TensorProduct index given the linear index of a support point.