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