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