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