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/types.h> 7 8 #include "setupgeo_helpers.h" 9 #include "utils.h" 10 11 // ***************************************************************************** 12 // This QFunction sets up the geometric factors required for integration and coordinate transformations 13 // 14 // Reference (parent) coordinates: X 15 // Physical (current) coordinates: x 16 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 17 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 18 // 19 // All quadrature data is stored in 10 field vector of quadrature data. 20 // 21 // We require the determinant of the Jacobian to properly compute integrals of the form: int( v u ) 22 // 23 // Determinant of Jacobian: 24 // detJ = J11*A11 + J21*A12 + J31*A13 25 // Jij = Jacobian entry ij 26 // Aij = Adjugate ij 27 // 28 // Stored: w detJ 29 // in q_data[0] 30 // 31 // We require the transpose of the inverse of the Jacobian to properly compute integrals of the form: int( gradv u ) 32 // 33 // Inverse of Jacobian: 34 // dXdx_i,j = Aij / detJ 35 // 36 // Stored: Aij / detJ 37 // in q_data[1:9] as 38 // (detJ^-1) * [A11 A12 A13] 39 // [A21 A22 A23] 40 // [A31 A32 A33] 41 // ***************************************************************************** 42 CEED_QFUNCTION(Setup)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 43 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 44 const CeedScalar(*w) = in[1]; 45 CeedScalar(*q_data) = out[0]; 46 47 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 48 CeedScalar detJ, dXdx[3][3]; 49 InvertMappingJacobian_3D(Q, i, J, dXdx, &detJ); 50 const CeedScalar wdetJ = w[i] * detJ; 51 52 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data); 53 StoredValuesPack(Q, i, 1, 9, (const CeedScalar *)dXdx, q_data); 54 } 55 return 0; 56 } 57 58 // ***************************************************************************** 59 // This QFunction sets up the geometric factor required for integration when reference coordinates are in 2D and the physical coordinates are in 3D 60 // 61 // Reference (parent) 2D coordinates: X 62 // Physical (current) 3D coordinates: x 63 // Change of coordinate matrix: 64 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 65 // Inverse change of coordinate matrix: 66 // dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 67 // 68 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 69 // 70 // detJb is the magnitude of (J1,J2,J3) 71 // 72 // dXdx is calculated via Moore–Penrose inverse: 73 // 74 // dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 75 // = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 76 // 77 // All quadrature data is stored in 10 field vector of quadrature data. 78 // 79 // We require the determinant of the Jacobian to properly compute integrals of 80 // the form: int( u v ) 81 // 82 // Stored: w detJb 83 // in q_data_sur[0] 84 // 85 // Normal vector = (J1,J2,J3) / detJb 86 // 87 // - TODO Could possibly remove normal vector, as it could be calculated in the Qfunction from dXdx 88 // See https://github.com/CEED/libCEED/pull/868#discussion_r871979484 89 // Stored: (J1,J2,J3) / detJb 90 // in q_data_sur[1:3] as 91 // (detJb^-1) * [ J1 ] 92 // [ J2 ] 93 // [ J3 ] 94 // 95 // Stored: dXdx_{i,j} 96 // in q_data_sur[4:9] as 97 // [dXdx_11 dXdx_12 dXdx_13] 98 // [dXdx_21 dXdx_22 dXdx_23] 99 // ***************************************************************************** 100 CEED_QFUNCTION(SetupBoundary)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 101 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 102 const CeedScalar(*w) = in[1]; 103 CeedScalar(*q_data_sur) = out[0]; 104 105 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 106 CeedScalar detJb, normal[3], dXdx[2][3]; 107 108 NormalVectorFromdxdX_3D(Q, i, J, normal, &detJb); 109 InvertBoundaryMappingJacobian_3D(Q, i, J, dXdx); 110 const CeedScalar wdetJ = w[i] * detJb; 111 112 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur); 113 StoredValuesPack(Q, i, 1, 3, normal, q_data_sur); 114 StoredValuesPack(Q, i, 4, 6, (const CeedScalar *)dXdx, q_data_sur); 115 } 116 return 0; 117 } 118 119 /** 120 @brief Compute geometric factors for integration, gradient transformations, and coordinate transformations on element faces. 121 122 Reference (parent) 2D coordinates are given by `X` and physical (current) 3D coordinates are given by `x`. 123 The change of coordinate matrix is given by`dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2]`. 124 125 `(N_1, N_2, N_3)` is given by the cross product of the columns of `dxdX_{i,j}`. 126 127 `detNb` is the magnitude of `(N_1, N_2, N_3)`. 128 129 @param[in] ctx QFunction context, unused 130 @param[in] Q Number of quadrature points 131 @param[in] in Input arrays 132 - 0 - Jacobian of cell coordinates 133 - 1 - Jacobian of face coordinates 134 - 2 - quadrature weights 135 @param[out] out Output array 136 - 0 - qdata, `w detNb`, `dXdx`, and `N` 137 138 @return An error code: 0 - success, otherwise - failure 139 **/ 140 CEED_QFUNCTION(SetupBoundaryGradient)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 141 const CeedScalar(*J_cell)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 142 const CeedScalar(*J_face)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1]; 143 const CeedScalar(*w) = in[2]; 144 CeedScalar(*q_data_sur) = out[0]; 145 146 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 147 CeedScalar detJ_face, normal[3], dXdx[3][3]; 148 149 NormalVectorFromdxdX_3D(Q, i, J_face, normal, &detJ_face); 150 const CeedScalar wdetJ = w[i] * detJ_face; 151 InvertMappingJacobian_3D(Q, i, J_cell, dXdx, NULL); 152 153 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur); 154 StoredValuesPack(Q, i, 1, 9, (CeedScalar *)dXdx, q_data_sur); 155 StoredValuesPack(Q, i, 10, 3, normal, q_data_sur); 156 } 157 return CEED_ERROR_SUCCESS; 158 } 159