1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7 #pragma once
8
9 #include <ceed/types.h>
10 #ifndef CEED_RUNNING_JIT_PASS
11 #include <math.h>
12 #endif
13
14 #ifndef M_PI
15 #define M_PI 3.14159265358979323846
16 #endif
17
Max(CeedScalar a,CeedScalar b)18 CEED_QFUNCTION_HELPER CeedScalar Max(CeedScalar a, CeedScalar b) { return a < b ? b : a; }
Min(CeedScalar a,CeedScalar b)19 CEED_QFUNCTION_HELPER CeedScalar Min(CeedScalar a, CeedScalar b) { return a < b ? a : b; }
20
SwapScalar(CeedScalar * a,CeedScalar * b)21 CEED_QFUNCTION_HELPER void SwapScalar(CeedScalar *a, CeedScalar *b) {
22 CeedScalar temp = *a;
23 *a = *b;
24 *b = temp;
25 }
26
Square(CeedScalar x)27 CEED_QFUNCTION_HELPER CeedScalar Square(CeedScalar x) { return x * x; }
Cube(CeedScalar x)28 CEED_QFUNCTION_HELPER CeedScalar Cube(CeedScalar x) { return x * x * x; }
29
30 // @brief Scale vector of length N by scalar alpha
ScaleN(CeedScalar * u,const CeedScalar alpha,const CeedInt N)31 CEED_QFUNCTION_HELPER void ScaleN(CeedScalar *u, const CeedScalar alpha, const CeedInt N) {
32 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) u[i] *= alpha;
33 }
34
35 // @brief Set vector of length N to a value alpha
SetValueN(CeedScalar * u,const CeedScalar alpha,const CeedInt N)36 CEED_QFUNCTION_HELPER void SetValueN(CeedScalar *u, const CeedScalar alpha, const CeedInt N) {
37 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) u[i] = alpha;
38 }
39
40 // @brief Copy N elements from x to y
CopyN(const CeedScalar * x,CeedScalar * y,const CeedInt N)41 CEED_QFUNCTION_HELPER void CopyN(const CeedScalar *x, CeedScalar *y, const CeedInt N) { CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) y[i] = x[i]; }
42
43 // @brief Copy 3x3 matrix from A to B
CopyMat3(const CeedScalar A[3][3],CeedScalar B[3][3])44 CEED_QFUNCTION_HELPER void CopyMat3(const CeedScalar A[3][3], CeedScalar B[3][3]) { CopyN((const CeedScalar *)A, (CeedScalar *)B, 9); }
45
46 // @brief Dot product of vectors with N elements
DotN(const CeedScalar * u,const CeedScalar * v,const CeedInt N)47 CEED_QFUNCTION_HELPER CeedScalar DotN(const CeedScalar *u, const CeedScalar *v, const CeedInt N) {
48 CeedScalar output = 0;
49 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) output += u[i] * v[i];
50 return output;
51 }
52
53 // @brief Dot product of 3 element vectors
Dot3(const CeedScalar * u,const CeedScalar * v)54 CEED_QFUNCTION_HELPER CeedScalar Dot3(const CeedScalar *u, const CeedScalar *v) { return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]; }
55
56 // @brief Cross product of vectors with 3 elements
Cross3(const CeedScalar u[3],const CeedScalar v[3],CeedScalar w[3])57 CEED_QFUNCTION_HELPER void Cross3(const CeedScalar u[3], const CeedScalar v[3], CeedScalar w[3]) {
58 w[0] = (u[1] * v[2]) - (u[2] * v[1]);
59 w[1] = (u[2] * v[0]) - (u[0] * v[2]);
60 w[2] = (u[0] * v[1]) - (u[1] * v[0]);
61 }
62
63 // @brief Curl of vector given its gradient
Curl3(const CeedScalar gradient[3][3],CeedScalar v[3])64 CEED_QFUNCTION_HELPER void Curl3(const CeedScalar gradient[3][3], CeedScalar v[3]) {
65 v[0] = gradient[2][1] - gradient[1][2];
66 v[1] = gradient[0][2] - gradient[2][0];
67 v[2] = gradient[1][0] - gradient[0][1];
68 }
69
70 // @brief Matrix vector product, b = Ax + b. A is NxM, x is M, b is N
MatVecNM(const CeedScalar * A,const CeedScalar * x,const CeedInt N,const CeedInt M,const CeedTransposeMode transpose_A,CeedScalar * b)71 CEED_QFUNCTION_HELPER void MatVecNM(const CeedScalar *A, const CeedScalar *x, const CeedInt N, const CeedInt M, const CeedTransposeMode transpose_A,
72 CeedScalar *b) {
73 switch (transpose_A) {
74 case CEED_NOTRANSPOSE:
75 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) b[i] += DotN(&A[i * M], x, M);
76 break;
77 case CEED_TRANSPOSE:
78 CeedPragmaSIMD for (CeedInt i = 0; i < M; i++) { CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) b[i] += A[j * M + i] * x[j]; }
79 break;
80 }
81 }
82
83 // @brief 3x3 Matrix vector product b = Ax + b.
MatVec3(const CeedScalar A[3][3],const CeedScalar x[3],const CeedTransposeMode transpose_A,CeedScalar b[3])84 CEED_QFUNCTION_HELPER void MatVec3(const CeedScalar A[3][3], const CeedScalar x[3], const CeedTransposeMode transpose_A, CeedScalar b[3]) {
85 MatVecNM((const CeedScalar *)A, (const CeedScalar *)x, 3, 3, transpose_A, (CeedScalar *)b);
86 }
87
88 // @brief Matrix-Matrix product, B = DA + B, where D is diagonal.
89 // @details A is NxM, D is diagonal NxN, represented by a vector of length N, and B is NxM. Optionally, A may be transposed.
MatDiagNM(const CeedScalar * A,const CeedScalar * D,const CeedInt N,const CeedInt M,const CeedTransposeMode transpose_A,CeedScalar * B)90 CEED_QFUNCTION_HELPER void MatDiagNM(const CeedScalar *A, const CeedScalar *D, const CeedInt N, const CeedInt M, const CeedTransposeMode transpose_A,
91 CeedScalar *B) {
92 switch (transpose_A) {
93 case CEED_NOTRANSPOSE:
94 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) { CeedPragmaSIMD for (CeedInt j = 0; j < M; j++) B[i * M + j] += D[i] * A[i * M + j]; }
95 break;
96 case CEED_TRANSPOSE:
97 CeedPragmaSIMD for (CeedInt i = 0; i < M; i++) { CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) B[i * N + j] += D[i] * A[j * M + i]; }
98 break;
99 }
100 }
101
102 // @brief 3x3 Matrix-Matrix product, B = DA + B, where D is diagonal.
103 // @details Optionally, A may be transposed.
MatDiag3(const CeedScalar A[3][3],const CeedScalar D[3],const CeedTransposeMode transpose_A,CeedScalar B[3][3])104 CEED_QFUNCTION_HELPER void MatDiag3(const CeedScalar A[3][3], const CeedScalar D[3], const CeedTransposeMode transpose_A, CeedScalar B[3][3]) {
105 MatDiagNM((const CeedScalar *)A, (const CeedScalar *)D, 3, 3, transpose_A, (CeedScalar *)B);
106 }
107 // @brief NxN Matrix-Matrix product, C = AB + C
MatMatN(const CeedScalar * A,const CeedScalar * B,const CeedInt N,const CeedTransposeMode transpose_A,const CeedTransposeMode transpose_B,CeedScalar * C)108 CEED_QFUNCTION_HELPER void MatMatN(const CeedScalar *A, const CeedScalar *B, const CeedInt N, const CeedTransposeMode transpose_A,
109 const CeedTransposeMode transpose_B, CeedScalar *C) {
110 switch (transpose_A) {
111 case CEED_NOTRANSPOSE:
112 switch (transpose_B) {
113 case CEED_NOTRANSPOSE:
114 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) {
115 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) {
116 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[i * N + k] * B[k * N + j];
117 }
118 }
119 break;
120 case CEED_TRANSPOSE:
121 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) {
122 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) {
123 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[i * N + k] * B[j * N + k];
124 }
125 }
126 break;
127 }
128 break;
129 case CEED_TRANSPOSE:
130 switch (transpose_B) {
131 case CEED_NOTRANSPOSE:
132 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) {
133 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) {
134 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[k * N + i] * B[k * N + j];
135 }
136 }
137 break;
138 case CEED_TRANSPOSE:
139 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) {
140 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) {
141 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[k * N + i] * B[j * N + k];
142 }
143 }
144 break;
145 }
146 break;
147 }
148 }
149
150 // @brief 3x3 Matrix-Matrix product, C = AB + C
MatMat3(const CeedScalar A[3][3],const CeedScalar B[3][3],const CeedTransposeMode transpose_A,const CeedTransposeMode transpose_B,CeedScalar C[3][3])151 CEED_QFUNCTION_HELPER void MatMat3(const CeedScalar A[3][3], const CeedScalar B[3][3], const CeedTransposeMode transpose_A,
152 const CeedTransposeMode transpose_B, CeedScalar C[3][3]) {
153 MatMatN((const CeedScalar *)A, (const CeedScalar *)B, 3, transpose_A, transpose_B, (CeedScalar *)C);
154 }
155
156 // @brief Unpack Kelvin-Mandel notation symmetric tensor into full tensor
KMUnpack(const CeedScalar v[6],CeedScalar A[3][3])157 CEED_QFUNCTION_HELPER void KMUnpack(const CeedScalar v[6], CeedScalar A[3][3]) {
158 const CeedScalar weight = 1 / sqrt(2.);
159 A[0][0] = v[0];
160 A[1][1] = v[1];
161 A[2][2] = v[2];
162 A[2][1] = A[1][2] = weight * v[3];
163 A[2][0] = A[0][2] = weight * v[4];
164 A[1][0] = A[0][1] = weight * v[5];
165 }
166
167 // @brief Pack full tensor into Kelvin-Mandel notation symmetric tensor
KMPack(const CeedScalar A[3][3],CeedScalar v[6])168 CEED_QFUNCTION_HELPER void KMPack(const CeedScalar A[3][3], CeedScalar v[6]) {
169 const CeedScalar weight = sqrt(2.);
170 v[0] = A[0][0];
171 v[1] = A[1][1];
172 v[2] = A[2][2];
173 v[3] = A[2][1] * weight;
174 v[4] = A[2][0] * weight;
175 v[5] = A[1][0] * weight;
176 }
177
178 // @brief Calculate metric tensor from mapping, g_{ij} = xi_{k,i} xi_{k,j} = dXdx^T dXdx
KMMetricTensor(const CeedScalar dXdx[3][3],CeedScalar km_g_ij[6])179 CEED_QFUNCTION_HELPER void KMMetricTensor(const CeedScalar dXdx[3][3], CeedScalar km_g_ij[6]) {
180 CeedScalar g_ij[3][3] = {{0.}};
181 MatMat3(dXdx, dXdx, CEED_TRANSPOSE, CEED_NOTRANSPOSE, g_ij);
182 KMPack(g_ij, km_g_ij);
183 }
184
185 // @brief Linear ramp evaluation
LinearRampCoefficient(CeedScalar amplitude,CeedScalar length,CeedScalar start,CeedScalar x)186 CEED_QFUNCTION_HELPER CeedScalar LinearRampCoefficient(CeedScalar amplitude, CeedScalar length, CeedScalar start, CeedScalar x) {
187 if (x < start) {
188 return amplitude;
189 } else if (x < start + length) {
190 return amplitude * ((x - start) * (-1 / length) + 1);
191 } else {
192 return 0;
193 }
194 }
195
196 /**
197 @brief Pack stored values at quadrature point
198
199 @param[in] Q Number of quadrature points
200 @param[in] i Current quadrature point
201 @param[in] start Starting index to store components
202 @param[in] num_comp Number of components to store
203 @param[in] values_at_qpnt Local values for quadrature point i
204 @param[out] stored Stored values
205
206 @return An error code: 0 - success, otherwise - failure
207 **/
StoredValuesPack(CeedInt Q,CeedInt i,CeedInt start,CeedInt num_comp,const CeedScalar * values_at_qpnt,CeedScalar * stored)208 CEED_QFUNCTION_HELPER int StoredValuesPack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *values_at_qpnt,
209 CeedScalar *stored) {
210 for (CeedInt j = 0; j < num_comp; j++) stored[(start + j) * Q + i] = values_at_qpnt[j];
211
212 return CEED_ERROR_SUCCESS;
213 }
214
215 /**
216 @brief Unpack stored values at quadrature point
217
218 @param[in] Q Number of quadrature points
219 @param[in] i Current quadrature point
220 @param[in] start Starting index to store components
221 @param[in] num_comp Number of components to store
222 @param[in] stored Stored values
223 @param[out] values_at_qpnt Local values for quadrature point i
224
225 @return An error code: 0 - success, otherwise - failure
226 **/
StoredValuesUnpack(CeedInt Q,CeedInt i,CeedInt start,CeedInt num_comp,const CeedScalar * stored,CeedScalar * values_at_qpnt)227 CEED_QFUNCTION_HELPER int StoredValuesUnpack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *stored,
228 CeedScalar *values_at_qpnt) {
229 for (CeedInt j = 0; j < num_comp; j++) values_at_qpnt[j] = stored[(start + j) * Q + i];
230
231 return CEED_ERROR_SUCCESS;
232 }
233
234 /**
235 @brief Unpack 3D element q_data at quadrature point
236
237 @param[in] Q Number of quadrature points
238 @param[in] i Current quadrature point
239 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`)
240 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian
241 @param[out] dXdx Inverse of the mapping Jacobian (shape [3][3])
242
243 @return An error code: 0 - success, otherwise - failure
244 **/
QdataUnpack_3D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar dXdx[3][3])245 CEED_QFUNCTION_HELPER int QdataUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[3][3]) {
246 StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
247 StoredValuesUnpack(Q, i, 1, 9, q_data, (CeedScalar *)dXdx);
248 return CEED_ERROR_SUCCESS;
249 }
250
251 /**
252 @brief Unpack boundary element q_data for 3D problem at quadrature point
253
254 @param[in] Q Number of quadrature points
255 @param[in] i Current quadrature point
256 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`)
257 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL`
258 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][3]), or `NULL`
259 @param[out] normal Components of the normal vector (shape [3]), or `NULL`
260
261 @return An error code: 0 - success, otherwise - failure
262 **/
QdataBoundaryUnpack_3D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar dXdx[2][3],CeedScalar normal[3])263 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][3],
264 CeedScalar normal[3]) {
265 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
266 if (normal) StoredValuesUnpack(Q, i, 1, 3, q_data, normal);
267 if (dXdx) StoredValuesUnpack(Q, i, 4, 6, q_data, (CeedScalar *)dXdx);
268 return CEED_ERROR_SUCCESS;
269 }
270
271 /**
272 @brief Unpack 2D element q_data at quadrature point
273
274 @param[in] Q Number of quadrature points
275 @param[in] i Current quadrature point
276 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`)
277 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian
278 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][2])
279
280 @return An error code: 0 - success, otherwise - failure
281 **/
QdataUnpack_2D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar dXdx[2][2])282 CEED_QFUNCTION_HELPER int QdataUnpack_2D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][2]) {
283 StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
284 StoredValuesUnpack(Q, i, 1, 4, q_data, (CeedScalar *)dXdx);
285 return CEED_ERROR_SUCCESS;
286 }
287
288 /**
289 @brief Unpack boundary element q_data for 2D problem at quadrature point
290
291 @param[in] Q Number of quadrature points
292 @param[in] i Current quadrature point
293 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary2d`)
294 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL`
295 @param[out] normal Components of the normal vector (shape [2]), or `NULL`
296
297 @return An error code: 0 - success, otherwise - failure
298 **/
QdataBoundaryUnpack_2D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar normal[2])299 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_2D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar normal[2]) {
300 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
301 if (normal) StoredValuesUnpack(Q, i, 1, 2, q_data, normal);
302 return CEED_ERROR_SUCCESS;
303 }
304