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