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