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 #pragma once 7 8 #include <ceed.h> 9 #include <math.h> 10 11 #include "utils.h" 12 13 /** 14 * @brief Calculate dXdx from dxdX for 3D elements 15 * 16 * Reference (parent) coordinates: X 17 * Physical (current) coordinates: x 18 * Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 19 * Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 20 * 21 * Determinant of Jacobian: 22 * detJ = J11*A11 + J21*A12 + J31*A13 23 * Jij = Jacobian entry ij 24 * Aij = Adjugate ij 25 * 26 * Inverse of Jacobian: 27 * dXdx_i,j = Aij / detJ 28 * 29 * @param[in] Q Number of quadrature points 30 * @param[in] i Current quadrature point 31 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 32 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 33 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 34 */ 35 CEED_QFUNCTION_HELPER void InvertMappingJacobian_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar dXdx[3][3], 36 CeedScalar *detJ_ptr) { 37 const CeedScalar dxdX_11 = dxdX_q[0][0][i]; 38 const CeedScalar dxdX_21 = dxdX_q[0][1][i]; 39 const CeedScalar dxdX_31 = dxdX_q[0][2][i]; 40 const CeedScalar dxdX_12 = dxdX_q[1][0][i]; 41 const CeedScalar dxdX_22 = dxdX_q[1][1][i]; 42 const CeedScalar dxdX_32 = dxdX_q[1][2][i]; 43 const CeedScalar dxdX_13 = dxdX_q[2][0][i]; 44 const CeedScalar dxdX_23 = dxdX_q[2][1][i]; 45 const CeedScalar dxdX_33 = dxdX_q[2][2][i]; 46 const CeedScalar A11 = dxdX_22 * dxdX_33 - dxdX_23 * dxdX_32; 47 const CeedScalar A12 = dxdX_13 * dxdX_32 - dxdX_12 * dxdX_33; 48 const CeedScalar A13 = dxdX_12 * dxdX_23 - dxdX_13 * dxdX_22; 49 const CeedScalar A21 = dxdX_23 * dxdX_31 - dxdX_21 * dxdX_33; 50 const CeedScalar A22 = dxdX_11 * dxdX_33 - dxdX_13 * dxdX_31; 51 const CeedScalar A23 = dxdX_13 * dxdX_21 - dxdX_11 * dxdX_23; 52 const CeedScalar A31 = dxdX_21 * dxdX_32 - dxdX_22 * dxdX_31; 53 const CeedScalar A32 = dxdX_12 * dxdX_31 - dxdX_11 * dxdX_32; 54 const CeedScalar A33 = dxdX_11 * dxdX_22 - dxdX_12 * dxdX_21; 55 const CeedScalar detJ = dxdX_11 * A11 + dxdX_21 * A12 + dxdX_31 * A13; 56 57 dXdx[0][0] = A11 / detJ; 58 dXdx[0][1] = A12 / detJ; 59 dXdx[0][2] = A13 / detJ; 60 dXdx[1][0] = A21 / detJ; 61 dXdx[1][1] = A22 / detJ; 62 dXdx[1][2] = A23 / detJ; 63 dXdx[2][0] = A31 / detJ; 64 dXdx[2][1] = A32 / detJ; 65 dXdx[2][2] = A33 / detJ; 66 if (detJ_ptr) *detJ_ptr = detJ; 67 } 68 69 /** 70 * @brief Calculate dXdx from dxdX for 3D elements 71 * 72 * Reference (parent) coordinates: X 73 * Physical (current) coordinates: x 74 * Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 75 * Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 76 * 77 * Determinant of Jacobian: 78 * detJ = J11*A11 + J21*A12 + J31*A13 79 * Jij = Jacobian entry ij 80 * Aij = Adjugate ij 81 * 82 * Inverse of Jacobian: 83 * dXdx_i,j = Aij / detJ 84 * 85 * @param[in] Q Number of quadrature points 86 * @param[in] i Current quadrature point 87 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 88 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 89 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 90 */ 91 CEED_QFUNCTION_HELPER void InvertMappingJacobian_2D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[2][CEED_Q_VLA], CeedScalar dXdx[2][2], 92 CeedScalar *detJ_ptr) { 93 const CeedScalar dxdX_11 = dxdX_q[0][0][i]; 94 const CeedScalar dxdX_21 = dxdX_q[0][1][i]; 95 const CeedScalar dxdX_12 = dxdX_q[1][0][i]; 96 const CeedScalar dxdX_22 = dxdX_q[1][1][i]; 97 const CeedScalar detJ = dxdX_11 * dxdX_22 - dxdX_21 * dxdX_12; 98 99 dXdx[0][0] = dxdX_22 / detJ; 100 dXdx[0][1] = -dxdX_12 / detJ; 101 dXdx[1][0] = -dxdX_21 / detJ; 102 dXdx[1][1] = dxdX_11 / detJ; 103 if (detJ_ptr) *detJ_ptr = detJ; 104 } 105 106 /** 107 * @brief Calculate face element's normal vector from dxdX 108 * 109 * Reference (parent) 2D coordinates: X 110 * Physical (current) 3D coordinates: x 111 * Change of coordinate matrix: 112 * dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 113 * Inverse change of coordinate matrix: 114 * dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 115 * 116 * (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 117 * 118 * detJb is the magnitude of (J1,J2,J3) 119 * 120 * Normal vector = (J1,J2,J3) / detJb 121 * 122 * @param[in] Q Number of quadrature points 123 * @param[in] i Current quadrature point 124 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 125 * @param[out] normal Inverse of mapping Jacobian at quadrature point i 126 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 127 */ 128 CEED_QFUNCTION_HELPER void NormalVectorFromdxdX_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar normal[3], 129 CeedScalar *detJ_ptr) { 130 const CeedScalar dxdX[3][2] = { 131 {dxdX_q[0][0][i], dxdX_q[1][0][i]}, 132 {dxdX_q[0][1][i], dxdX_q[1][1][i]}, 133 {dxdX_q[0][2][i], dxdX_q[1][2][i]} 134 }; 135 // J1, J2, and J3 are given by the cross product of the columns of dxdX 136 const CeedScalar J1 = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 137 const CeedScalar J2 = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 138 const CeedScalar J3 = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 139 140 const CeedScalar detJ = sqrt(J1 * J1 + J2 * J2 + J3 * J3); 141 142 normal[0] = J1 / detJ; 143 normal[1] = J2 / detJ; 144 normal[2] = J3 / detJ; 145 if (detJ_ptr) *detJ_ptr = detJ; 146 } 147 148 /** 149 * This QFunction sets up the geometric factor required for integration when reference coordinates are in 1D and the physical coordinates are in 2D 150 * 151 * Reference (parent) 1D coordinates: X 152 * Physical (current) 2D coordinates: x 153 * Change of coordinate vector: 154 * J1 = dx_1/dX 155 * J2 = dx_2/dX 156 * 157 * detJb is the magnitude of (J1,J2) 158 * 159 * We require the determinant of the Jacobian to properly compute integrals of the form: int( u v ) 160 * 161 * Normal vector is given by the cross product of (J1,J2)/detJ and ẑ 162 * 163 * @param[in] Q Number of quadrature points 164 * @param[in] i Current quadrature point 165 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 166 * @param[out] normal Inverse of mapping Jacobian at quadrature point i 167 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 168 */ 169 CEED_QFUNCTION_HELPER void NormalVectorFromdxdX_2D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[CEED_Q_VLA], CeedScalar normal[2], 170 CeedScalar *detJ_ptr) { 171 const CeedScalar J1 = dxdX_q[0][i]; 172 const CeedScalar J2 = dxdX_q[1][i]; 173 174 CeedScalar detJb = sqrt(J1 * J1 + J2 * J2); 175 normal[0] = J2 / detJb; 176 normal[1] = -J1 / detJb; 177 if (detJ_ptr) *detJ_ptr = detJb; 178 } 179 180 /** 181 * @brief Calculate inverse of mapping Jacobian, (dxdX)^-1 182 * 183 * Reference (parent) 2D coordinates: X 184 * Physical (current) 3D coordinates: x 185 * Change of coordinate matrix: 186 * dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 187 * Inverse change of coordinate matrix: 188 * dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 189 * 190 * dXdx is calculated via Moore–Penrose inverse: 191 * 192 * dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 193 * = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 194 * 195 * @param[in] Q Number of quadrature points 196 * @param[in] i Current quadrature point 197 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 198 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 199 */ 200 CEED_QFUNCTION_HELPER void InvertBoundaryMappingJacobian_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar dXdx[2][3]) { 201 const CeedScalar dxdX[3][2] = { 202 {dxdX_q[0][0][i], dxdX_q[1][0][i]}, 203 {dxdX_q[0][1][i], dxdX_q[1][1][i]}, 204 {dxdX_q[0][2][i], dxdX_q[1][2][i]} 205 }; 206 207 // dxdX_k,j * dxdX_j,k 208 CeedScalar dxdXTdxdX[2][2] = {{0.}}; 209 for (CeedInt j = 0; j < 2; j++) { 210 for (CeedInt k = 0; k < 2; k++) { 211 for (CeedInt l = 0; l < 3; l++) dxdXTdxdX[j][k] += dxdX[l][j] * dxdX[l][k]; 212 } 213 } 214 215 const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1] - dxdXTdxdX[1][0] * dxdXTdxdX[0][1]; 216 217 // Compute inverse of dxdXTdxdX 218 CeedScalar dxdXTdxdX_inv[2][2]; 219 dxdXTdxdX_inv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX; 220 dxdXTdxdX_inv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX; 221 dxdXTdxdX_inv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX; 222 dxdXTdxdX_inv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX; 223 224 // Compute dXdx from dxdXTdxdX^-1 and dxdX 225 for (CeedInt j = 0; j < 2; j++) { 226 for (CeedInt k = 0; k < 3; k++) { 227 dXdx[j][k] = 0; 228 for (CeedInt l = 0; l < 2; l++) dXdx[j][k] += dxdXTdxdX_inv[l][j] * dxdX[k][l]; 229 } 230 } 231 } 232