1 // Copyright (c) 2017-2024, 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 #include <ceed.h> 11 #include "setupgeo_helpers.h" 12 #include "utils.h" 13 14 // ***************************************************************************** 15 // This QFunction sets up the geometric factors required for integration and coordinate transformations 16 // 17 // Reference (parent) coordinates: X 18 // Physical (current) coordinates: x 19 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 20 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 21 // 22 // All quadrature data is stored in 10 field vector of quadrature data. 23 // 24 // We require the determinant of the Jacobian to properly compute integrals of the form: int( v u ) 25 // 26 // Determinant of Jacobian: 27 // detJ = J11*J22 - J21*J12 28 // Jij = Jacobian entry ij 29 // 30 // Stored: w detJ 31 // in q_data[0] 32 // 33 // We require the transpose of the inverse of the Jacobian to properly compute integrals of the form: int( gradv u ) 34 // 35 // Inverse of Jacobian: 36 // dXdx_i,j = Aij / detJ 37 // Aij = Adjugate ij 38 // 39 // Stored: Aij / detJ 40 // in q_data[1:4] as 41 // (detJ^-1) * [A11 A12] 42 // [A21 A22] 43 // ***************************************************************************** 44 CEED_QFUNCTION(Setup2d)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 45 const CeedScalar(*J)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[0]; 46 const CeedScalar(*w) = in[1]; 47 CeedScalar(*q_data) = out[0]; 48 49 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 50 CeedScalar dXdx[2][2], detJ; 51 InvertMappingJacobian_2D(Q, i, J, dXdx, &detJ); 52 const CeedScalar wdetJ = w[i] * detJ; 53 54 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data); 55 StoredValuesPack(Q, i, 1, 4, (const CeedScalar *)dXdx, q_data); 56 } 57 return 0; 58 } 59 60 // ***************************************************************************** 61 // This QFunction sets up the geometric factor required for integration when reference coordinates are in 1D and the physical coordinates are in 2D 62 // 63 // Reference (parent) 1D coordinates: X 64 // Physical (current) 2D coordinates: x 65 // Change of coordinate vector: 66 // J1 = dx_1/dX 67 // J2 = dx_2/dX 68 // 69 // detJb is the magnitude of (J1,J2) 70 // 71 // All quadrature data is stored in 3 field vector of quadrature data. 72 // 73 // We require the determinant of the Jacobian to properly compute integrals of the form: int( u v ) 74 // 75 // Stored: w detJb 76 // in q_data_sur[0] 77 // 78 // Normal vector is given by the cross product of (J1,J2)/detJ and ẑ 79 // 80 // Stored: (J1,J2,0) x (0,0,1) / detJb 81 // in q_data_sur[1:2] as 82 // (detJb^-1) * [ J2 ] 83 // [-J1 ] 84 // ***************************************************************************** 85 CEED_QFUNCTION(SetupBoundary2d)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 86 const CeedScalar(*J)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 87 const CeedScalar(*w) = in[1]; 88 CeedScalar(*q_data_sur) = out[0]; 89 90 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 91 CeedScalar normal[2], detJb; 92 NormalVectorFromdxdX_2D(Q, i, J, normal, &detJb); 93 const CeedScalar wdetJ = w[i] * detJb; 94 95 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur); 96 StoredValuesPack(Q, i, 1, 2, normal, q_data_sur); 97 } 98 return 0; 99 } 100 101 // ***************************************************************************** 102 // This QFunction sets up the geometric factor required for integration when reference coordinates are in 2D and the physical coordinates are in 3D 103 // 104 // Reference (parent) 2D coordinates: X 105 // Physical (current) 3D coordinates: x 106 // Change of coordinate matrix: 107 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 108 // Inverse change of coordinate matrix: 109 // dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 110 // 111 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 112 // 113 // detJb is the magnitude of (J1,J2,J3) 114 // 115 // dXdx is calculated via Moore–Penrose inverse: 116 // 117 // dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 118 // = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 119 // 120 // All quadrature data is stored in 10 field vector of quadrature data. 121 // 122 // We require the determinant of the Jacobian to properly compute integrals of 123 // the form: int( u v ) 124 // 125 // Stored: w detJb 126 // in q_data_sur[0] 127 // 128 // Normal vector = (J1,J2,J3) / detJb 129 // 130 // Stored: (J1,J2,J3) / detJb 131 // 132 // Stored: dXdx_{i,j} 133 // in q_data_sur[1:6] as 134 // [dXdx_11 dXdx_12 dXdx_13] 135 // [dXdx_21 dXdx_22 dXdx_23] 136 // ***************************************************************************** 137 CEED_QFUNCTION(Setup2D_3Dcoords)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 138 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 139 const CeedScalar(*w) = in[1]; 140 CeedScalar(*q_data_sur) = out[0]; 141 142 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 143 CeedScalar detJb, normal[3], dXdx[2][3]; 144 145 NormalVectorFromdxdX_3D(Q, i, J, normal, &detJb); 146 InvertBoundaryMappingJacobian_3D(Q, i, J, dXdx); 147 const CeedScalar wdetJ = w[i] * detJb; 148 149 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur); 150 StoredValuesPack(Q, i, 1, 6, (const CeedScalar *)dXdx, q_data_sur); 151 } 152 return 0; 153 } 154