1 // Copyright (c) 2017-2022, 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 8 #ifndef utils_h 9 #define utils_h 10 11 #include <ceed.h> 12 #include <math.h> 13 14 #ifndef M_PI 15 #define M_PI 3.14159265358979323846 16 #endif 17 18 CEED_QFUNCTION_HELPER CeedScalar Max(CeedScalar a, CeedScalar b) { return a < b ? b : a; } 19 CEED_QFUNCTION_HELPER CeedScalar Min(CeedScalar a, CeedScalar b) { return a < b ? a : b; } 20 21 CEED_QFUNCTION_HELPER void SwapScalar(CeedScalar *a, CeedScalar *b) { 22 CeedScalar temp = *a; 23 *a = *b; 24 *b = temp; 25 } 26 27 CEED_QFUNCTION_HELPER CeedScalar Square(CeedScalar x) { return x * x; } 28 CEED_QFUNCTION_HELPER CeedScalar Cube(CeedScalar x) { return x * x * x; } 29 30 // @brief Scale vector of length N by scalar alpha 31 CEED_QFUNCTION_HELPER void ScaleN(CeedScalar *u, const CeedScalar alpha, const CeedInt N) { 32 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) u[i] *= alpha; 33 } 34 35 // @brief Set vector of length N to a value alpha 36 CEED_QFUNCTION_HELPER void SetValueN(CeedScalar *u, const CeedScalar alpha, const CeedInt N) { 37 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) u[i] = alpha; 38 } 39 40 // @brief Copy N elements from x to y 41 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]; } 42 43 // @brief Copy 3x3 matrix from A to B 44 CEED_QFUNCTION_HELPER void CopyMat3(const CeedScalar A[3][3], CeedScalar B[3][3]) { CopyN((const CeedScalar *)A, (CeedScalar *)B, 9); } 45 46 // @brief Dot product of vectors with N elements 47 CEED_QFUNCTION_HELPER CeedScalar DotN(const CeedScalar *u, const CeedScalar *v, const CeedInt N) { 48 CeedScalar output = 0; 49 CeedPragmaSIMD for (CeedInt i = 0; i < N; i++) output += u[i] * v[i]; 50 return output; 51 } 52 53 // @brief Dot product of 3 element vectors 54 CEED_QFUNCTION_HELPER CeedScalar Dot3(const CeedScalar *u, const CeedScalar *v) { return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]; } 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 108 // @brief 3x3 Matrix-Matrix product, C = AB + C 109 CEED_QFUNCTION_HELPER void MatMat3(const CeedScalar A[3][3], const CeedScalar B[3][3], const CeedTransposeMode transpose_A, 110 const CeedTransposeMode transpose_B, CeedScalar C[3][3]) { 111 switch (transpose_A) { 112 case CEED_NOTRANSPOSE: 113 switch (transpose_B) { 114 case CEED_NOTRANSPOSE: 115 CeedPragmaSIMD for (CeedInt i = 0; i < 3; i++) { 116 CeedPragmaSIMD for (CeedInt j = 0; j < 3; j++) { CeedPragmaSIMD for (CeedInt k = 0; k < 3; k++) C[i][j] += A[i][k] * B[k][j]; } 117 } 118 break; 119 case CEED_TRANSPOSE: 120 CeedPragmaSIMD for (CeedInt i = 0; i < 3; i++) { 121 CeedPragmaSIMD for (CeedInt j = 0; j < 3; j++) { CeedPragmaSIMD for (CeedInt k = 0; k < 3; k++) C[i][j] += A[i][k] * B[j][k]; } 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 < 3; i++) { 130 CeedPragmaSIMD for (CeedInt j = 0; j < 3; j++) { CeedPragmaSIMD for (CeedInt k = 0; k < 3; k++) C[i][j] += A[k][i] * B[k][j]; } 131 } 132 break; 133 case CEED_TRANSPOSE: 134 CeedPragmaSIMD for (CeedInt i = 0; i < 3; i++) { 135 CeedPragmaSIMD for (CeedInt j = 0; j < 3; j++) { CeedPragmaSIMD for (CeedInt k = 0; k < 3; k++) C[i][j] += A[k][i] * B[j][k]; } 136 } 137 break; 138 } 139 break; 140 } 141 } 142 143 // @brief Unpack Kelvin-Mandel notation symmetric tensor into full tensor 144 CEED_QFUNCTION_HELPER void KMUnpack(const CeedScalar v[6], CeedScalar A[3][3]) { 145 const CeedScalar weight = 1 / sqrt(2.); 146 A[0][0] = v[0]; 147 A[1][1] = v[1]; 148 A[2][2] = v[2]; 149 A[2][1] = A[1][2] = weight * v[3]; 150 A[2][0] = A[0][2] = weight * v[4]; 151 A[1][0] = A[0][1] = weight * v[5]; 152 } 153 154 // @brief Pack full tensor into Kelvin-Mandel notation symmetric tensor 155 CEED_QFUNCTION_HELPER void KMPack(const CeedScalar A[3][3], CeedScalar v[6]) { 156 const CeedScalar weight = sqrt(2.); 157 v[0] = A[0][0]; 158 v[1] = A[1][1]; 159 v[2] = A[2][2]; 160 v[3] = A[2][1] * weight; 161 v[4] = A[2][0] * weight; 162 v[5] = A[1][0] * weight; 163 } 164 165 // @brief Calculate metric tensor from mapping, g_{ij} = xi_{k,i} xi_{k,j} = dXdx^T dXdx 166 CEED_QFUNCTION_HELPER void KMMetricTensor(const CeedScalar dXdx[3][3], CeedScalar km_g_ij[6]) { 167 CeedScalar g_ij[3][3] = {{0.}}; 168 MatMat3(dXdx, dXdx, CEED_TRANSPOSE, CEED_NOTRANSPOSE, g_ij); 169 KMPack(g_ij, km_g_ij); 170 } 171 172 // @brief Linear ramp evaluation 173 CEED_QFUNCTION_HELPER CeedScalar LinearRampCoefficient(CeedScalar amplitude, CeedScalar length, CeedScalar start, CeedScalar x) { 174 if (x < start) { 175 return amplitude; 176 } else if (x < start + length) { 177 return amplitude * ((x - start) * (-1 / length) + 1); 178 } else { 179 return 0; 180 } 181 } 182 183 /** 184 @brief Pack stored values at quadrature point 185 186 @param[in] Q Number of quadrature points 187 @param[in] i Current quadrature point 188 @param[in] start Starting index to store components 189 @param[in] num_comp Number of components to store 190 @param[in] values_at_qpnt Local values for quadrature point i 191 @param[out] stored Stored values 192 193 @return An error code: 0 - success, otherwise - failure 194 **/ 195 CEED_QFUNCTION_HELPER int StoredValuesPack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *values_at_qpnt, 196 CeedScalar *stored) { 197 for (CeedInt j = 0; j < num_comp; j++) stored[(start + j) * Q + i] = values_at_qpnt[j]; 198 199 return CEED_ERROR_SUCCESS; 200 } 201 202 /** 203 @brief Unpack stored values at quadrature point 204 205 @param[in] Q Number of quadrature points 206 @param[in] i Current quadrature point 207 @param[in] start Starting index to store components 208 @param[in] num_comp Number of components to store 209 @param[in] stored Stored values 210 @param[out] values_at_qpnt Local values for quadrature point i 211 212 @return An error code: 0 - success, otherwise - failure 213 **/ 214 CEED_QFUNCTION_HELPER int StoredValuesUnpack(CeedInt Q, CeedInt i, CeedInt start, CeedInt num_comp, const CeedScalar *stored, 215 CeedScalar *values_at_qpnt) { 216 for (CeedInt j = 0; j < num_comp; j++) values_at_qpnt[j] = stored[(start + j) * Q + i]; 217 218 return CEED_ERROR_SUCCESS; 219 } 220 221 /** 222 @brief Unpack 3D element q_data at quadrature point 223 224 @param[in] Q Number of quadrature points 225 @param[in] i Current quadrature point 226 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`) 227 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian 228 @param[out] dXdx Inverse of the mapping Jacobian (shape [3][3]) 229 230 @return An error code: 0 - success, otherwise - failure 231 **/ 232 CEED_QFUNCTION_HELPER int QdataUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[3][3]) { 233 StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 234 StoredValuesUnpack(Q, i, 1, 9, q_data, (CeedScalar *)dXdx); 235 return CEED_ERROR_SUCCESS; 236 } 237 238 /** 239 @brief Unpack boundary element q_data for 3D problem at quadrature point 240 241 @param[in] Q Number of quadrature points 242 @param[in] i Current quadrature point 243 @param[in] q_data Pointer to q_data (generated by `setupgeo.h:Setup`) 244 @param[out] wdetJ Quadrature weight times determinant of the mapping Jacobian, or `NULL` 245 @param[out] dXdx Inverse of the mapping Jacobian (shape [2][3]), or `NULL` 246 @param[out] normal Components of the normal vector (shape [3]), or `NULL` 247 248 @return An error code: 0 - success, otherwise - failure 249 **/ 250 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_3D(CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar dXdx[2][3], 251 CeedScalar normal[3]) { 252 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 253 if (normal) StoredValuesUnpack(Q, i, 1, 3, q_data, normal); 254 if (dXdx) StoredValuesUnpack(Q, i, 4, 6, q_data, (CeedScalar *)dXdx); 255 return CEED_ERROR_SUCCESS; 256 } 257 258 #endif // utils_h 259