1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 #pragma once 8 9 #include <ceed.h> 10 #include <math.h> 11 12 #ifndef M_PI 13 #define M_PI 3.14159265358979323846 14 #endif 15 16 CEED_QFUNCTION_HELPER CeedScalar Max(CeedScalar a, CeedScalar b) { return a < b ? b : a; } 17 CEED_QFUNCTION_HELPER CeedScalar Min(CeedScalar a, CeedScalar b) { return a < b ? a : b; } 18 19 CEED_QFUNCTION_HELPER void SwapScalar(CeedScalar *a, CeedScalar *b) { 20 CeedScalar temp = *a; 21 *a = *b; 22 *b = temp; 23 } 24 25 CEED_QFUNCTION_HELPER CeedScalar Square(CeedScalar x) { return x * x; } 26 CEED_QFUNCTION_HELPER CeedScalar Cube(CeedScalar x) { return x * x * x; } 27 28 // @brief Scale vector of length N by scalar alpha 29 CEED_QFUNCTION_HELPER void ScaleN(CeedScalar *u, const CeedScalar alpha, const CeedInt N) { 30 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) u[i] *= alpha; 31 } 32 33 // @brief Set vector of length N to a value alpha 34 CEED_QFUNCTION_HELPER void SetValueN(CeedScalar *u, const CeedScalar alpha, const CeedInt N) { 35 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) u[i] = alpha; 36 } 37 38 // @brief Copy N elements from x to y 39 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]; } 40 41 // @brief Copy 3x3 matrix from A to B 42 CEED_QFUNCTION_HELPER void CopyMat3(const CeedScalar A[3][3], CeedScalar B[3][3]) { CopyN((const CeedScalar *)A, (CeedScalar *)B, 9); } 43 44 // @brief Dot product of vectors with N elements 45 CEED_QFUNCTION_HELPER CeedScalar DotN(const CeedScalar *u, const CeedScalar *v, const CeedInt N) { 46 CeedScalar output = 0; 47 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) output += u[i] * v[i]; 48 return output; 49 } 50 51 // @brief Dot product of 3 element vectors 52 CEED_QFUNCTION_HELPER CeedScalar Dot3(const CeedScalar *u, const CeedScalar *v) { return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]; } 53 54 // @brief \ell^2 norm of 3 element vectors 55 CEED_QFUNCTION_HELPER CeedScalar Norm3(const CeedScalar *u) { return sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]); } 56 57 // @brief Cross product of vectors with 3 elements 58 CEED_QFUNCTION_HELPER void Cross3(const CeedScalar u[3], const CeedScalar v[3], CeedScalar w[3]) { 59 w[0] = (u[1] * v[2]) - (u[2] * v[1]); 60 w[1] = (u[2] * v[0]) - (u[0] * v[2]); 61 w[2] = (u[0] * v[1]) - (u[1] * v[0]); 62 } 63 64 // @brief Curl of vector given its gradient 65 CEED_QFUNCTION_HELPER void Curl3(const CeedScalar gradient[3][3], CeedScalar v[3]) { 66 v[0] = gradient[2][1] - gradient[1][2]; 67 v[1] = gradient[0][2] - gradient[2][0]; 68 v[2] = gradient[1][0] - gradient[0][1]; 69 } 70 71 // @brief Matrix vector product, b = Ax + b. A is NxM, x is M, b is N 72 CEED_QFUNCTION_HELPER void MatVecNM(const CeedScalar *A, const CeedScalar *x, const CeedInt N, const CeedInt M, const CeedTransposeMode transpose_A, 73 CeedScalar *b) { 74 switch (transpose_A) { 75 case CEED_NOTRANSPOSE: 76 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) b[i] += DotN(&A[i * M], x, M); 77 break; 78 case CEED_TRANSPOSE: 79 CeedPragmaSIMD for (CeedInt i = 0; i < M; i++) { CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) b[i] += A[j * M + i] * x[j]; } 80 break; 81 } 82 } 83 84 // @brief 3x3 Matrix vector product b = Ax + b. 85 CEED_QFUNCTION_HELPER void MatVec3(const CeedScalar A[3][3], const CeedScalar x[3], const CeedTransposeMode transpose_A, CeedScalar b[3]) { 86 MatVecNM((const CeedScalar *)A, (const CeedScalar *)x, 3, 3, transpose_A, (CeedScalar *)b); 87 } 88 89 // @brief Matrix-Matrix product, B = DA + B, where D is diagonal. 90 // @details A is NxM, D is diagonal NxN, represented by a vector of length N, and B is NxM. Optionally, A may be transposed. 91 CEED_QFUNCTION_HELPER void MatDiagNM(const CeedScalar *A, const CeedScalar *D, const CeedInt N, const CeedInt M, const CeedTransposeMode transpose_A, 92 CeedScalar *B) { 93 switch (transpose_A) { 94 case CEED_NOTRANSPOSE: 95 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]; } 96 break; 97 case CEED_TRANSPOSE: 98 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]; } 99 break; 100 } 101 } 102 103 // @brief 3x3 Matrix-Matrix product, B = DA + B, where D is diagonal. 104 // @details Optionally, A may be transposed. 105 CEED_QFUNCTION_HELPER void MatDiag3(const CeedScalar A[3][3], const CeedScalar D[3], const CeedTransposeMode transpose_A, CeedScalar B[3][3]) { 106 MatDiagNM((const CeedScalar *)A, (const CeedScalar *)D, 3, 3, transpose_A, (CeedScalar *)B); 107 } 108 // @brief NxN Matrix-Matrix product, C = AB + C 109 CEED_QFUNCTION_HELPER void MatMatN(const CeedScalar *A, const CeedScalar *B, const CeedInt N, const CeedTransposeMode transpose_A, 110 const CeedTransposeMode transpose_B, CeedScalar *C) { 111 switch (transpose_A) { 112 case CEED_NOTRANSPOSE: 113 switch (transpose_B) { 114 case CEED_NOTRANSPOSE: 115 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) { 116 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { 117 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[i * N + k] * B[k * N + j]; 118 } 119 } 120 break; 121 case CEED_TRANSPOSE: 122 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) { 123 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { 124 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[i * N + k] * B[j * N + k]; 125 } 126 } 127 break; 128 } 129 break; 130 case CEED_TRANSPOSE: 131 switch (transpose_B) { 132 case CEED_NOTRANSPOSE: 133 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) { 134 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { 135 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[k * N + i] * B[k * N + j]; 136 } 137 } 138 break; 139 case CEED_TRANSPOSE: 140 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) { 141 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { 142 CeedPragmaSIMD for (CeedInt k = 0; k < N; k++) C[i * N + j] += A[k * N + i] * B[j * N + k]; 143 } 144 } 145 break; 146 } 147 break; 148 } 149 } 150 151 // @brief 3x3 Matrix-Matrix product, C = AB + C 152 CEED_QFUNCTION_HELPER void MatMat3(const CeedScalar A[3][3], const CeedScalar B[3][3], const CeedTransposeMode transpose_A, 153 const CeedTransposeMode transpose_B, CeedScalar C[3][3]) { 154 MatMatN((const CeedScalar *)A, (const CeedScalar *)B, 3, transpose_A, transpose_B, (CeedScalar *)C); 155 } 156 157 // @brief Unpack Kelvin-Mandel notation symmetric tensor into full tensor 158 CEED_QFUNCTION_HELPER void KMUnpack(const CeedScalar v[6], CeedScalar A[3][3]) { 159 const CeedScalar weight = 1 / sqrt(2.); 160 A[0][0] = v[0]; 161 A[1][1] = v[1]; 162 A[2][2] = v[2]; 163 A[2][1] = A[1][2] = weight * v[3]; 164 A[2][0] = A[0][2] = weight * v[4]; 165 A[1][0] = A[0][1] = weight * v[5]; 166 } 167 168 // @brief Pack full tensor into Kelvin-Mandel notation symmetric tensor 169 CEED_QFUNCTION_HELPER void KMPack(const CeedScalar A[3][3], CeedScalar v[6]) { 170 const CeedScalar weight = sqrt(2.); 171 v[0] = A[0][0]; 172 v[1] = A[1][1]; 173 v[2] = A[2][2]; 174 v[3] = A[2][1] * weight; 175 v[4] = A[2][0] * weight; 176 v[5] = A[1][0] * weight; 177 } 178 179 // @brief Calculate metric tensor from mapping, g_{ij} = xi_{k,i} xi_{k,j} = dXdx^T dXdx 180 CEED_QFUNCTION_HELPER void KMMetricTensor(const CeedScalar dXdx[3][3], CeedScalar km_g_ij[6]) { 181 CeedScalar g_ij[3][3] = {{0.}}; 182 MatMat3(dXdx, dXdx, CEED_TRANSPOSE, CEED_NOTRANSPOSE, g_ij); 183 KMPack(g_ij, km_g_ij); 184 } 185 186 // @brief Linear ramp evaluation 187 CEED_QFUNCTION_HELPER CeedScalar LinearRampCoefficient(CeedScalar amplitude, CeedScalar length, CeedScalar start, CeedScalar x) { 188 if (x < start) { 189 return amplitude; 190 } else if (x < start + length) { 191 return amplitude * ((x - start) * (-1 / length) + 1); 192 } else { 193 return 0; 194 } 195 } 196 197 /** 198 @brief Pack stored values at quadrature point 199 200 @param[in] Q Number of quadrature points 201 @param[in] i Current quadrature point 202 @param[in] start Starting index to store components 203 @param[in] num_comp Number of components to store 204 @param[in] values_at_qpnt Local values for quadrature point i 205 @param[out] stored Stored values 206 207 @return An error code: 0 - success, otherwise - failure 208 **/ 209 CEED_QFUNCTION_HELPER int StoredValuesPack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *values_at_qpnt, 210 CeedScalar *stored) { 211 for (CeedInt j = 0; j < num_comp; j++) stored[(start + j) * Q + i] = values_at_qpnt[j]; 212 213 return CEED_ERROR_SUCCESS; 214 } 215 216 /** 217 @brief Unpack 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] stored Stored values 224 @param[out] values_at_qpnt Local values for quadrature point i 225 226 @return An error code: 0 - success, otherwise - failure 227 **/ 228 CEED_QFUNCTION_HELPER int StoredValuesUnpack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *stored, 229 CeedScalar *values_at_qpnt) { 230 for (CeedInt j = 0; j < num_comp; j++) values_at_qpnt[j] = stored[(start + j) * Q + i]; 231 232 return CEED_ERROR_SUCCESS; 233 } 234 235 /** 236 @brief Unpack 3D element q_data at quadrature point 237 238 @param[in] Q Number of quadrature points 239 @param[in] i Current quadrature point 240 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`) 241 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian 242 @param[out] dXdx Inverse of the mapping Jacobian (shape [3][3]) 243 244 @return An error code: 0 - success, otherwise - failure 245 **/ 246 CEED_QFUNCTION_HELPER int QdataUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[3][3]) { 247 StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 248 StoredValuesUnpack(Q, i, 1, 9, q_data, (CeedScalar *)dXdx); 249 return CEED_ERROR_SUCCESS; 250 } 251 252 /** 253 @brief Unpack boundary element q_data for 3D problem at quadrature point 254 255 @param[in] Q Number of quadrature points 256 @param[in] i Current quadrature point 257 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`) 258 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL` 259 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][3]), or `NULL` 260 @param[out] normal Components of the normal vector (shape [3]), or `NULL` 261 262 @return An error code: 0 - success, otherwise - failure 263 **/ 264 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][3], 265 CeedScalar normal[3]) { 266 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 267 if (normal) StoredValuesUnpack(Q, i, 1, 3, q_data, normal); 268 if (dXdx) StoredValuesUnpack(Q, i, 4, 6, q_data, (CeedScalar *)dXdx); 269 return CEED_ERROR_SUCCESS; 270 } 271 272 /** 273 @brief Unpack boundary element q_data for 3D problem at quadrature point 274 275 @param[in] Q Number of quadrature points 276 @param[in] i Current quadrature point 277 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary`) 278 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL` 279 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][3]), or `NULL` 280 @param[out] normal Components of the normal vector (shape [3]), or `NULL` 281 282 @return An error code: 0 - success, otherwise - failure 283 **/ 284 CEED_QFUNCTION_HELPER int QdataBoundaryGradientUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][3], 285 CeedScalar normal[3]) { 286 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 287 if (dXdx) StoredValuesUnpack(Q, i, 1, 9, q_data, (CeedScalar *)dXdx); 288 if (normal) StoredValuesUnpack(Q, i, 10, 3, q_data, normal); 289 return CEED_ERROR_SUCCESS; 290 } 291 292 /** 293 @brief Unpack 2D element q_data at quadrature point 294 295 @param[in] Q Number of quadrature points 296 @param[in] i Current quadrature point 297 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`) 298 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian 299 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][2]) 300 301 @return An error code: 0 - success, otherwise - failure 302 **/ 303 CEED_QFUNCTION_HELPER int QdataUnpack_2D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][2]) { 304 StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 305 StoredValuesUnpack(Q, i, 1, 4, q_data, (CeedScalar *)dXdx); 306 return CEED_ERROR_SUCCESS; 307 } 308 309 /** 310 @brief Unpack boundary element q_data for 2D problem at quadrature point 311 312 @param[in] Q Number of quadrature points 313 @param[in] i Current quadrature point 314 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:SetupBoundary2d`) 315 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL` 316 @param[out] normal Components of the normal vector (shape [2]), or `NULL` 317 318 @return An error code: 0 - success, otherwise - failure 319 **/ 320 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_2D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar normal[2]) { 321 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 322 if (normal) StoredValuesUnpack(Q, i, 1, 2, q_data, normal); 323 return CEED_ERROR_SUCCESS; 324 } 325