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