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