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 CeedScalar dxdX[3][3]; 38 39 GradUnpack3(Q, i, 3, (CeedScalar *)dxdX_q, dxdX); 40 // Compute Adjugate of dxdX 41 dXdx[0][0] = dxdX[1][1] * dxdX[2][2] - dxdX[1][2] * dxdX[2][1]; 42 dXdx[0][1] = dxdX[0][2] * dxdX[2][1] - dxdX[0][1] * dxdX[2][2]; 43 dXdx[0][2] = dxdX[0][1] * dxdX[1][2] - dxdX[0][2] * dxdX[1][1]; 44 dXdx[1][0] = dxdX[1][2] * dxdX[2][0] - dxdX[1][0] * dxdX[2][2]; 45 dXdx[1][1] = dxdX[0][0] * dxdX[2][2] - dxdX[0][2] * dxdX[2][0]; 46 dXdx[1][2] = dxdX[0][2] * dxdX[1][0] - dxdX[0][0] * dxdX[1][2]; 47 dXdx[2][0] = dxdX[1][0] * dxdX[2][1] - dxdX[1][1] * dxdX[2][0]; 48 dXdx[2][1] = dxdX[0][1] * dxdX[2][0] - dxdX[0][0] * dxdX[2][1]; 49 dXdx[2][2] = dxdX[0][0] * dxdX[1][1] - dxdX[0][1] * dxdX[1][0]; 50 51 const CeedScalar detJ = dxdX[0][0] * dXdx[0][0] + dxdX[1][0] * dXdx[0][1] + dxdX[2][0] * dXdx[0][2]; 52 ScaleN((CeedScalar *)dXdx, 1 / detJ, 9); 53 if (detJ_ptr) *detJ_ptr = detJ; 54 } 55 56 /** 57 * @brief Calculate dXdx from dxdX for 2D elements 58 * 59 * Reference (parent) coordinates: X 60 * Physical (current) coordinates: x 61 * Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 62 * Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 63 * 64 * Determinant of Jacobian: 65 * detJ = J11*J22 - J21*J12 66 * Jij = Jacobian entry ij 67 * Aij = Adjugate ij 68 * 69 * Inverse of Jacobian: 70 * dXdx_i,j = Aij / detJ 71 * 72 * @param[in] Q Number of quadrature points 73 * @param[in] i Current quadrature point 74 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 75 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 76 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 77 */ 78 CEED_QFUNCTION_HELPER void InvertMappingJacobian_2D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[2][CEED_Q_VLA], CeedScalar dXdx[2][2], 79 CeedScalar *detJ_ptr) { 80 CeedScalar dxdX[2][2]; 81 82 GradUnpack2(Q, i, 2, (CeedScalar *)dxdX_q, dxdX); 83 const CeedScalar detJ = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 84 85 dXdx[0][0] = dxdX[1][1] / detJ; 86 dXdx[0][1] = -dxdX[0][1] / detJ; 87 dXdx[1][0] = -dxdX[1][0] / detJ; 88 dXdx[1][1] = dxdX[0][0] / detJ; 89 if (detJ_ptr) *detJ_ptr = detJ; 90 } 91 92 /** 93 * @brief Calculate face element's normal vector from dxdX 94 * 95 * Reference (parent) 2D coordinates: X 96 * Physical (current) 3D coordinates: x 97 * Change of coordinate matrix: 98 * dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 99 * Inverse change of coordinate matrix: 100 * dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 101 * 102 * (N1,N2,N3) is given by the cross product of the columns of dxdX_{i,j} 103 * 104 * detJb is the magnitude of (N1,N2,N3) 105 * 106 * Normal vector = (N1,N2,N3) / detJb 107 * 108 * @param[in] Q Number of quadrature points 109 * @param[in] i Current quadrature point 110 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 111 * @param[out] normal Inverse of mapping Jacobian at quadrature point i 112 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 113 */ 114 CEED_QFUNCTION_HELPER void NormalVectorFromdxdX_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar normal[3], 115 CeedScalar *detJ_ptr) { 116 CeedScalar dxdX[3][2]; 117 118 GradUnpack2(Q, i, 3, (CeedScalar *)dxdX_q, dxdX); 119 // N1, N2, and N3 are given by the cross product of the columns of dxdX 120 normal[0] = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 121 normal[1] = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 122 normal[2] = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 123 124 const CeedScalar detJ = Norm3(normal); 125 ScaleN(normal, 1 / detJ, 3); 126 if (detJ_ptr) *detJ_ptr = detJ; 127 } 128 129 /** 130 * This QFunction sets up the geometric factor required for integration when reference coordinates are in 1D and the physical coordinates are in 2D 131 * 132 * Reference (parent) 1D coordinates: X 133 * Physical (current) 2D coordinates: x 134 * Change of coordinate vector: 135 * N1 = dx_1/dX 136 * N2 = dx_2/dX 137 * 138 * detJb is the magnitude of (N1,N2) 139 * 140 * We require the determinant of the Jacobian to properly compute integrals of the form: int( u v ) 141 * 142 * Normal vector is given by the cross product of (N1,N2)/detJ and ẑ 143 * 144 * @param[in] Q Number of quadrature points 145 * @param[in] i Current quadrature point 146 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 147 * @param[out] normal Inverse of mapping Jacobian at quadrature point i 148 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 149 */ 150 CEED_QFUNCTION_HELPER void NormalVectorFromdxdX_2D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[CEED_Q_VLA], CeedScalar normal[2], 151 CeedScalar *detJ_ptr) { 152 normal[0] = dxdX_q[1][i]; 153 normal[1] = -dxdX_q[0][i]; 154 const CeedScalar detJb = Norm2(normal); 155 ScaleN(normal, 1 / detJb, 2); 156 if (detJ_ptr) *detJ_ptr = detJb; 157 } 158 159 /** 160 * @brief Calculate inverse of mapping Jacobian, (dxdX)^-1 161 * 162 * Reference (parent) 2D coordinates: X 163 * Physical (current) 3D coordinates: x 164 * Change of coordinate matrix: 165 * dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 166 * Inverse change of coordinate matrix: 167 * dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 168 * 169 * dXdx is calculated via Moore–Penrose inverse: 170 * 171 * dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 172 * = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 173 * 174 * @param[in] Q Number of quadrature points 175 * @param[in] i Current quadrature point 176 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 177 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 178 */ 179 CEED_QFUNCTION_HELPER void InvertBoundaryMappingJacobian_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar dXdx[2][3]) { 180 CeedScalar dxdX[3][2]; 181 GradUnpack2(Q, i, 3, (CeedScalar *)dxdX_q, dxdX); 182 183 // dxdX_k,j * dxdX_j,k 184 CeedScalar dxdXTdxdX[2][2] = {{0.}}; 185 for (CeedInt j = 0; j < 2; j++) { 186 for (CeedInt k = 0; k < 2; k++) { 187 for (CeedInt l = 0; l < 3; l++) dxdXTdxdX[j][k] += dxdX[l][j] * dxdX[l][k]; 188 } 189 } 190 191 const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1] - dxdXTdxdX[1][0] * dxdXTdxdX[0][1]; 192 193 // Compute inverse of dxdXTdxdX 194 CeedScalar dxdXTdxdX_inv[2][2]; 195 dxdXTdxdX_inv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX; 196 dxdXTdxdX_inv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX; 197 dxdXTdxdX_inv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX; 198 dxdXTdxdX_inv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX; 199 200 // Compute dXdx from dxdXTdxdX^-1 and dxdX 201 for (CeedInt j = 0; j < 2; j++) { 202 for (CeedInt k = 0; k < 3; k++) { 203 dXdx[j][k] = 0; 204 for (CeedInt l = 0; l < 2; l++) dXdx[j][k] += dxdXTdxdX_inv[l][j] * dxdX[k][l]; 205 } 206 } 207 } 208