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