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