GEOS
H1_Tetrahedron_Lagrange1_Gauss1.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 Total, S.A
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
20 #ifndef GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS1_HPP_
21 #define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS1_HPP_
22 
23 #include "FiniteElementBase.hpp"
24 #include "LagrangeBasis1.hpp"
25 
26 
27 namespace geos
28 {
29 namespace finiteElement
30 {
31 
50 {
51 public:
52 
55 
57  constexpr static localIndex numNodes = 4;
58 
60  constexpr static localIndex numFaces = 4;
61 
63  constexpr static localIndex maxSupportPoints = numNodes;
64 
66  constexpr static localIndex numQuadraturePoints = 1;
67 
70 
72  virtual ~H1_Tetrahedron_Lagrange1_Gauss1() override
73  {}
74 
76  virtual localIndex getNumQuadraturePoints() const override
77  {
78  return numQuadraturePoints;
79  }
80 
87  static localIndex
89  {
90  GEOS_UNUSED_VAR( stack );
91  return numQuadraturePoints;
92  }
93 
95  virtual localIndex getNumSupportPoints() const override
96  {
97  return numNodes;
98  }
99 
101  virtual localIndex getMaxSupportPoints() const override
102  {
103  return maxSupportPoints;
104  }
105 
113  {
114  GEOS_UNUSED_VAR( stack );
115  return numNodes;
116  }
117 
126  static void getSamplingPointCoordInParentSpace( int const & linearIndex,
127  real64 (& samplingPointCoord)[3] )
128  {
129  int const i0 = linearIndex % numSamplingPointsPerDirection;
130  int const i1 = ( (linearIndex - i0)/numSamplingPointsPerDirection ) % numSamplingPointsPerDirection;
131  int const i2 = ( (linearIndex - i0)/numSamplingPointsPerDirection - i1 ) / numSamplingPointsPerDirection;
132 
133  real64 const step = 1. / ( numSamplingPointsPerDirection - 1 );
134 
135  real64 const r = i0 * step;
136  real64 const s = i1 * step;
137  real64 const t = i2 * step;
138  if( (r+s) <= t )
139  {
140  samplingPointCoord[0] = r;
141  samplingPointCoord[1] = s;
142  samplingPointCoord[2] = t;
143  }
144  else
145  {
146  // if outside of the triangle need to reproject it. Points will be doubled though.
147  samplingPointCoord[0] = 1 - t - r;
148  samplingPointCoord[1] = 1 - t - s;
149  samplingPointCoord[2] = t;
150  }
151  }
152 
162  static void calcN( real64 const (&pointCoord)[3],
163  real64 ( &N )[numNodes] );
164 
173  static void calcN( localIndex const q,
174  real64 ( &N )[numNodes] );
175 
185  inline
186  static void calcN( localIndex const q,
187  StackVariables const & stack,
188  real64 ( &N )[numNodes] );
189 
199  static void calcFaceBubbleN( real64 const (&pointCoord)[3],
200  real64 (& N)[numFaces] )
201  {
202 
203  real64 const r = pointCoord[0];
204  real64 const s = pointCoord[1];
205  real64 const t = pointCoord[2];
206 
207  N[0] = (1 - r - s - t) * r * t;
208  N[1] = (1 - r - s - t) * r * s;
209  N[2] = (1 - r - s - t) * t * s;
210  N[3] = r * s * t;
211  }
212 
221  inline
222  static void calcFaceBubbleN( localIndex const q,
223  real64 (& N)[numFaces] )
224  {
225  GEOS_UNUSED_VAR( q );
226 
227  // single quadrature point (centroid), i.e. r = s = t = 1/4
228  real64 const pointCoord[3] = {0.25, 0.25, 0.25};
229 
230  calcFaceBubbleN( pointCoord, N );
231  }
232 
243  static real64 calcGradN( localIndex const q,
244  real64 const (&X)[numNodes][3],
245  real64 ( &gradN )[numNodes][3] );
246 
258  inline
259  static real64 calcGradN( localIndex const q,
260  real64 const (&X)[numNodes][3],
261  StackVariables const & stack,
262  real64 ( &gradN )[numNodes][3] );
263 
275  real64 const (&X)[numNodes][3],
276  real64 ( &gradN )[numFaces][3] );
277 
287  real64 const (&X)[numNodes][3] );
288 
299  real64 const (&X)[numNodes][3],
300  StackVariables const & stack )
301  { GEOS_UNUSED_VAR( stack ); return transformedQuadratureWeight( q, X ); }
302 
312  static real64 invJacobianTransformation( int const q,
313  real64 const (&X)[numNodes][3],
314  real64 ( & J )[3][3] )
315  {
316  GEOS_UNUSED_VAR( q, X );
317 // jacobianTransformation( q, X, J );
318  return LvArray::tensorOps::invert< 3 >( J );
319  }
320 private:
322  constexpr static real64 parentVolume = 1.0 / 6.0;
323 
325  constexpr static real64 weight = parentVolume / numQuadraturePoints;
326 
334  static real64 determinantJacobianTransformation( real64 const (&X)[numNodes][3] );
335 
336 };
337 
339 
341 inline
342 real64
343 H1_Tetrahedron_Lagrange1_Gauss1::
344  determinantJacobianTransformation( real64 const (&X)[numNodes][3] )
345 {
346  return ( X[1][0] - X[0][0] )*( ( X[2][1] - X[0][1] )*( X[3][2] - X[0][2] ) - ( X[3][1] - X[0][1] )*( X[2][2] - X[0][2] ) )
347  + ( X[1][1] - X[0][1] )*( ( X[3][0] - X[0][0] )*( X[2][2] - X[0][2] ) - ( X[2][0] - X[0][0] )*( X[3][2] - X[0][2] ) )
348  + ( X[1][2] - X[0][2] )*( ( X[2][0] - X[0][0] )*( X[3][1] - X[0][1] ) - ( X[3][0] - X[0][0] )*( X[2][1] - X[0][1] ) );
349 }
350 
351 //*************************************************************************************************
352 
354 inline
355 void
357  calcN( real64 const (&coords)[3],
358  real64 (& N)[numNodes] )
359 {
360  real64 const r = coords[0];
361  real64 const s = coords[1];
362  real64 const t = coords[2];
363 
364  // single quadrature point (centroid), i.e. r = s = t = 1/4
365  N[0] = 1 - r - s - t;
366  N[1] = r;
367  N[2] = s;
368  N[3] = t;
369 }
370 
371 
373 inline
374 void
376  calcN( localIndex const q,
377  real64 (& N)[numNodes] )
378 {
379  GEOS_UNUSED_VAR( q );
380 
381  // single quadrature point (centroid), i.e. r = s = t = 1/4
382  real64 const pointCoord[3] = {0.25, 0.25, 0.25};
383  calcN( pointCoord, N );
384 }
385 
387 inline
389  calcN( localIndex const q,
390  StackVariables const & GEOS_UNUSED_PARAM( stack ),
391  real64 ( & N )[numNodes] )
392 {
393  return calcN( q, N );
394 }
395 
396 //*************************************************************************************************
397 
399 inline
400 real64
402  calcGradN( localIndex const q,
403  real64 const (&X)[numNodes][3],
404  real64 (& gradN)[numNodes][3] )
405 {
406  GEOS_UNUSED_VAR( q );
407 
408  gradN[0][0] = X[1][1]*( X[3][2] - X[2][2] ) - X[2][1]*( X[3][2] - X[1][2] ) + X[3][1]*( X[2][2] - X[1][2] );
409  gradN[0][1] = -X[1][0]*( X[3][2] - X[2][2] ) + X[2][0]*( X[3][2] - X[1][2] ) - X[3][0]*( X[2][2] - X[1][2] );
410  gradN[0][2] = X[1][0]*( X[3][1] - X[2][1] ) - X[2][0]*( X[3][1] - X[1][1] ) + X[3][0]*( X[2][1] - X[1][1] );
411 
412  gradN[1][0] = -X[0][1]*( X[3][2] - X[2][2] ) + X[2][1]*( X[3][2] - X[0][2] ) - X[3][1]*( X[2][2] - X[0][2] );
413  gradN[1][1] = X[0][0]*( X[3][2] - X[2][2] ) - X[2][0]*( X[3][2] - X[0][2] ) + X[3][0]*( X[2][2] - X[0][2] );
414  gradN[1][2] = -X[0][0]*( X[3][1] - X[2][1] ) + X[2][0]*( X[3][1] - X[0][1] ) - X[3][0]*( X[2][1] - X[0][1] );
415 
416  gradN[2][0] = X[0][1]*( X[3][2] - X[1][2] ) - X[1][1]*( X[3][2] - X[0][2] ) + X[3][1]*( X[1][2] - X[0][2] );
417  gradN[2][1] = -X[0][0]*( X[3][2] - X[1][2] ) + X[1][0]*( X[3][2] - X[0][2] ) - X[3][0]*( X[1][2] - X[0][2] );
418  gradN[2][2] = X[0][0]*( X[3][1] - X[1][1] ) - X[1][0]*( X[3][1] - X[0][1] ) + X[3][0]*( X[1][1] - X[0][1] );
419 
420  gradN[3][0] = -X[0][1]*( X[2][2] - X[1][2] ) + X[1][1]*( X[2][2] - X[0][2] ) - X[2][1]*( X[1][2] - X[0][2] );
421  gradN[3][1] = X[0][0]*( X[2][2] - X[1][2] ) - X[1][0]*( X[2][2] - X[0][2] ) + X[2][0]*( X[1][2] - X[0][2] );
422  gradN[3][2] = -X[0][0]*( X[2][1] - X[1][1] ) + X[1][0]*( X[2][1] - X[0][1] ) - X[2][0]*( X[1][1] - X[0][1] );
423 
424  real64 detJ = determinantJacobianTransformation( X );
425  real64 factor = 1.0 / ( detJ );
426 
427  for( int i = 0; i < numNodes; ++i )
428  {
429  for( int j = 0; j < 3; ++j )
430  {
431  gradN[i][j] *= factor;
432  }
433  }
434 
435  return detJ * weight;
436 }
437 
439 inline
441  calcGradN( localIndex const q,
442  real64 const (&X)[numNodes][3],
443  StackVariables const & GEOS_UNUSED_PARAM( stack ),
444  real64 ( & gradN )[numNodes][3] )
445 {
446  return calcGradN( q, X, gradN );
447 }
448 
450 inline
451 real64
453  real64 const (&X)[numNodes][3],
454  real64 (& gradN)[numFaces][3] )
455 {
456 
457  GEOS_UNUSED_VAR( q );
458 
459  real64 detJ = determinantJacobianTransformation( X );
460  real64 factor = 1.0 / ( detJ );
461 
462  real64 J[3][3] = {{0}};
463 
464  J[0][0] = (-X[0][1]*( X[3][2] - X[2][2] ) + X[2][1]*( X[3][2] - X[0][2] ) - X[3][1]*( X[2][2] - X[0][2] ))*factor;
465  J[0][1] = ( X[0][0]*( X[3][2] - X[2][2] ) - X[2][0]*( X[3][2] - X[0][2] ) + X[3][0]*( X[2][2] - X[0][2] ))*factor;
466  J[0][2] = (-X[0][0]*( X[3][1] - X[2][1] ) + X[2][0]*( X[3][1] - X[0][1] ) - X[3][0]*( X[2][1] - X[0][1] ))*factor;
467 
468  J[1][0] = ( X[0][1]*( X[3][2] - X[1][2] ) - X[1][1]*( X[3][2] - X[0][2] ) + X[3][1]*( X[1][2] - X[0][2] ))*factor;
469  J[1][1] = (-X[0][0]*( X[3][2] - X[1][2] ) + X[1][0]*( X[3][2] - X[0][2] ) - X[3][0]*( X[1][2] - X[0][2] ))*factor;
470  J[1][2] = ( X[0][0]*( X[3][1] - X[1][1] ) - X[1][0]*( X[3][1] - X[0][1] ) + X[3][0]*( X[1][1] - X[0][1] ))*factor;
471 
472  J[2][0] = (-X[0][1]*( X[2][2] - X[1][2] ) + X[1][1]*( X[2][2] - X[0][2] ) - X[2][1]*( X[1][2] - X[0][2] ))*factor;
473  J[2][1] = ( X[0][0]*( X[2][2] - X[1][2] ) - X[1][0]*( X[2][2] - X[0][2] ) + X[2][0]*( X[1][2] - X[0][2] ))*factor;
474  J[2][2] = (-X[0][0]*( X[2][1] - X[1][1] ) + X[1][0]*( X[2][1] - X[0][1] ) - X[2][0]*( X[1][1] - X[0][1] ))*factor;
475 
476  real64 dNdXi[numFaces][3] = {{0}};
477  // single quadrature point (centroid), i.e. r = s = t = 1/4
478  real64 const r = 1.0 / 4.0;
479  real64 const s = 1.0 / 4.0;
480  real64 const t = 1.0 / 4.0;
481 
482  dNdXi[0][0] = ( 1 - 2 * r - s - t ) * t; // dN0/dr
483  dNdXi[0][1] = -r * t; // dN0/ds
484  dNdXi[0][2] = ( 1 - r - s - 2 * t ) * r; // dN0/dt
485 
486  dNdXi[1][0] = ( 1 - 2 * r - s - t ) * s; // dN1/dr
487  dNdXi[1][1] = ( 1 - r - 2 * s - t ) * r; // dN1/ds
488  dNdXi[1][2] = -r * s; // dN1/dt
489 
490  dNdXi[2][0] = -t * s; // dN2/dr
491  dNdXi[2][1] = ( 1 - r - 2 * s - t ) * t; // dN2/ds
492  dNdXi[2][2] = ( 1 - r - s - 2 * t ) * s; // dN2/dt
493 
494  dNdXi[3][0] = t * s; // dN3/dr
495  dNdXi[3][1] = r * t; // dN3/ds
496  dNdXi[3][2] = r * s; // dN3/dt
497 
498  for( int fi=0; fi<numFaces; ++fi )
499  {
500  gradN[fi][0] = dNdXi[fi][0] * J[0][0] + dNdXi[fi][1] * J[1][0] + dNdXi[fi][2] * J[2][0];
501  gradN[fi][1] = dNdXi[fi][0] * J[0][1] + dNdXi[fi][1] * J[1][1] + dNdXi[fi][2] * J[2][1];
502  gradN[fi][2] = dNdXi[fi][0] * J[0][2] + dNdXi[fi][1] * J[1][2] + dNdXi[fi][2] * J[2][2];
503  }
504 
505  return detJ * weight;
506 }
507 
508 //*************************************************************************************************
509 
511 inline
512 real64
515  real64 const (&X)[numNodes][3] )
516 {
517  GEOS_UNUSED_VAR( q );
518 
519  real64 detJ = determinantJacobianTransformation( X );
520 
521  return detJ * weight;
522 }
523 
525 
526 }
527 }
528 
529 #endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS1_HPP_
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
#define GEOS_FORCE_INLINE
Marks a function or lambda for inlining.
Definition: GeosxMacros.hpp:51
#define GEOS_UNUSED_PARAM(X)
Mark an unused argument and silence compiler warnings.
Definition: GeosxMacros.hpp:72
Base class for FEM element implementations.
constexpr static int numSamplingPointsPerDirection
Number of sampling points.
static GEOS_HOST_DEVICE 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.
static GEOS_HOST_DEVICE localIndex getNumSupportPoints(StackVariables const &stack)
Get the number of support points.
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void getSamplingPointCoordInParentSpace(int const &linearIndex, real64(&samplingPointCoord)[3])
Get the Sampling Point Coord In the Parent Space.
static GEOS_HOST_DEVICE real64 transformedQuadratureWeight(localIndex const q, real64 const (&X)[numNodes][3])
Calculate the integration weights for a quadrature point.
static GEOS_HOST_DEVICE void calcN(localIndex const q, real64(&N)[numNodes])
Calculate shape functions values for each support point at a quadrature point.
constexpr static localIndex numFaces
The number of faces/support points per element.
constexpr static localIndex numNodes
The number of nodes/support points per element.
constexpr static int numSamplingPoints
The number of sampling points per element.
static GEOS_HOST_DEVICE real64 calcGradFaceBubbleN(localIndex const q, real64 const (&X)[numNodes][3], real64(&gradN)[numFaces][3])
Calculate the bubble function derivatives wrt the physical coordinates.
static GEOS_HOST_DEVICE void calcFaceBubbleN(localIndex const q, real64(&N)[numFaces])
Calculate face bubble functions values for each face at a quadrature point.
constexpr static localIndex numQuadraturePoints
The number of quadrature points per element.
constexpr static localIndex maxSupportPoints
The maximum number of support points per element.
static GEOS_HOST_DEVICE real64 transformedQuadratureWeight(localIndex const q, real64 const (&X)[numNodes][3], StackVariables const &stack)
Calculate the integration weights for a quadrature point.
virtual GEOS_HOST_DEVICE localIndex getNumQuadraturePoints() const override
Virtual getter for the number of quadrature points per element.
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcFaceBubbleN(real64 const (&pointCoord)[3], real64(&N)[numFaces])
Calculate face bubble functions values for each face at a given point in the parent space.
static GEOS_HOST_DEVICE real64 calcGradN(localIndex const q, real64 const (&X)[numNodes][3], real64(&gradN)[numNodes][3])
Calculate the shape functions derivatives wrt the physical coordinates.
virtual GEOS_HOST_DEVICE localIndex getNumSupportPoints() const override
Virtual getter for the number of support points per element.
GEOS_HOST_DEVICE static GEOS_FORCE_INLINE void calcN(real64 const (&pointCoord)[3], real64(&N)[numNodes])
Calculate shape functions values for each support point at a given point in the parent space.
static GEOS_HOST_DEVICE localIndex getNumQuadraturePoints(StackVariables const &stack)
Get the number of quadrature points.
virtual GEOS_HOST_DEVICE localIndex getMaxSupportPoints() const override
Get the maximum number of support points for this element.
static GEOS_HOST_DEVICE real64 invJacobianTransformation(int const q, real64 const (&X)[numNodes][3], real64(&J)[3][3])
Calculates the isoparametric "Jacobian" transformation matrix/mapping from the parent space to the ph...
static GEOS_HOST_DEVICE void calcN(localIndex const q, StackVariables const &stack, real64(&N)[numNodes])
Calculate shape functions values for each support point at a quadrature 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