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 face element's normal vector from dxdX 77 * 78 * Reference (parent) 2D coordinates: X 79 * Physical (current) 3D coordinates: x 80 * Change of coordinate matrix: 81 * dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 82 * Inverse change of coordinate matrix: 83 * dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 84 * 85 * (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 86 * 87 * detJb is the magnitude of (J1,J2,J3) 88 * 89 * Normal vector = (J1,J2,J3) / detJb 90 * 91 * Stored: (J1,J2,J3) / detJb 92 * in q_data_sur[1:3] as 93 * (detJb^-1) * [ J1 ] 94 * [ J2 ] 95 * [ J3 ] 96 * 97 * @param[in] Q Number of quadrature points 98 * @param[in] i Current quadrature point 99 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 100 * @param[out] normal Inverse of mapping Jacobian at quadrature point i 101 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 102 */ 103 CEED_QFUNCTION_HELPER void NormalVectorFromdxdX_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar normal[3], 104 CeedScalar *detJ_ptr) { 105 const CeedScalar dxdX[3][2] = { 106 {dxdX_q[0][0][i], dxdX_q[1][0][i]}, 107 {dxdX_q[0][1][i], dxdX_q[1][1][i]}, 108 {dxdX_q[0][2][i], dxdX_q[1][2][i]} 109 }; 110 // J1, J2, and J3 are given by the cross product of the columns of dxdX 111 const CeedScalar J1 = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 112 const CeedScalar J2 = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 113 const CeedScalar J3 = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 114 115 const CeedScalar detJ = sqrt(J1 * J1 + J2 * J2 + J3 * J3); 116 117 normal[0] = J1 / detJ; 118 normal[1] = J2 / detJ; 119 normal[2] = J3 / detJ; 120 if (detJ_ptr) *detJ_ptr = detJ; 121 } 122 123 /** 124 * @brief Calculate inverse of mapping Jacobian, (dxdX)^-1 125 * 126 * Reference (parent) 2D coordinates: X 127 * Physical (current) 3D coordinates: x 128 * Change of coordinate matrix: 129 * dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 130 * Inverse change of coordinate matrix: 131 * dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 132 * 133 * dXdx is calculated via Moore–Penrose inverse: 134 * 135 * dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 136 * = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 137 * 138 * @param[in] Q Number of quadrature points 139 * @param[in] i Current quadrature point 140 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 141 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 142 */ 143 CEED_QFUNCTION_HELPER void InvertBoundaryMappingJacobian_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar dXdx[2][3]) { 144 const CeedScalar dxdX[3][2] = { 145 {dxdX_q[0][0][i], dxdX_q[1][0][i]}, 146 {dxdX_q[0][1][i], dxdX_q[1][1][i]}, 147 {dxdX_q[0][2][i], dxdX_q[1][2][i]} 148 }; 149 150 // dxdX_k,j * dxdX_j,k 151 CeedScalar dxdXTdxdX[2][2] = {{0.}}; 152 for (CeedInt j = 0; j < 2; j++) { 153 for (CeedInt k = 0; k < 2; k++) { 154 for (CeedInt l = 0; l < 3; l++) dxdXTdxdX[j][k] += dxdX[l][j] * dxdX[l][k]; 155 } 156 } 157 158 const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1] - dxdXTdxdX[1][0] * dxdXTdxdX[0][1]; 159 160 // Compute inverse of dxdXTdxdX 161 CeedScalar dxdXTdxdX_inv[2][2]; 162 dxdXTdxdX_inv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX; 163 dxdXTdxdX_inv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX; 164 dxdXTdxdX_inv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX; 165 dxdXTdxdX_inv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX; 166 167 // Compute dXdx from dxdXTdxdX^-1 and dxdX 168 for (CeedInt j = 0; j < 2; j++) { 169 for (CeedInt k = 0; k < 3; k++) { 170 dXdx[j][k] = 0; 171 for (CeedInt l = 0; l < 2; l++) dXdx[j][k] += dxdXTdxdX_inv[l][j] * dxdX[k][l]; 172 } 173 } 174 } 175 176 #endif // setupgeo_helpers_h 177