1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3 #pragma once
4
5 #include <ceed/types.h>
6 #ifndef CEED_RUNNING_JIT_PASS
7 #include <math.h>
8 #endif
9
10 #ifndef M_PI
11 #define M_PI 3.14159265358979323846
12 #endif
13
Max(CeedScalar a,CeedScalar b)14 CEED_QFUNCTION_HELPER CeedScalar Max(CeedScalar a, CeedScalar b) { return a < b ? b : a; }
Min(CeedScalar a,CeedScalar b)15 CEED_QFUNCTION_HELPER CeedScalar Min(CeedScalar a, CeedScalar b) { return a < b ? a : b; }
16
SwapScalar(CeedScalar * a,CeedScalar * b)17 CEED_QFUNCTION_HELPER void SwapScalar(CeedScalar *a, CeedScalar *b) {
18 CeedScalar temp = *a;
19 *a = *b;
20 *b = temp;
21 }
22
Square(CeedScalar x)23 CEED_QFUNCTION_HELPER CeedScalar Square(CeedScalar x) { return x * x; }
Cube(CeedScalar x)24 CEED_QFUNCTION_HELPER CeedScalar Cube(CeedScalar x) { return x * x * x; }
25
26 // @brief Scale vector of length N by scalar alpha
ScaleN(CeedScalar * u,const CeedScalar alpha,const CeedInt N)27 CEED_QFUNCTION_HELPER void ScaleN(CeedScalar *u, const CeedScalar alpha, const CeedInt N) {
28 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) u[i] *= alpha;
29 }
30
31 // @brief Set vector of length N to a value alpha
SetValueN(CeedScalar * u,const CeedScalar alpha,const CeedInt N)32 CEED_QFUNCTION_HELPER void SetValueN(CeedScalar *u, const CeedScalar alpha, const CeedInt N) {
33 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) u[i] = alpha;
34 }
35
36 // @brief Copy N elements from x to y
CopyN(const CeedScalar * x,CeedScalar * y,const CeedInt N)37 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]; }
38
39 // @brief Copy 3x3 matrix from A to B
CopyMat3(const CeedScalar A[3][3],CeedScalar B[3][3])40 CEED_QFUNCTION_HELPER void CopyMat3(const CeedScalar A[3][3], CeedScalar B[3][3]) { CopyN((const CeedScalar *)A, (CeedScalar *)B, 9); }
41
42 // @brief Dot product of vectors with N elements
DotN(const CeedScalar * u,const CeedScalar * v,const CeedInt N)43 CEED_QFUNCTION_HELPER CeedScalar DotN(const CeedScalar *u, const CeedScalar *v, const CeedInt N) {
44 CeedScalar output = 0;
45 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) output += u[i] * v[i];
46 return output;
47 }
48
49 // @brief y = \alpha x + y
AXPY(CeedScalar alpha,const CeedScalar * x,CeedScalar * y,CeedInt N)50 CEED_QFUNCTION_HELPER void AXPY(CeedScalar alpha, const CeedScalar *x, CeedScalar *y, CeedInt N) {
51 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) y[i] += alpha * x[i];
52 }
53
54 // @brief Dot product of 3 element vectors
Dot3(const CeedScalar * u,const CeedScalar * v)55 CEED_QFUNCTION_HELPER CeedScalar Dot3(const CeedScalar *u, const CeedScalar *v) { return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]; }
56
57 // @brief Dot product of 2 element vectors
Dot2(const CeedScalar * u,const CeedScalar * v)58 CEED_QFUNCTION_HELPER CeedScalar Dot2(const CeedScalar *u, const CeedScalar *v) { return u[0] * v[0] + u[1] * v[1]; }
59
60 // @brief \ell^2 norm of 3 element vectors
Norm3(const CeedScalar * u)61 CEED_QFUNCTION_HELPER CeedScalar Norm3(const CeedScalar *u) { return sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]); }
62
63 // @brief \ell^2 norm of 2 element vectors
Norm2(const CeedScalar * u)64 CEED_QFUNCTION_HELPER CeedScalar Norm2(const CeedScalar *u) { return sqrt(u[0] * u[0] + u[1] * u[1]); }
65
66 // @brief Cross product of vectors with 3 elements
Cross3(const CeedScalar u[3],const CeedScalar v[3],CeedScalar w[3])67 CEED_QFUNCTION_HELPER void Cross3(const CeedScalar u[3], const CeedScalar v[3], CeedScalar w[3]) {
68 w[0] = (u[1] * v[2]) - (u[2] * v[1]);
69 w[1] = (u[2] * v[0]) - (u[0] * v[2]);
70 w[2] = (u[0] * v[1]) - (u[1] * v[0]);
71 }
72
73 // @brief Curl of vector given its gradient
Curl3(const CeedScalar gradient[3][3],CeedScalar v[3])74 CEED_QFUNCTION_HELPER void Curl3(const CeedScalar gradient[3][3], CeedScalar v[3]) {
75 v[0] = gradient[2][1] - gradient[1][2];
76 v[1] = gradient[0][2] - gradient[2][0];
77 v[2] = gradient[1][0] - gradient[0][1];
78 }
79
80 // @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)81 CEED_QFUNCTION_HELPER void MatVecNM(const CeedScalar *A, const CeedScalar *x, const CeedInt N, const CeedInt M, const CeedTransposeMode transpose_A,
82 CeedScalar *b) {
83 switch (transpose_A) {
84 case CEED_NOTRANSPOSE:
85 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) b[i] += DotN(&A[i * M], x, M);
86 break;
87 case CEED_TRANSPOSE:
88 CeedPragmaSIMD for (CeedInt i = 0; i < M; i++) { CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) b[i] += A[j * M + i] * x[j]; }
89 break;
90 }
91 }
92
93 // @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])94 CEED_QFUNCTION_HELPER void MatVec3(const CeedScalar A[3][3], const CeedScalar x[3], const CeedTransposeMode transpose_A, CeedScalar b[3]) {
95 MatVecNM((const CeedScalar *)A, (const CeedScalar *)x, 3, 3, transpose_A, (CeedScalar *)b);
96 }
97
98 // @brief 2x2 Matrix vector product b = Ax + b.
MatVec2(const CeedScalar A[2][2],const CeedScalar x[2],const CeedTransposeMode transpose_A,CeedScalar b[2])99 CEED_QFUNCTION_HELPER void MatVec2(const CeedScalar A[2][2], const CeedScalar x[2], const CeedTransposeMode transpose_A, CeedScalar b[2]) {
100 MatVecNM((const CeedScalar *)A, (const CeedScalar *)x, 2, 2, transpose_A, (CeedScalar *)b);
101 }
102
103 // @brief Matrix-Matrix product, B = DA + B, where D is diagonal.
104 // @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)105 CEED_QFUNCTION_HELPER void MatDiagNM(const CeedScalar *A, const CeedScalar *D, const CeedInt N, const CeedInt M, const CeedTransposeMode transpose_A,
106 CeedScalar *B) {
107 switch (transpose_A) {
108 case CEED_NOTRANSPOSE:
109 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]; }
110 break;
111 case CEED_TRANSPOSE:
112 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]; }
113 break;
114 }
115 }
116
117 // @brief 3x3 Matrix-Matrix product, B = DA + B, where D is diagonal.
118 // @details Optionally, A may be transposed.
MatDiag3(const CeedScalar A[3][3],const CeedScalar D[3],const CeedTransposeMode transpose_A,CeedScalar B[3][3])119 CEED_QFUNCTION_HELPER void MatDiag3(const CeedScalar A[3][3], const CeedScalar D[3], const CeedTransposeMode transpose_A, CeedScalar B[3][3]) {
120 MatDiagNM((const CeedScalar *)A, (const CeedScalar *)D, 3, 3, transpose_A, (CeedScalar *)B);
121 }
122 // @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)123 CEED_QFUNCTION_HELPER void MatMatN(const CeedScalar *A, const CeedScalar *B, const CeedInt N, const CeedTransposeMode transpose_A,
124 const CeedTransposeMode transpose_B, CeedScalar *C) {
125 switch (transpose_A) {
126 case CEED_NOTRANSPOSE:
127 switch (transpose_B) {
128 case CEED_NOTRANSPOSE:
129 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) {
130 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) {
131 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[i * N + k] * B[k * N + j];
132 }
133 }
134 break;
135 case CEED_TRANSPOSE:
136 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) {
137 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) {
138 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[i * N + k] * B[j * N + k];
139 }
140 }
141 break;
142 }
143 break;
144 case CEED_TRANSPOSE:
145 switch (transpose_B) {
146 case CEED_NOTRANSPOSE:
147 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) {
148 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) {
149 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[k * N + i] * B[k * N + j];
150 }
151 }
152 break;
153 case CEED_TRANSPOSE:
154 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) {
155 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) {
156 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[k * N + i] * B[j * N + k];
157 }
158 }
159 break;
160 }
161 break;
162 }
163 }
164
165 // @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])166 CEED_QFUNCTION_HELPER void MatMat3(const CeedScalar A[3][3], const CeedScalar B[3][3], const CeedTransposeMode transpose_A,
167 const CeedTransposeMode transpose_B, CeedScalar C[3][3]) {
168 MatMatN((const CeedScalar *)A, (const CeedScalar *)B, 3, transpose_A, transpose_B, (CeedScalar *)C);
169 }
170
171 // @brief 2x2 Matrix-Matrix product, C = AB + C
MatMat2(const CeedScalar A[2][2],const CeedScalar B[2][2],const CeedTransposeMode transpose_A,const CeedTransposeMode transpose_B,CeedScalar C[2][2])172 CEED_QFUNCTION_HELPER void MatMat2(const CeedScalar A[2][2], const CeedScalar B[2][2], const CeedTransposeMode transpose_A,
173 const CeedTransposeMode transpose_B, CeedScalar C[2][2]) {
174 MatMatN((const CeedScalar *)A, (const CeedScalar *)B, 2, transpose_A, transpose_B, (CeedScalar *)C);
175 }
176
177 /**
178 * @brief Calculate inverse of 2x2 matrix
179 *
180 * @param[in] A Input matrix
181 * @param[out] detJ_ptr Determinate of A, may be NULL is not desired
182 * @param[out] A_inv Output matrix inverse
183 */
MatInv2(const CeedScalar A[2][2],CeedScalar A_inv[2][2],CeedScalar * detJ_ptr)184 CEED_QFUNCTION_HELPER void MatInv2(const CeedScalar A[2][2], CeedScalar A_inv[2][2], CeedScalar *detJ_ptr) {
185 const CeedScalar detJ = A[0][0] * A[1][1] - A[1][0] * A[0][1];
186
187 A_inv[0][0] = A[1][1] / detJ;
188 A_inv[0][1] = -A[0][1] / detJ;
189 A_inv[1][0] = -A[1][0] / detJ;
190 A_inv[1][1] = A[0][0] / detJ;
191 if (detJ_ptr) *detJ_ptr = detJ;
192 }
193
194 /**
195 * @brief Calculate inverse of 3x3 matrix
196 *
197 * @param[in] A Input matrix
198 * @param[out] detJ_ptr Determinate of A, may be NULL is not desired
199 * @param[out] A_inv Output matrix inverse
200 */
MatInv3(const CeedScalar A[3][3],CeedScalar A_inv[3][3],CeedScalar * detJ_ptr)201 CEED_QFUNCTION_HELPER void MatInv3(const CeedScalar A[3][3], CeedScalar A_inv[3][3], CeedScalar *detJ_ptr) {
202 // Compute Adjugate of dxdX
203 A_inv[0][0] = A[1][1] * A[2][2] - A[1][2] * A[2][1];
204 A_inv[0][1] = A[0][2] * A[2][1] - A[0][1] * A[2][2];
205 A_inv[0][2] = A[0][1] * A[1][2] - A[0][2] * A[1][1];
206 A_inv[1][0] = A[1][2] * A[2][0] - A[1][0] * A[2][2];
207 A_inv[1][1] = A[0][0] * A[2][2] - A[0][2] * A[2][0];
208 A_inv[1][2] = A[0][2] * A[1][0] - A[0][0] * A[1][2];
209 A_inv[2][0] = A[1][0] * A[2][1] - A[1][1] * A[2][0];
210 A_inv[2][1] = A[0][1] * A[2][0] - A[0][0] * A[2][1];
211 A_inv[2][2] = A[0][0] * A[1][1] - A[0][1] * A[1][0];
212
213 const CeedScalar detJ = A[0][0] * A_inv[0][0] + A[1][0] * A_inv[0][1] + A[2][0] * A_inv[0][2];
214 ScaleN((CeedScalar *)A_inv, 1 / detJ, 9);
215 if (detJ_ptr) *detJ_ptr = detJ;
216 }
217
218 /**
219 @brief MxN Matrix-Matrix product, C = AB + C
220
221 C is NxM, A is NxP, B is PxM
222
223 @param[in] mat_A Row-major matrix `A`
224 @param[in] mat_B Row-major matrix `B`
225 @param[out] mat_C Row-major output matrix `C`
226 @param[in] N Number of rows of `C`
227 @param[in] M Number of columns of `C`
228 @param[in] P Number of columns of `A`/rows of `B`
229 **/
MatMatNM(const CeedScalar * mat_A,const CeedScalar * mat_B,CeedScalar * mat_C,CeedInt N,CeedInt M,CeedInt P)230 CEED_QFUNCTION_HELPER void MatMatNM(const CeedScalar *mat_A, const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt N, CeedInt M, CeedInt P) {
231 for (CeedInt i = 0; i < N; i++) {
232 for (CeedInt j = 0; j < M; j++) {
233 for (CeedInt k = 0; k < P; k++) mat_C[i * M + j] += mat_A[i * P + k] * mat_B[k * M + j];
234 }
235 }
236 }
237
238 // @brief Unpack Kelvin-Mandel notation symmetric tensor into full tensor
KMUnpack(const CeedScalar v[6],CeedScalar A[3][3])239 CEED_QFUNCTION_HELPER void KMUnpack(const CeedScalar v[6], CeedScalar A[3][3]) {
240 const CeedScalar weight = 1 / sqrt(2.);
241 A[0][0] = v[0];
242 A[1][1] = v[1];
243 A[2][2] = v[2];
244 A[2][1] = A[1][2] = weight * v[3];
245 A[2][0] = A[0][2] = weight * v[4];
246 A[1][0] = A[0][1] = weight * v[5];
247 }
248
249 // @brief Pack full tensor into Kelvin-Mandel notation symmetric tensor
KMPack(const CeedScalar A[3][3],CeedScalar v[6])250 CEED_QFUNCTION_HELPER void KMPack(const CeedScalar A[3][3], CeedScalar v[6]) {
251 const CeedScalar weight = sqrt(2.);
252 v[0] = A[0][0];
253 v[1] = A[1][1];
254 v[2] = A[2][2];
255 v[3] = A[2][1] * weight;
256 v[4] = A[2][0] * weight;
257 v[5] = A[1][0] * weight;
258 }
259
260 // @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])261 CEED_QFUNCTION_HELPER void KMMetricTensor(const CeedScalar dXdx[3][3], CeedScalar km_g_ij[6]) {
262 CeedScalar g_ij[3][3] = {{0.}};
263 MatMat3(dXdx, dXdx, CEED_TRANSPOSE, CEED_NOTRANSPOSE, g_ij);
264 KMPack(g_ij, km_g_ij);
265 }
266
267 /**
268 @brief Linear ramp evaluation from set amplitude to zero
269
270 ▲
271 │
272 a│-------┬.
273 │ ┊ `-.
274 │ ┊ `-.
275 │ ┊ `-.______
276 └───────┴─────────┴────────> x
277 s s+l
278
279 where "a" is `amplitude`, "s" is `start`, and "l" is `length`.
280
281 @param[in] amplitude Maximum value of the ramp
282 @param[in] length Length of the ramp
283 @param[in] start Location where ramp begins to reduce from `amplitude` to 0
284 @param[in] x Input location
285 @return Value of linear ramp function
286 **/
LinearRampCoefficient(CeedScalar amplitude,CeedScalar length,CeedScalar start,CeedScalar x)287 CEED_QFUNCTION_HELPER CeedScalar LinearRampCoefficient(CeedScalar amplitude, CeedScalar length, CeedScalar start, CeedScalar x) {
288 if (x < start) {
289 return amplitude;
290 } else if (x < start + length) {
291 return amplitude * ((x - start) * (-1 / length) + 1);
292 } else {
293 return 0;
294 }
295 }
296
297 /**
298 @brief Pack stored values at quadrature point
299
300 @param[in] Q Number of quadrature points
301 @param[in] i Current quadrature point
302 @param[in] start Starting index to store components
303 @param[in] num_comp Number of components to store
304 @param[in] values_at_qpnt Local values for quadrature point i
305 @param[out] stored Stored values
306
307 @return An error code: 0 - success, otherwise - failure
308 **/
StoredValuesPack(CeedInt Q,CeedInt i,CeedInt start,CeedInt num_comp,const CeedScalar * values_at_qpnt,CeedScalar * stored)309 CEED_QFUNCTION_HELPER int StoredValuesPack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *values_at_qpnt,
310 CeedScalar *stored) {
311 for (CeedInt j = 0; j < num_comp; j++) stored[(start + j) * Q + i] = values_at_qpnt[j];
312
313 return CEED_ERROR_SUCCESS;
314 }
315
316 /**
317 @brief Unpack stored values at quadrature point
318
319 @param[in] Q Number of quadrature points
320 @param[in] i Current quadrature point
321 @param[in] start Starting index to store components
322 @param[in] num_comp Number of components to store
323 @param[in] stored Stored values
324 @param[out] values_at_qpnt Local values for quadrature point i
325
326 @return An error code: 0 - success, otherwise - failure
327 **/
StoredValuesUnpack(CeedInt Q,CeedInt i,CeedInt start,CeedInt num_comp,const CeedScalar * stored,CeedScalar * values_at_qpnt)328 CEED_QFUNCTION_HELPER int StoredValuesUnpack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *stored,
329 CeedScalar *values_at_qpnt) {
330 for (CeedInt j = 0; j < num_comp; j++) values_at_qpnt[j] = stored[(start + j) * Q + i];
331
332 return CEED_ERROR_SUCCESS;
333 }
334
335 /**
336 @brief Unpack N-D element q_data at quadrature point
337
338 @param[in] dim Dimension of the element
339 @param[in] Q Number of quadrature points
340 @param[in] i Current quadrature point
341 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`)
342 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL`
343 @param[out] dXdx Inverse of the mapping Jacobian (shape [dim][dim]), or `NULL`
344
345 @return An error code: 0 - success, otherwise - failure
346 **/
QdataUnpack_ND(CeedInt dim,CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar * dXdx)347 CEED_QFUNCTION_HELPER int QdataUnpack_ND(CeedInt dim, CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar *dXdx) {
348 switch (dim) {
349 case 2:
350 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
351 if (dXdx) StoredValuesUnpack(Q, i, 1, 4, q_data, dXdx);
352 break;
353 case 3:
354 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
355 if (dXdx) StoredValuesUnpack(Q, i, 1, 9, q_data, dXdx);
356 break;
357 }
358 return CEED_ERROR_SUCCESS;
359 }
360
361 /**
362 @brief Unpack boundary element q_data for N-D problem at quadrature point
363
364 @param[in] dim Dimension of the element
365 @param[in] Q Number of quadrature points
366 @param[in] i Current quadrature point
367 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`)
368 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL`
369 @param[out] dXdx Inverse of the mapping Jacobian (shape [dim - 1][dim]), or `NULL`
370 @param[out] normal Components of the normal vector (shape [dim]), or `NULL`
371
372 @return An error code: 0 - success, otherwise - failure
373 **/
QdataBoundaryUnpack_ND(CeedInt dim,CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar * dXdx,CeedScalar * normal)374 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_ND(CeedInt dim, CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar *dXdx,
375 CeedScalar *normal) {
376 switch (dim) {
377 case 2:
378 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
379 if (normal) StoredValuesUnpack(Q, i, 1, 2, q_data, normal);
380 break;
381 case 3:
382 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
383 if (normal) StoredValuesUnpack(Q, i, 1, 3, q_data, normal);
384 if (dXdx) StoredValuesUnpack(Q, i, 4, 6, q_data, (CeedScalar *)dXdx);
385 break;
386 }
387 return CEED_ERROR_SUCCESS;
388 }
389
390 /**
391 @brief Unpack boundary element q_data for N-D problem at quadrature point
392
393 @param[in] dim Dimension of the element
394 @param[in] Q Number of quadrature points
395 @param[in] i Current quadrature point
396 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundaryGradient`)
397 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL`
398 @param[out] dXdx Inverse of the mapping Jacobian (shape [dim][dim]), or `NULL`
399 @param[out] normal Components of the normal vector (shape [dim]), or `NULL`
400
401 @return An error code: 0 - success, otherwise - failure
402 **/
QdataBoundaryGradientUnpack_ND(CeedInt dim,CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar * dXdx,CeedScalar * normal)403 CEED_QFUNCTION_HELPER int QdataBoundaryGradientUnpack_ND(CeedInt dim, CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ,
404 CeedScalar *dXdx, CeedScalar *normal) {
405 switch (dim) {
406 case 2:
407 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
408 if (dXdx) StoredValuesUnpack(Q, i, 1, 4, q_data, dXdx);
409 if (normal) StoredValuesUnpack(Q, i, 5, 2, q_data, normal);
410 break;
411 case 3:
412 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ);
413 if (dXdx) StoredValuesUnpack(Q, i, 1, 9, q_data, dXdx);
414 if (normal) StoredValuesUnpack(Q, i, 10, 3, q_data, normal);
415 break;
416 }
417 return CEED_ERROR_SUCCESS;
418 }
419
420 /**
421 @brief Unpack 3D element q_data at quadrature point
422
423 @param[in] Q Number of quadrature points
424 @param[in] i Current quadrature point
425 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`)
426 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian
427 @param[out] dXdx Inverse of the mapping Jacobian (shape [3][3])
428
429 @return An error code: 0 - success, otherwise - failure
430 **/
QdataUnpack_3D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar dXdx[3][3])431 CEED_QFUNCTION_HELPER int QdataUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[3][3]) {
432 return QdataUnpack_ND(3, Q, i, q_data, wdetJ, (CeedScalar *)dXdx);
433 }
434
435 /**
436 @brief Unpack boundary element q_data for 3D problem at quadrature point
437
438 @param[in] Q Number of quadrature points
439 @param[in] i Current quadrature point
440 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`)
441 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL`
442 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][3]), or `NULL`
443 @param[out] normal Components of the normal vector (shape [3]), or `NULL`
444
445 @return An error code: 0 - success, otherwise - failure
446 **/
QdataBoundaryUnpack_3D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar dXdx[2][3],CeedScalar normal[3])447 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][3],
448 CeedScalar normal[3]) {
449 return QdataBoundaryUnpack_ND(3, Q, i, q_data, wdetJ, (CeedScalar *)dXdx, normal);
450 }
451
452 /**
453 @brief Unpack boundary element q_data for 3D problem at quadrature point
454
455 @param[in] Q Number of quadrature points
456 @param[in] i Current quadrature point
457 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`)
458 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL`
459 @param[out] dXdx Inverse of the mapping Jacobian (shape [3][3]), or `NULL`
460 @param[out] normal Components of the normal vector (shape [3]), or `NULL`
461
462 @return An error code: 0 - success, otherwise - failure
463 **/
QdataBoundaryGradientUnpack_3D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar dXdx[3][3],CeedScalar normal[3])464 CEED_QFUNCTION_HELPER int QdataBoundaryGradientUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[3][3],
465 CeedScalar normal[3]) {
466 return QdataBoundaryGradientUnpack_ND(3, Q, i, q_data, wdetJ, (CeedScalar *)dXdx, normal);
467 }
468
469 /**
470 @brief Unpack 2D element q_data at quadrature point
471
472 @param[in] Q Number of quadrature points
473 @param[in] i Current quadrature point
474 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`)
475 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian
476 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][2])
477
478 @return An error code: 0 - success, otherwise - failure
479 **/
QdataUnpack_2D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar dXdx[2][2])480 CEED_QFUNCTION_HELPER int QdataUnpack_2D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][2]) {
481 QdataUnpack_ND(2, Q, i, q_data, wdetJ, (CeedScalar *)dXdx);
482 return CEED_ERROR_SUCCESS;
483 }
484
485 /**
486 @brief Unpack boundary element q_data for 2D problem at quadrature point
487
488 @param[in] Q Number of quadrature points
489 @param[in] i Current quadrature point
490 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary2d`)
491 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL`
492 @param[out] normal Components of the normal vector (shape [2]), or `NULL`
493
494 @return An error code: 0 - success, otherwise - failure
495 **/
QdataBoundaryUnpack_2D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar normal[2])496 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_2D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar normal[2]) {
497 QdataBoundaryUnpack_ND(3, Q, i, q_data, wdetJ, NULL, normal);
498 return CEED_ERROR_SUCCESS;
499 }
500
501 /**
502 @brief Unpack boundary element q_data for 2D problem at quadrature point
503
504 @param[in] Q Number of quadrature points
505 @param[in] i Current quadrature point
506 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`)
507 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL`
508 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][2]), or `NULL`
509 @param[out] normal Components of the normal vector (shape [2]), or `NULL`
510
511 @return An error code: 0 - success, otherwise - failure
512 **/
QdataBoundaryGradientUnpack_2D(CeedInt Q,CeedInt i,const CeedScalar * q_data,CeedScalar * wdetJ,CeedScalar dXdx[2][2],CeedScalar normal[2])513 CEED_QFUNCTION_HELPER int QdataBoundaryGradientUnpack_2D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][2],
514 CeedScalar normal[2]) {
515 return QdataBoundaryGradientUnpack_ND(2, Q, i, q_data, wdetJ, (CeedScalar *)dXdx, normal);
516 }
517
518 /**
519 @brief Unpack `CEED_EVAL_GRAD` QF input into quadrature-point local array
520
521 @param[in] Q Number of quadrature points
522 @param[in] i Current quadrature point
523 @param[in] num_comp Number of components of the input
524 @param[in] dim Topological dimension of the element (ie. number of derivative terms per component)
525 @param[in] grad QF gradient input, shape `[dim][num_comp][Q]`
526 @param[out] grad_local Gradient array at quadrature point Q, shape `[num_comp][dim]`
527 **/
GradUnpackND(CeedInt Q,CeedInt i,CeedInt num_comp,CeedInt dim,const CeedScalar * grad,CeedScalar * grad_local)528 CEED_QFUNCTION_HELPER void GradUnpackND(CeedInt Q, CeedInt i, CeedInt num_comp, CeedInt dim, const CeedScalar *grad, CeedScalar *grad_local) {
529 for (CeedInt d = 0; d < dim; d++) {
530 for (CeedInt c = 0; c < num_comp; c++) {
531 grad_local[dim * c + d] = grad[(Q * num_comp) * d + Q * c + i];
532 }
533 }
534 }
535
536 /**
537 @brief Unpack `CEED_EVAL_GRAD` QF input into quadrature-point local array for 3D elements
538
539 @param[in] Q Number of quadrature points
540 @param[in] i Current quadrature point
541 @param[in] num_comp Number of components of the input
542 @param[in] grad QF gradient input, shape `[3][num_comp][Q]`
543 @param[out] grad_local Gradient array at quadrature point Q, shape `[num_comp][3]`
544 **/
GradUnpack3D(CeedInt Q,CeedInt i,CeedInt num_comp,const CeedScalar * grad,CeedScalar (* grad_local)[3])545 CEED_QFUNCTION_HELPER void GradUnpack3D(CeedInt Q, CeedInt i, CeedInt num_comp, const CeedScalar *grad, CeedScalar (*grad_local)[3]) {
546 GradUnpackND(Q, i, num_comp, 3, grad, (CeedScalar *)grad_local);
547 }
548
549 /**
550 @brief Unpack `CEED_EVAL_GRAD` QF input into quadrature-point local array for 2D elements
551
552 @param[in] Q Number of quadrature points
553 @param[in] i Current quadrature point
554 @param[in] num_comp Number of components of the input
555 @param[in] grad QF gradient input, shape `[2][num_comp][Q]`
556 @param[out] grad_local Gradient array at quadrature point Q, shape `[num_comp][2]`
557 **/
GradUnpack2D(CeedInt Q,CeedInt i,CeedInt num_comp,const CeedScalar * grad,CeedScalar (* grad_local)[2])558 CEED_QFUNCTION_HELPER void GradUnpack2D(CeedInt Q, CeedInt i, CeedInt num_comp, const CeedScalar *grad, CeedScalar (*grad_local)[2]) {
559 GradUnpackND(Q, i, num_comp, 2, grad, (CeedScalar *)grad_local);
560 }
561
562 /**
563 @brief Calculate divergence from reference gradient
564
565 Given gradient array G_{ij} and inverse element mapping X_{ij}, then the divergence is
566
567 G_{ij} X{ji}
568
569 @param[in] grad_qn Gradient array, orientation [vector component][gradient direction]
570 @param[in] dXdx Inverse of the mapping Jacobian (shape [dim][dim])
571 @param[in] dim Dimension of the problem
572 @param[out] divergence The divergence
573 **/
DivergenceND(const CeedScalar * grad_qn,const CeedScalar * dXdx,const CeedInt dim,CeedScalar * divergence)574 CEED_QFUNCTION_HELPER void DivergenceND(const CeedScalar *grad_qn, const CeedScalar *dXdx, const CeedInt dim, CeedScalar *divergence) {
575 for (CeedInt i = 0; i < dim; i++) {
576 for (CeedInt j = 0; j < dim; j++) {
577 *divergence += grad_qn[i * dim + j] * dXdx[j * dim + i];
578 }
579 }
580 }
581
582 /**
583 @brief Calculate divergence from reference gradient for 3D problem
584
585 Given gradient array G_{ij} and inverse element mapping X_{ij}, then the divergence is
586
587 G_{ij} X{ji}
588
589 @param[in] grad_qn Gradient array, orientation [vector component][gradient direction]
590 @param[in] dXdx Inverse of the mapping Jacobian (shape [3][3])
591 @param[out] divergence The divergence
592 **/
Divergence3D(const CeedScalar grad_qn[3][3],const CeedScalar dXdx[3][3],CeedScalar * divergence)593 CEED_QFUNCTION_HELPER void Divergence3D(const CeedScalar grad_qn[3][3], const CeedScalar dXdx[3][3], CeedScalar *divergence) {
594 DivergenceND((const CeedScalar *)grad_qn, (const CeedScalar *)dXdx, 3, divergence);
595 }
596