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