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 14 CEED_QFUNCTION_HELPER CeedScalar Max(CeedScalar a, CeedScalar b) { return a < b ? b : a; } 15 CEED_QFUNCTION_HELPER CeedScalar Min(CeedScalar a, CeedScalar b) { return a < b ? a : b; } 16 17 CEED_QFUNCTION_HELPER void SwapScalar(CeedScalar *a, CeedScalar *b) { 18 CeedScalar temp = *a; 19 *a = *b; 20 *b = temp; 21 } 22 23 CEED_QFUNCTION_HELPER CeedScalar Square(CeedScalar x) { return x * x; } 24 CEED_QFUNCTION_HELPER CeedScalar Cube(CeedScalar x) { return x * x * x; } 25 26 // @brief Scale vector of length N by scalar alpha 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 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 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 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 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 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 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 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 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 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 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 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 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. 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. 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. 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. 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 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 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 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 */ 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 */ 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 **/ 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 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 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 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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 **/ 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