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 (2D) for Navier-Stokes example using PETSc 19 20 #ifndef setup_geo_2d_h 21 #define setup_geo_2d_h 22 23 #include <math.h> 24 25 // ***************************************************************************** 26 // This QFunction sets up the geometric factors required for integration and 27 // coordinate transformations 28 // 29 // Reference (parent) coordinates: X 30 // Physical (current) coordinates: x 31 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 32 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 33 // 34 // All quadrature data is stored in 10 field vector of quadrature data. 35 // 36 // We require the determinant of the Jacobian to properly compute integrals of 37 // the form: int( v u ) 38 // 39 // Determinant of Jacobian: 40 // detJ = J11*J22 - J21*J12 41 // Jij = Jacobian entry ij 42 // 43 // Stored: w detJ 44 // in q_data[0] 45 // 46 // We require the transpose of the inverse of the Jacobian to properly compute 47 // integrals of the form: int( gradv u ) 48 // 49 // Inverse of Jacobian: 50 // dXdx_i,j = Aij / detJ 51 // Aij = Adjoint ij 52 // 53 // Stored: Aij / detJ 54 // in q_data[1:4] as 55 // (detJ^-1) * [A11 A12] 56 // [A21 A22] 57 // 58 // ***************************************************************************** 59 CEED_QFUNCTION(Setup2d)(void *ctx, CeedInt Q, 60 const CeedScalar *const *in, CeedScalar *const *out) { 61 // *INDENT-OFF* 62 // Inputs 63 const CeedScalar (*J)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[0], 64 (*w) = in[1]; 65 // Outputs 66 CeedScalar (*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 67 // *INDENT-ON* 68 69 CeedPragmaSIMD 70 // Quadrature Point Loop 71 for (CeedInt i=0; i<Q; i++) { 72 // Setup 73 const CeedScalar J11 = J[0][0][i]; 74 const CeedScalar J21 = J[0][1][i]; 75 const CeedScalar J12 = J[1][0][i]; 76 const CeedScalar J22 = J[1][1][i]; 77 const CeedScalar detJ = J11*J22 - J21*J12; 78 79 // Qdata 80 // -- Interp-to-Interp q_data 81 q_data[0][i] = w[i] * detJ; 82 // -- Interp-to-Grad q_data 83 // Inverse of change of coordinate matrix: X_i,j 84 q_data[1][i] = J22 / detJ; 85 q_data[2][i] = -J21 / detJ; 86 q_data[3][i] = -J12 / detJ; 87 q_data[4][i] = J11 / detJ; 88 } // End of Quadrature Point Loop 89 90 // Return 91 return 0; 92 } 93 94 // ***************************************************************************** 95 // This QFunction sets up the geometric factor required for integration when 96 // reference coordinates are in 1D and the physical coordinates are in 2D 97 // 98 // Reference (parent) 1D coordinates: X 99 // Physical (current) 2D coordinates: x 100 // Change of coordinate vector: 101 // J1 = dx_1/dX 102 // J2 = dx_2/dX 103 // 104 // detJb is the magnitude of (J1,J2) 105 // 106 // All quadrature data is stored in 3 field vector of quadrature data. 107 // 108 // We require the determinant of the Jacobian to properly compute integrals of 109 // the form: int( u v ) 110 // 111 // Stored: w detJb 112 // in q_data_sur[0] 113 // 114 // Normal vector is given by the cross product of (J1,J2)/detJ and ẑ 115 // 116 // Stored: (J1,J2,0) x (0,0,1) / detJb 117 // in q_data_sur[1:2] as 118 // (detJb^-1) * [ J2 ] 119 // [-J1 ] 120 // 121 // ***************************************************************************** 122 CEED_QFUNCTION(SetupBoundary2d)(void *ctx, CeedInt Q, 123 const CeedScalar *const *in, CeedScalar *const *out) { 124 // *INDENT-OFF* 125 // Inputs 126 const CeedScalar (*J)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 127 (*w) = in[1]; 128 // Outputs 129 CeedScalar (*q_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 130 // *INDENT-ON* 131 132 CeedPragmaSIMD 133 // Quadrature Point Loop 134 for (CeedInt i=0; i<Q; i++) { 135 // Setup 136 const CeedScalar J1 = J[0][i]; 137 const CeedScalar J2 = J[1][i]; 138 139 const CeedScalar detJb = sqrt(J1*J1 + J2*J2); 140 141 q_data_sur[0][i] = w[i] * detJb; 142 q_data_sur[1][i] = J2 / detJb; 143 q_data_sur[2][i] = -J1 / detJb; 144 } // End of Quadrature Point Loop 145 146 // Return 147 return 0; 148 } 149 150 // ***************************************************************************** 151 152 #endif // setup_geo_2d_h 153