1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 /// @file 5 /// Geometric factors (3D) for Navier-Stokes example using PETSc 6 #include <ceed.h> 7 #include <math.h> 8 9 #include "setupgeo_helpers.h" 10 #include "utils.h" 11 12 // ***************************************************************************** 13 // This QFunction sets up the geometric factors required for integration and coordinate transformations 14 // 15 // Reference (parent) coordinates: X 16 // Physical (current) coordinates: x 17 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 18 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 19 // 20 // All quadrature data is stored in 10 field vector of quadrature data. 21 // 22 // We require the determinant of the Jacobian to properly compute integrals of the form: int( v u ) 23 // 24 // Determinant of Jacobian: 25 // detJ = J11*A11 + J21*A12 + J31*A13 26 // Jij = Jacobian entry ij 27 // Aij = Adjugate ij 28 // 29 // Stored: w detJ 30 // in q_data[0] 31 // 32 // We require the transpose of the inverse of the Jacobian to properly compute integrals of the form: int( gradv u ) 33 // 34 // Inverse of Jacobian: 35 // dXdx_i,j = Aij / detJ 36 // 37 // Stored: Aij / detJ 38 // in q_data[1:9] as 39 // (detJ^-1) * [A11 A12 A13] 40 // [A21 A22 A23] 41 // [A31 A32 A33] 42 // ***************************************************************************** 43 CEED_QFUNCTION(Setup)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 44 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 45 const CeedScalar(*w) = in[1]; 46 CeedScalar(*q_data) = out[0]; 47 48 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 49 CeedScalar detJ, dXdx[3][3]; 50 InvertMappingJacobian_3D(Q, i, J, dXdx, &detJ); 51 const CeedScalar wdetJ = w[i] * detJ; 52 53 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data); 54 StoredValuesPack(Q, i, 1, 9, (const CeedScalar *)dXdx, q_data); 55 } 56 return 0; 57 } 58 59 // ***************************************************************************** 60 // This QFunction sets up the geometric factor required for integration when reference coordinates are in 2D and the physical coordinates are in 3D 61 // 62 // Reference (parent) 2D coordinates: X 63 // Physical (current) 3D coordinates: x 64 // Change of coordinate matrix: 65 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 66 // Inverse change of coordinate matrix: 67 // dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 68 // 69 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 70 // 71 // detJb is the magnitude of (J1,J2,J3) 72 // 73 // dXdx is calculated via Moore–Penrose inverse: 74 // 75 // dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 76 // = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 77 // 78 // All quadrature data is stored in 10 field vector of quadrature data. 79 // 80 // We require the determinant of the Jacobian to properly compute integrals of 81 // the form: int( u v ) 82 // 83 // Stored: w detJb 84 // in q_data_sur[0] 85 // 86 // Normal vector = (J1,J2,J3) / detJb 87 // 88 // - TODO Could possibly remove normal vector, as it could be calculated in the Qfunction from dXdx 89 // See https://github.com/CEED/libCEED/pull/868#discussion_r871979484 90 // Stored: (J1,J2,J3) / detJb 91 // in q_data_sur[1:3] as 92 // (detJb^-1) * [ J1 ] 93 // [ J2 ] 94 // [ J3 ] 95 // 96 // Stored: dXdx_{i,j} 97 // in q_data_sur[4:9] as 98 // [dXdx_11 dXdx_12 dXdx_13] 99 // [dXdx_21 dXdx_22 dXdx_23] 100 // ***************************************************************************** 101 CEED_QFUNCTION(SetupBoundary)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 102 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 103 const CeedScalar(*w) = in[1]; 104 CeedScalar(*q_data_sur) = out[0]; 105 106 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 107 CeedScalar detJb, normal[3], dXdx[2][3]; 108 109 NormalVectorFromdxdX_3D(Q, i, J, normal, &detJb); 110 InvertBoundaryMappingJacobian_3D(Q, i, J, dXdx); 111 const CeedScalar wdetJ = w[i] * detJb; 112 113 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur); 114 StoredValuesPack(Q, i, 1, 3, normal, q_data_sur); 115 StoredValuesPack(Q, i, 4, 6, (const CeedScalar *)dXdx, q_data_sur); 116 } 117 return 0; 118 } 119