xref: /libCEED/examples/fluids/qfunctions/utils.h (revision dc9b5c4a033c329154210fb4ea036e15fdbd7186)
113fa47b2SJames Wright // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
213fa47b2SJames Wright // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
313fa47b2SJames Wright //
413fa47b2SJames Wright // SPDX-License-Identifier: BSD-2-Clause
513fa47b2SJames Wright //
613fa47b2SJames Wright // This file is part of CEED:  http://github.com/ceed
713fa47b2SJames Wright 
813fa47b2SJames Wright #ifndef utils_h
913fa47b2SJames Wright #define utils_h
1013fa47b2SJames Wright 
1113fa47b2SJames Wright #include <ceed.h>
12c9c2c079SJeremy L Thompson #include <math.h>
1313fa47b2SJames Wright 
1413fa47b2SJames Wright #ifndef M_PI
1513fa47b2SJames Wright #define M_PI 3.14159265358979323846
1613fa47b2SJames Wright #endif
1713fa47b2SJames Wright 
1813fa47b2SJames Wright CEED_QFUNCTION_HELPER CeedScalar Max(CeedScalar a, CeedScalar b) { return a < b ? b : a; }
1913fa47b2SJames Wright CEED_QFUNCTION_HELPER CeedScalar Min(CeedScalar a, CeedScalar b) { return a < b ? a : b; }
2013fa47b2SJames Wright 
21*dc9b5c4aSJames Wright CEED_QFUNCTION_HELPER void SwapScalar(CeedScalar *a, CeedScalar *b) {
22*dc9b5c4aSJames Wright   CeedScalar temp = *a;
23*dc9b5c4aSJames Wright   *a              = *b;
24*dc9b5c4aSJames Wright   *b              = temp;
25*dc9b5c4aSJames Wright }
26*dc9b5c4aSJames Wright 
2713fa47b2SJames Wright CEED_QFUNCTION_HELPER CeedScalar Square(CeedScalar x) { return x * x; }
2813fa47b2SJames Wright CEED_QFUNCTION_HELPER CeedScalar Cube(CeedScalar x) { return x * x * x; }
2913fa47b2SJames Wright 
30530ad8c4SKenneth E. Jansen // @brief Scale vector of length N by scalar alpha
31530ad8c4SKenneth E. Jansen CEED_QFUNCTION_HELPER void ScaleN(CeedScalar *u, const CeedScalar alpha, const CeedInt N) {
32530ad8c4SKenneth E. Jansen   CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) { u[i] *= alpha; }
33530ad8c4SKenneth E. Jansen }
34530ad8c4SKenneth E. Jansen 
3513fa47b2SJames Wright // @brief Dot product of 3 element vectors
362b730f8bSJeremy L Thompson CEED_QFUNCTION_HELPER CeedScalar Dot3(const CeedScalar u[3], const CeedScalar v[3]) { return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]; }
3713fa47b2SJames Wright 
3813fa47b2SJames Wright // @brief Unpack Kelvin-Mandel notation symmetric tensor into full tensor
3913fa47b2SJames Wright CEED_QFUNCTION_HELPER void KMUnpack(const CeedScalar v[6], CeedScalar A[3][3]) {
4013fa47b2SJames Wright   const CeedScalar weight = 1 / sqrt(2.);
4113fa47b2SJames Wright   A[0][0]                 = v[0];
4213fa47b2SJames Wright   A[1][1]                 = v[1];
4313fa47b2SJames Wright   A[2][2]                 = v[2];
4413fa47b2SJames Wright   A[2][1] = A[1][2] = weight * v[3];
4513fa47b2SJames Wright   A[2][0] = A[0][2] = weight * v[4];
4613fa47b2SJames Wright   A[1][0] = A[0][1] = weight * v[5];
4713fa47b2SJames Wright }
4813fa47b2SJames Wright 
49530ad8c4SKenneth E. Jansen // @brief Linear ramp evaluation
50530ad8c4SKenneth E. Jansen CEED_QFUNCTION_HELPER CeedScalar LinearRampCoefficient(CeedScalar amplitude, CeedScalar length, CeedScalar start, CeedScalar x) {
51530ad8c4SKenneth E. Jansen   if (x < start) {
52530ad8c4SKenneth E. Jansen     return amplitude;
53530ad8c4SKenneth E. Jansen   } else if (x < start + length) {
54530ad8c4SKenneth E. Jansen     return amplitude * ((x - start) * (-1 / length) + 1);
55530ad8c4SKenneth E. Jansen   } else {
56530ad8c4SKenneth E. Jansen     return 0;
57530ad8c4SKenneth E. Jansen   }
58530ad8c4SKenneth E. Jansen }
59530ad8c4SKenneth E. Jansen 
6013fa47b2SJames Wright #endif  // utils_h
61