1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Geometric factors (3D) for Navier-Stokes example using PETSc 19 20 #ifndef setup_geo_h 21 #define setup_geo_h 22 23 #ifndef __CUDACC__ 24 # include <math.h> 25 #endif 26 27 // ***************************************************************************** 28 // This QFunction sets up the geometric factors required for integration and 29 // coordinate transformations 30 // 31 // Reference (parent) coordinates: X 32 // Physical (current) coordinates: x 33 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 34 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 35 // 36 // All quadrature data is stored in 10 field vector of quadrature data. 37 // 38 // We require the determinant of the Jacobian to properly compute integrals of 39 // the form: int( v u ) 40 // 41 // Determinant of Jacobian: 42 // detJ = J11*A11 + J21*A12 + J31*A13 43 // Jij = Jacobian entry ij 44 // Aij = Adjoint ij 45 // 46 // Stored: w detJ 47 // in q_data[0] 48 // 49 // We require the transpose of the inverse of the Jacobian to properly compute 50 // integrals of the form: int( gradv u ) 51 // 52 // Inverse of Jacobian: 53 // dXdx_i,j = Aij / detJ 54 // 55 // Stored: Aij / detJ 56 // in q_data[1:9] as 57 // (detJ^-1) * [A11 A12 A13] 58 // [A21 A22 A23] 59 // [A31 A32 A33] 60 // 61 // ***************************************************************************** 62 CEED_QFUNCTION(Setup)(void *ctx, CeedInt Q, 63 const CeedScalar *const *in, CeedScalar *const *out) { 64 // *INDENT-OFF* 65 // Inputs 66 const CeedScalar (*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 67 (*w) = in[1]; 68 69 // Outputs 70 CeedScalar (*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 71 // *INDENT-ON* 72 73 CeedPragmaSIMD 74 // Quadrature Point Loop 75 for (CeedInt i=0; i<Q; i++) { 76 // Setup 77 const CeedScalar J11 = J[0][0][i]; 78 const CeedScalar J21 = J[0][1][i]; 79 const CeedScalar J31 = J[0][2][i]; 80 const CeedScalar J12 = J[1][0][i]; 81 const CeedScalar J22 = J[1][1][i]; 82 const CeedScalar J32 = J[1][2][i]; 83 const CeedScalar J13 = J[2][0][i]; 84 const CeedScalar J23 = J[2][1][i]; 85 const CeedScalar J33 = J[2][2][i]; 86 const CeedScalar A11 = J22*J33 - J23*J32; 87 const CeedScalar A12 = J13*J32 - J12*J33; 88 const CeedScalar A13 = J12*J23 - J13*J22; 89 const CeedScalar A21 = J23*J31 - J21*J33; 90 const CeedScalar A22 = J11*J33 - J13*J31; 91 const CeedScalar A23 = J13*J21 - J11*J23; 92 const CeedScalar A31 = J21*J32 - J22*J31; 93 const CeedScalar A32 = J12*J31 - J11*J32; 94 const CeedScalar A33 = J11*J22 - J12*J21; 95 const CeedScalar detJ = J11*A11 + J21*A12 + J31*A13; 96 97 // Qdata 98 // -- Interp-to-Interp q_data 99 q_data[0][i] = w[i] * detJ; 100 // -- Interp-to-Grad q_data 101 // Inverse of change of coordinate matrix: X_i,j 102 q_data[1][i] = A11 / detJ; 103 q_data[2][i] = A12 / detJ; 104 q_data[3][i] = A13 / detJ; 105 q_data[4][i] = A21 / detJ; 106 q_data[5][i] = A22 / detJ; 107 q_data[6][i] = A23 / detJ; 108 q_data[7][i] = A31 / detJ; 109 q_data[8][i] = A32 / detJ; 110 q_data[9][i] = A33 / detJ; 111 112 } // End of Quadrature Point Loop 113 114 // Return 115 return 0; 116 } 117 118 // ***************************************************************************** 119 // This QFunction sets up the geometric factor required for integration when 120 // reference coordinates are in 2D and the physical coordinates are in 3D 121 // 122 // Reference (parent) 2D coordinates: X 123 // Physical (current) 3D coordinates: x 124 // Change of coordinate matrix: 125 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 126 // 127 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 128 // 129 // detJb is the magnitude of (J1,J2,J3) 130 // 131 // All quadrature data is stored in 4 field vector of quadrature data. 132 // 133 // We require the determinant of the Jacobian to properly compute integrals of 134 // the form: int( u v ) 135 // 136 // Stored: w detJb 137 // in q_data_sur[0] 138 // 139 // Normal vector = (J1,J2,J3) / detJb 140 // 141 // Stored: (J1,J2,J3) / detJb 142 // in q_data_sur[1:3] as 143 // (detJb^-1) * [ J1 ] 144 // [ J2 ] 145 // [ J3 ] 146 // 147 // ***************************************************************************** 148 CEED_QFUNCTION(SetupBoundary)(void *ctx, CeedInt Q, 149 const CeedScalar *const *in, CeedScalar *const *out) { 150 // *INDENT-OFF* 151 // Inputs 152 const CeedScalar (*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 153 (*w) = in[1]; 154 // Outputs 155 CeedScalar (*q_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 156 157 CeedPragmaSIMD 158 // Quadrature Point Loop 159 for (CeedInt i=0; i<Q; i++) { 160 // Setup 161 const CeedScalar dxdX[3][2] = {{J[0][0][i], 162 J[1][0][i]}, 163 {J[0][1][i], 164 J[1][1][i]}, 165 {J[0][2][i], 166 J[1][2][i]} 167 }; 168 // *INDENT-ON* 169 // J1, J2, and J3 are given by the cross product of the columns of dxdX 170 const CeedScalar J1 = dxdX[1][0]*dxdX[2][1] - dxdX[2][0]*dxdX[1][1]; 171 const CeedScalar J2 = dxdX[2][0]*dxdX[0][1] - dxdX[0][0]*dxdX[2][1]; 172 const CeedScalar J3 = dxdX[0][0]*dxdX[1][1] - dxdX[1][0]*dxdX[0][1]; 173 174 const CeedScalar detJb = sqrt(J1*J1 + J2*J2 + J3*J3); 175 176 // q_data_sur 177 // -- Interp-to-Interp q_data_sur 178 q_data_sur[0][i] = w[i] * detJb; 179 q_data_sur[1][i] = J1 / detJb; 180 q_data_sur[2][i] = J2 / detJb; 181 q_data_sur[3][i] = J3 / detJb; 182 183 } // End of Quadrature Point Loop 184 185 // Return 186 return 0; 187 } 188 189 // ***************************************************************************** 190 191 #endif // setup_geo_h 192