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