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 /// @file 9 /// Geometric factors (3D) for Navier-Stokes example using PETSc 10 11 #ifndef setup_geo_h 12 #define setup_geo_h 13 14 #include <ceed.h> 15 #include <math.h> 16 17 // ***************************************************************************** 18 // This QFunction sets up the geometric factors required for integration and coordinate transformations 19 // 20 // Reference (parent) coordinates: X 21 // Physical (current) coordinates: x 22 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 23 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 24 // 25 // All quadrature data is stored in 10 field vector of quadrature data. 26 // 27 // We require the determinant of the Jacobian to properly compute integrals of the form: int( v u ) 28 // 29 // Determinant of Jacobian: 30 // detJ = J11*A11 + J21*A12 + J31*A13 31 // Jij = Jacobian entry ij 32 // Aij = Adjoint ij 33 // 34 // Stored: w detJ 35 // in q_data[0] 36 // 37 // We require the transpose of the inverse of the Jacobian to properly compute integrals of the form: int( gradv u ) 38 // 39 // Inverse of Jacobian: 40 // dXdx_i,j = Aij / detJ 41 // 42 // Stored: Aij / detJ 43 // in q_data[1:9] as 44 // (detJ^-1) * [A11 A12 A13] 45 // [A21 A22 A23] 46 // [A31 A32 A33] 47 // ***************************************************************************** 48 CEED_QFUNCTION(Setup)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 49 // Inputs 50 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 51 const CeedScalar(*w) = in[1]; 52 53 // Outputs 54 CeedScalar(*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 55 56 CeedPragmaSIMD 57 // Quadrature Point Loop 58 for (CeedInt i = 0; i < Q; i++) { 59 // Setup 60 const CeedScalar J11 = J[0][0][i]; 61 const CeedScalar J21 = J[0][1][i]; 62 const CeedScalar J31 = J[0][2][i]; 63 const CeedScalar J12 = J[1][0][i]; 64 const CeedScalar J22 = J[1][1][i]; 65 const CeedScalar J32 = J[1][2][i]; 66 const CeedScalar J13 = J[2][0][i]; 67 const CeedScalar J23 = J[2][1][i]; 68 const CeedScalar J33 = J[2][2][i]; 69 const CeedScalar A11 = J22 * J33 - J23 * J32; 70 const CeedScalar A12 = J13 * J32 - J12 * J33; 71 const CeedScalar A13 = J12 * J23 - J13 * J22; 72 const CeedScalar A21 = J23 * J31 - J21 * J33; 73 const CeedScalar A22 = J11 * J33 - J13 * J31; 74 const CeedScalar A23 = J13 * J21 - J11 * J23; 75 const CeedScalar A31 = J21 * J32 - J22 * J31; 76 const CeedScalar A32 = J12 * J31 - J11 * J32; 77 const CeedScalar A33 = J11 * J22 - J12 * J21; 78 const CeedScalar detJ = J11 * A11 + J21 * A12 + J31 * A13; 79 80 // Qdata 81 // -- Interp-to-Interp q_data 82 q_data[0][i] = w[i] * detJ; 83 // -- Interp-to-Grad q_data 84 // Inverse of change of coordinate matrix: X_i,j 85 q_data[1][i] = A11 / detJ; 86 q_data[2][i] = A12 / detJ; 87 q_data[3][i] = A13 / detJ; 88 q_data[4][i] = A21 / detJ; 89 q_data[5][i] = A22 / detJ; 90 q_data[6][i] = A23 / detJ; 91 q_data[7][i] = A31 / detJ; 92 q_data[8][i] = A32 / detJ; 93 q_data[9][i] = A33 / detJ; 94 95 } // End of Quadrature Point Loop 96 97 // Return 98 return 0; 99 } 100 101 // ***************************************************************************** 102 // This QFunction sets up the geometric factor required for integration when reference coordinates are in 2D and the physical coordinates are in 3D 103 // 104 // Reference (parent) 2D coordinates: X 105 // Physical (current) 3D coordinates: x 106 // Change of coordinate matrix: 107 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 108 // Inverse change of coordinate matrix: 109 // dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 110 // 111 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 112 // 113 // detJb is the magnitude of (J1,J2,J3) 114 // 115 // dXdx is calculated via Moore–Penrose inverse: 116 // 117 // dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 118 // = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 119 // 120 // All quadrature data is stored in 10 field vector of quadrature data. 121 // 122 // We require the determinant of the Jacobian to properly compute integrals of 123 // the form: int( u v ) 124 // 125 // Stored: w detJb 126 // in q_data_sur[0] 127 // 128 // Normal vector = (J1,J2,J3) / detJb 129 // 130 // - TODO Could possibly remove normal vector, as it could be calculated in the Qfunction from dXdx 131 // Stored: (J1,J2,J3) / detJb 132 // in q_data_sur[1:3] as 133 // (detJb^-1) * [ J1 ] 134 // [ J2 ] 135 // [ J3 ] 136 // 137 // Stored: dXdx_{i,j} 138 // in q_data_sur[4:9] as 139 // [dXdx_11 dXdx_12 dXdx_13] 140 // [dXdx_21 dXdx_22 dXdx_23] 141 // ***************************************************************************** 142 CEED_QFUNCTION(SetupBoundary)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 143 // Inputs 144 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 145 const CeedScalar(*w) = in[1]; 146 147 // Outputs 148 CeedScalar(*q_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 149 150 CeedPragmaSIMD 151 // Quadrature Point Loop 152 for (CeedInt i = 0; i < Q; i++) { 153 // Setup 154 const CeedScalar dxdX[3][2] = { 155 {J[0][0][i], J[1][0][i]}, 156 {J[0][1][i], J[1][1][i]}, 157 {J[0][2][i], J[1][2][i]} 158 }; 159 // J1, J2, and J3 are given by the cross product of the columns of dxdX 160 const CeedScalar J1 = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 161 const CeedScalar J2 = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 162 const CeedScalar J3 = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 163 164 const CeedScalar detJb = sqrt(J1 * J1 + J2 * J2 + J3 * J3); 165 166 // q_data_sur 167 // -- Interp-to-Interp q_data_sur 168 q_data_sur[0][i] = w[i] * detJb; 169 q_data_sur[1][i] = J1 / detJb; 170 q_data_sur[2][i] = J2 / detJb; 171 q_data_sur[3][i] = J3 / detJb; 172 173 // dxdX_k,j * dxdX_j,k 174 CeedScalar dxdXTdxdX[2][2] = {{0.}}; 175 for (CeedInt j = 0; j < 2; j++) { 176 for (CeedInt k = 0; k < 2; k++) { 177 for (CeedInt l = 0; l < 3; l++) dxdXTdxdX[j][k] += dxdX[l][j] * dxdX[l][k]; 178 } 179 } 180 181 const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1] - dxdXTdxdX[1][0] * dxdXTdxdX[0][1]; 182 183 // Compute inverse of dxdXTdxdX 184 CeedScalar dxdXTdxdX_inv[2][2]; 185 dxdXTdxdX_inv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX; 186 dxdXTdxdX_inv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX; 187 dxdXTdxdX_inv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX; 188 dxdXTdxdX_inv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX; 189 190 // Compute dXdx from dxdXTdxdX^-1 and dxdX 191 CeedScalar dXdx[2][3] = {{0.}}; 192 for (CeedInt j = 0; j < 2; j++) { 193 for (CeedInt k = 0; k < 3; k++) { 194 for (CeedInt l = 0; l < 2; l++) dXdx[j][k] += dxdXTdxdX_inv[l][j] * dxdX[k][l]; 195 } 196 } 197 198 q_data_sur[4][i] = dXdx[0][0]; 199 q_data_sur[5][i] = dXdx[0][1]; 200 q_data_sur[6][i] = dXdx[0][2]; 201 q_data_sur[7][i] = dXdx[1][0]; 202 q_data_sur[8][i] = dXdx[1][1]; 203 q_data_sur[9][i] = dXdx[1][2]; 204 205 } // End of Quadrature Point Loop 206 207 // Return 208 return 0; 209 } 210 211 // ***************************************************************************** 212 213 #endif // setup_geo_h 214