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 HONEE 6 #pragma once 7 8 #include <ceed/types.h> 9 #include "utils.h" 10 11 /** 12 * @brief Calculate dXdx from dxdX for 3D elements 13 * 14 * Reference (parent) coordinates: X 15 * Physical (current) coordinates: x 16 * Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 17 * Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 18 * 19 * @param[in] Q Number of quadrature points 20 * @param[in] i Current quadrature point 21 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 22 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 23 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 24 */ 25 CEED_QFUNCTION_HELPER void InvertMappingJacobian_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar dXdx[3][3], 26 CeedScalar *detJ_ptr) { 27 CeedScalar dxdX[3][3]; 28 29 GradUnpack3D(Q, i, 3, (CeedScalar *)dxdX_q, dxdX); 30 MatInv3(dxdX, dXdx, detJ_ptr); 31 } 32 33 /** 34 * @brief Calculate dXdx from dxdX for 2D elements 35 * 36 * Reference (parent) coordinates: X 37 * Physical (current) coordinates: x 38 * Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 39 * Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 40 * 41 * @param[in] Q Number of quadrature points 42 * @param[in] i Current quadrature point 43 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 44 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 45 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 46 */ 47 CEED_QFUNCTION_HELPER void InvertMappingJacobian_2D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[2][CEED_Q_VLA], CeedScalar dXdx[2][2], 48 CeedScalar *detJ_ptr) { 49 CeedScalar dxdX[2][2]; 50 51 GradUnpack2D(Q, i, 2, (CeedScalar *)dxdX_q, dxdX); 52 MatInv2(dxdX, dXdx, detJ_ptr); 53 } 54 55 /** 56 * @brief Calculate face element's normal vector from dxdX 57 * 58 * Reference (parent) 2D coordinates: X 59 * Physical (current) 3D coordinates: x 60 * Change of coordinate matrix: 61 * dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 62 * Inverse change of coordinate matrix: 63 * dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 64 * 65 * (N1,N2,N3) is given by the cross product of the columns of dxdX_{i,j} 66 * 67 * detJb is the magnitude of (N1,N2,N3) 68 * 69 * Normal vector = (N1,N2,N3) / detJb 70 * 71 * @param[in] Q Number of quadrature points 72 * @param[in] i Current quadrature point 73 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 74 * @param[out] normal Inverse of mapping Jacobian at quadrature point i 75 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 76 */ 77 CEED_QFUNCTION_HELPER void NormalVectorFromdxdX_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar normal[3], 78 CeedScalar *detJ_ptr) { 79 CeedScalar dxdX[3][2]; 80 81 GradUnpack2D(Q, i, 3, (CeedScalar *)dxdX_q, dxdX); 82 // N1, N2, and N3 are given by the cross product of the columns of dxdX 83 normal[0] = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 84 normal[1] = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 85 normal[2] = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 86 87 const CeedScalar detJ = Norm3(normal); 88 ScaleN(normal, 1 / detJ, 3); 89 if (detJ_ptr) *detJ_ptr = detJ; 90 } 91 92 /** 93 * This QFunction sets up the geometric factor required for integration when reference coordinates are in 1D and the physical coordinates are in 2D 94 * 95 * Reference (parent) 1D coordinates: X 96 * Physical (current) 2D coordinates: x 97 * Change of coordinate vector: 98 * N1 = dx_1/dX 99 * N2 = dx_2/dX 100 * 101 * detJb is the magnitude of (N1,N2) 102 * 103 * We require the determinant of the Jacobian to properly compute integrals of the form: int( u v ) 104 * 105 * Normal vector is given by the cross product of (N1,N2)/detJ and ẑ 106 * 107 * @param[in] Q Number of quadrature points 108 * @param[in] i Current quadrature point 109 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 110 * @param[out] normal Inverse of mapping Jacobian at quadrature point i 111 * @param[out] detJ_ptr Determinate of the Jacobian, may be NULL is not desired 112 */ 113 CEED_QFUNCTION_HELPER void NormalVectorFromdxdX_2D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[CEED_Q_VLA], CeedScalar normal[2], 114 CeedScalar *detJ_ptr) { 115 normal[0] = dxdX_q[1][i]; 116 normal[1] = -dxdX_q[0][i]; 117 const CeedScalar detJb = Norm2(normal); 118 ScaleN(normal, 1 / detJb, 2); 119 if (detJ_ptr) *detJ_ptr = detJb; 120 } 121 122 /** 123 * @brief Calculate inverse of mapping Jacobian, (dxdX)^-1 124 * 125 * Reference (parent) 2D coordinates: X 126 * Physical (current) 3D coordinates: x 127 * Change of coordinate matrix: 128 * dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 129 * Inverse change of coordinate matrix: 130 * dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 131 * 132 * dXdx is calculated via Moore–Penrose inverse: 133 * 134 * dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 135 * = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 136 * 137 * @param[in] Q Number of quadrature points 138 * @param[in] i Current quadrature point 139 * @param[in] dxdX_q Mapping Jacobian (gradient of the coordinate space) 140 * @param[out] dXdx Inverse of mapping Jacobian at quadrature point i 141 */ 142 CEED_QFUNCTION_HELPER void InvertBoundaryMappingJacobian_3D(CeedInt Q, CeedInt i, const CeedScalar (*dxdX_q)[3][CEED_Q_VLA], CeedScalar dXdx[2][3]) { 143 CeedScalar dxdX[3][2]; 144 GradUnpack2D(Q, i, 3, (CeedScalar *)dxdX_q, dxdX); 145 146 // dxdX_k,j * dxdX_j,k 147 CeedScalar dxdXTdxdX[2][2] = {{0.}}; 148 for (CeedInt j = 0; j < 2; j++) { 149 for (CeedInt k = 0; k < 2; k++) { 150 for (CeedInt l = 0; l < 3; l++) dxdXTdxdX[j][k] += dxdX[l][j] * dxdX[l][k]; 151 } 152 } 153 154 const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1] - dxdXTdxdX[1][0] * dxdXTdxdX[0][1]; 155 156 // Compute inverse of dxdXTdxdX 157 CeedScalar dxdXTdxdX_inv[2][2]; 158 dxdXTdxdX_inv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX; 159 dxdXTdxdX_inv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX; 160 dxdXTdxdX_inv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX; 161 dxdXTdxdX_inv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX; 162 163 // Compute dXdx from dxdXTdxdX^-1 and dxdX 164 for (CeedInt j = 0; j < 2; j++) { 165 for (CeedInt k = 0; k < 3; k++) { 166 dXdx[j][k] = 0; 167 for (CeedInt l = 0; l < 2; l++) dXdx[j][k] += dxdXTdxdX_inv[l][j] * dxdX[k][l]; 168 } 169 } 170 } 171