1 // Copyright (c) 2017-2025, 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/types.h> 13 #ifndef CEED_RUNNING_JIT_PASS 14 #include <math.h> 15 #endif 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 * @param[in] Q Number of quadrature points 129 * @param[in] i Current quadrature point 130 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 131 * @param[out] normal Inverse of mapping Jacobian at quadrature point i 132 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 133 */ 134 CEED_QFUNCTION_HELPER void NormalVectorFromdxdX_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar normal[3], 135 CeedScalar *detJ_ptr) { 136 const CeedScalar dxdX[3][2] = { 137 {dxdX_q[0][0][i], dxdX_q[1][0][i]}, 138 {dxdX_q[0][1][i], dxdX_q[1][1][i]}, 139 {dxdX_q[0][2][i], dxdX_q[1][2][i]} 140 }; 141 // J1, J2, and J3 are given by the cross product of the columns of dxdX 142 const CeedScalar J1 = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 143 const CeedScalar J2 = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 144 const CeedScalar J3 = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 145 146 const CeedScalar detJ = sqrt(J1 * J1 + J2 * J2 + J3 * J3); 147 148 normal[0] = J1 / detJ; 149 normal[1] = J2 / detJ; 150 normal[2] = J3 / detJ; 151 if (detJ_ptr) *detJ_ptr = detJ; 152 } 153 154 /** 155 * This QFunction sets up the geometric factor required for integration when reference coordinates are in 1D and the physical coordinates are in 2D 156 * 157 * Reference (parent) 1D coordinates: X 158 * Physical (current) 2D coordinates: x 159 * Change of coordinate vector: 160 * J1 = dx_1/dX 161 * J2 = dx_2/dX 162 * 163 * detJb is the magnitude of (J1,J2) 164 * 165 * We require the determinant of the Jacobian to properly compute integrals of the form: int( u v ) 166 * 167 * Normal vector is given by the cross product of (J1,J2)/detJ and ẑ 168 * 169 * @param[in] Q Number of quadrature points 170 * @param[in] i Current quadrature point 171 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 172 * @param[out] normal Inverse of mapping Jacobian at quadrature point i 173 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 174 */ 175 CEED_QFUNCTION_HELPER void NormalVectorFromdxdX_2D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[CEED_Q_VLA], CeedScalar normal[2], 176 CeedScalar *detJ_ptr) { 177 const CeedScalar J1 = dxdX_q[0][i]; 178 const CeedScalar J2 = dxdX_q[1][i]; 179 180 CeedScalar detJb = sqrt(J1 * J1 + J2 * J2); 181 normal[0] = J2 / detJb; 182 normal[1] = -J1 / detJb; 183 if (detJ_ptr) *detJ_ptr = detJb; 184 } 185 186 /** 187 * @brief Calculate inverse of mapping Jacobian, (dxdX)^-1 188 * 189 * Reference (parent) 2D coordinates: X 190 * Physical (current) 3D coordinates: x 191 * Change of coordinate matrix: 192 * dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 193 * Inverse change of coordinate matrix: 194 * dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 195 * 196 * dXdx is calculated via Moore–Penrose inverse: 197 * 198 * dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 199 * = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 200 * 201 * @param[in] Q Number of quadrature points 202 * @param[in] i Current quadrature point 203 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 204 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 205 */ 206 CEED_QFUNCTION_HELPER void InvertBoundaryMappingJacobian_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar dXdx[2][3]) { 207 const CeedScalar dxdX[3][2] = { 208 {dxdX_q[0][0][i], dxdX_q[1][0][i]}, 209 {dxdX_q[0][1][i], dxdX_q[1][1][i]}, 210 {dxdX_q[0][2][i], dxdX_q[1][2][i]} 211 }; 212 213 // dxdX_k,j * dxdX_j,k 214 CeedScalar dxdXTdxdX[2][2] = {{0.}}; 215 for (CeedInt j = 0; j < 2; j++) { 216 for (CeedInt k = 0; k < 2; k++) { 217 for (CeedInt l = 0; l < 3; l++) dxdXTdxdX[j][k] += dxdX[l][j] * dxdX[l][k]; 218 } 219 } 220 221 const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1] - dxdXTdxdX[1][0] * dxdXTdxdX[0][1]; 222 223 // Compute inverse of dxdXTdxdX 224 CeedScalar dxdXTdxdX_inv[2][2]; 225 dxdXTdxdX_inv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX; 226 dxdXTdxdX_inv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX; 227 dxdXTdxdX_inv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX; 228 dxdXTdxdX_inv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX; 229 230 // Compute dXdx from dxdXTdxdX^-1 and dxdX 231 for (CeedInt j = 0; j < 2; j++) { 232 for (CeedInt k = 0; k < 3; k++) { 233 dXdx[j][k] = 0; 234 for (CeedInt l = 0; l < 2; l++) dXdx[j][k] += dxdXTdxdX_inv[l][j] * dxdX[k][l]; 235 } 236 } 237 } 238