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