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 Cross product of vectors with 3 elements 54 CEED_QFUNCTION_HELPER void Cross3(const CeedScalar u[3], const CeedScalar v[3], CeedScalar w[3]) { 55 w[0] = (u[1] * v[2]) - (u[2] * v[1]); 56 w[1] = (u[2] * v[0]) - (u[0] * v[2]); 57 w[2] = (u[0] * v[1]) - (u[1] * v[0]); 58 } 59 60 // @brief Curl of vector given its gradient 61 CEED_QFUNCTION_HELPER void Curl3(const CeedScalar gradient[3][3], CeedScalar v[3]) { 62 v[0] = gradient[2][1] - gradient[1][2]; 63 v[1] = gradient[0][2] - gradient[2][0]; 64 v[2] = gradient[1][0] - gradient[0][1]; 65 } 66 67 // @brief Matrix vector product, b = Ax + b. A is NxM, x is M, b is N 68 CEED_QFUNCTION_HELPER void MatVecNM(const CeedScalar *A, const CeedScalar *x, const CeedInt N, const CeedInt M, const CeedTransposeMode transpose_A, 69 CeedScalar *b) { 70 switch (transpose_A) { 71 case CEED_NOTRANSPOSE: 72 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) b[i] += DotN(&A[i * M], x, M); 73 break; 74 case CEED_TRANSPOSE: 75 CeedPragmaSIMD for (CeedInt i = 0; i < M; i++) { CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) b[i] += A[j * M + i] * x[j]; } 76 break; 77 } 78 } 79 80 // @brief 3x3 Matrix vector product b = Ax + b. 81 CEED_QFUNCTION_HELPER void MatVec3(const CeedScalar A[3][3], const CeedScalar x[3], const CeedTransposeMode transpose_A, CeedScalar b[3]) { 82 MatVecNM((const CeedScalar *)A, (const CeedScalar *)x, 3, 3, transpose_A, (CeedScalar *)b); 83 } 84 85 // @brief Matrix-Matrix product, B = DA + B, where D is diagonal. 86 // @details A is NxM, D is diagonal NxN, represented by a vector of length N, and B is NxM. Optionally, A may be transposed. 87 CEED_QFUNCTION_HELPER void MatDiagNM(const CeedScalar *A, const CeedScalar *D, const CeedInt N, const CeedInt M, const CeedTransposeMode transpose_A, 88 CeedScalar *B) { 89 switch (transpose_A) { 90 case CEED_NOTRANSPOSE: 91 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]; } 92 break; 93 case CEED_TRANSPOSE: 94 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]; } 95 break; 96 } 97 } 98 99 // @brief 3x3 Matrix-Matrix product, B = DA + B, where D is diagonal. 100 // @details Optionally, A may be transposed. 101 CEED_QFUNCTION_HELPER void MatDiag3(const CeedScalar A[3][3], const CeedScalar D[3], const CeedTransposeMode transpose_A, CeedScalar B[3][3]) { 102 MatDiagNM((const CeedScalar *)A, (const CeedScalar *)D, 3, 3, transpose_A, (CeedScalar *)B); 103 } 104 // @brief NxN Matrix-Matrix product, C = AB + C 105 CEED_QFUNCTION_HELPER void MatMatN(const CeedScalar *A, const CeedScalar *B, const CeedInt N, const CeedTransposeMode transpose_A, 106 const CeedTransposeMode transpose_B, CeedScalar *C) { 107 switch (transpose_A) { 108 case CEED_NOTRANSPOSE: 109 switch (transpose_B) { 110 case CEED_NOTRANSPOSE: 111 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) { 112 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { 113 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[i * N + k] * B[k * N + j]; 114 } 115 } 116 break; 117 case CEED_TRANSPOSE: 118 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) { 119 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { 120 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[i * N + k] * B[j * N + k]; 121 } 122 } 123 break; 124 } 125 break; 126 case CEED_TRANSPOSE: 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[k * N + i] * 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[k * N + i] * B[j * N + k]; 139 } 140 } 141 break; 142 } 143 break; 144 } 145 } 146 147 // @brief 3x3 Matrix-Matrix product, C = AB + C 148 CEED_QFUNCTION_HELPER void MatMat3(const CeedScalar A[3][3], const CeedScalar B[3][3], const CeedTransposeMode transpose_A, 149 const CeedTransposeMode transpose_B, CeedScalar C[3][3]) { 150 MatMatN((const CeedScalar *)A, (const CeedScalar *)B, 3, transpose_A, transpose_B, (CeedScalar *)C); 151 } 152 153 // @brief Unpack Kelvin-Mandel notation symmetric tensor into full tensor 154 CEED_QFUNCTION_HELPER void KMUnpack(const CeedScalar v[6], CeedScalar A[3][3]) { 155 const CeedScalar weight = 1 / sqrt(2.); 156 A[0][0] = v[0]; 157 A[1][1] = v[1]; 158 A[2][2] = v[2]; 159 A[2][1] = A[1][2] = weight * v[3]; 160 A[2][0] = A[0][2] = weight * v[4]; 161 A[1][0] = A[0][1] = weight * v[5]; 162 } 163 164 // @brief Pack full tensor into Kelvin-Mandel notation symmetric tensor 165 CEED_QFUNCTION_HELPER void KMPack(const CeedScalar A[3][3], CeedScalar v[6]) { 166 const CeedScalar weight = sqrt(2.); 167 v[0] = A[0][0]; 168 v[1] = A[1][1]; 169 v[2] = A[2][2]; 170 v[3] = A[2][1] * weight; 171 v[4] = A[2][0] * weight; 172 v[5] = A[1][0] * weight; 173 } 174 175 // @brief Calculate metric tensor from mapping, g_{ij} = xi_{k,i} xi_{k,j} = dXdx^T dXdx 176 CEED_QFUNCTION_HELPER void KMMetricTensor(const CeedScalar dXdx[3][3], CeedScalar km_g_ij[6]) { 177 CeedScalar g_ij[3][3] = {{0.}}; 178 MatMat3(dXdx, dXdx, CEED_TRANSPOSE, CEED_NOTRANSPOSE, g_ij); 179 KMPack(g_ij, km_g_ij); 180 } 181 182 // @brief Linear ramp evaluation 183 CEED_QFUNCTION_HELPER CeedScalar LinearRampCoefficient(CeedScalar amplitude, CeedScalar length, CeedScalar start, CeedScalar x) { 184 if (x < start) { 185 return amplitude; 186 } else if (x < start + length) { 187 return amplitude * ((x - start) * (-1 / length) + 1); 188 } else { 189 return 0; 190 } 191 } 192 193 /** 194 @brief Pack stored values at quadrature point 195 196 @param[in] Q Number of quadrature points 197 @param[in] i Current quadrature point 198 @param[in] start Starting index to store components 199 @param[in] num_comp Number of components to store 200 @param[in] values_at_qpnt Local values for quadrature point i 201 @param[out] stored Stored values 202 203 @return An error code: 0 - success, otherwise - failure 204 **/ 205 CEED_QFUNCTION_HELPER int StoredValuesPack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *values_at_qpnt, 206 CeedScalar *stored) { 207 for (CeedInt j = 0; j < num_comp; j++) stored[(start + j) * Q + i] = values_at_qpnt[j]; 208 209 return CEED_ERROR_SUCCESS; 210 } 211 212 /** 213 @brief Unpack stored values at quadrature point 214 215 @param[in] Q Number of quadrature points 216 @param[in] i Current quadrature point 217 @param[in] start Starting index to store components 218 @param[in] num_comp Number of components to store 219 @param[in] stored Stored values 220 @param[out] values_at_qpnt Local values for quadrature point i 221 222 @return An error code: 0 - success, otherwise - failure 223 **/ 224 CEED_QFUNCTION_HELPER int StoredValuesUnpack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *stored, 225 CeedScalar *values_at_qpnt) { 226 for (CeedInt j = 0; j < num_comp; j++) values_at_qpnt[j] = stored[(start + j) * Q + i]; 227 228 return CEED_ERROR_SUCCESS; 229 } 230 231 /** 232 @brief Unpack N-D element q_data at quadrature point 233 234 @param[in] dim Dimension of the element 235 @param[in] Q Number of quadrature points 236 @param[in] i Current quadrature point 237 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`) 238 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL` 239 @param[out] dXdx Inverse of the mapping Jacobian (shape [dim][dim]), or `NULL` 240 **/ 241 CEED_QFUNCTION_HELPER void QdataUnpack_ND(CeedInt dim, CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar *dXdx) { 242 switch (dim) { 243 case 2: 244 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 245 if (dXdx) StoredValuesUnpack(Q, i, 1, 4, q_data, dXdx); 246 break; 247 case 3: 248 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 249 if (dXdx) StoredValuesUnpack(Q, i, 1, 9, q_data, dXdx); 250 break; 251 } 252 } 253 254 /** 255 @brief Unpack boundary element q_data for N-D problem at quadrature point 256 257 @param[in] Q Number of quadrature points 258 @param[in] i Current quadrature point 259 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`) 260 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL` 261 @param[out] dXdx Inverse of the mapping Jacobian (shape [dim - 1][dim]), or `NULL` 262 @param[out] normal Components of the normal vector (shape [dim]), or `NULL` 263 **/ 264 CEED_QFUNCTION_HELPER void QdataBoundaryUnpack_ND(CeedInt dim, CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar *dXdx, 265 CeedScalar *normal) { 266 switch (dim) { 267 case 2: 268 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 269 if (normal) StoredValuesUnpack(Q, i, 1, 2, q_data, normal); 270 break; 271 case 3: 272 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 273 if (normal) StoredValuesUnpack(Q, i, 1, 3, q_data, normal); 274 if (dXdx) StoredValuesUnpack(Q, i, 4, 6, q_data, (CeedScalar *)dXdx); 275 break; 276 } 277 } 278 279 /** 280 @brief Unpack 3D element q_data at quadrature point 281 282 @param[in] Q Number of quadrature points 283 @param[in] i Current quadrature point 284 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`) 285 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian 286 @param[out] dXdx Inverse of the mapping Jacobian (shape [3][3]) 287 288 @return An error code: 0 - success, otherwise - failure 289 **/ 290 CEED_QFUNCTION_HELPER int QdataUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[3][3]) { 291 QdataUnpack_ND(3, Q, i, q_data, wdetJ, (CeedScalar *)dXdx); 292 return CEED_ERROR_SUCCESS; 293 } 294 295 /** 296 @brief Unpack boundary element q_data for 3D problem at quadrature point 297 298 @param[in] Q Number of quadrature points 299 @param[in] i Current quadrature point 300 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`) 301 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL` 302 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][3]), or `NULL` 303 @param[out] normal Components of the normal vector (shape [3]), or `NULL` 304 305 @return An error code: 0 - success, otherwise - failure 306 **/ 307 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][3], 308 CeedScalar normal[3]) { 309 QdataBoundaryUnpack_ND(3, Q, i, q_data, wdetJ, (CeedScalar *)dXdx, normal); 310 return CEED_ERROR_SUCCESS; 311 } 312 313 /** 314 @brief Unpack boundary element q_data for 3D problem at quadrature point 315 316 @param[in] Q Number of quadrature points 317 @param[in] i Current quadrature point 318 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`) 319 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL` 320 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][3]), or `NULL` 321 @param[out] normal Components of the normal vector (shape [3]), or `NULL` 322 323 @return An error code: 0 - success, otherwise - failure 324 **/ 325 CEED_QFUNCTION_HELPER int QdataBoundaryGradientUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][3], 326 CeedScalar normal[3]) { 327 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 328 if (dXdx) StoredValuesUnpack(Q, i, 1, 9, q_data, (CeedScalar *)dXdx); 329 if (normal) StoredValuesUnpack(Q, i, 10, 3, q_data, normal); 330 return CEED_ERROR_SUCCESS; 331 } 332 333 /** 334 @brief Unpack 2D element q_data at quadrature point 335 336 @param[in] Q Number of quadrature points 337 @param[in] i Current quadrature point 338 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`) 339 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian 340 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][2]) 341 342 @return An error code: 0 - success, otherwise - failure 343 **/ 344 CEED_QFUNCTION_HELPER int QdataUnpack_2D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][2]) { 345 QdataUnpack_ND(2, Q, i, q_data, wdetJ, (CeedScalar *)dXdx); 346 return CEED_ERROR_SUCCESS; 347 } 348 349 /** 350 @brief Unpack boundary element q_data for 2D problem at quadrature point 351 352 @param[in] Q Number of quadrature points 353 @param[in] i Current quadrature point 354 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary2d`) 355 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL` 356 @param[out] normal Components of the normal vector (shape [2]), or `NULL` 357 358 @return An error code: 0 - success, otherwise - failure 359 **/ 360 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_2D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar normal[2]) { 361 QdataBoundaryUnpack_ND(3, Q, i, q_data, wdetJ, NULL, normal); 362 return CEED_ERROR_SUCCESS; 363 } 364