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