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 */
InvertMappingJacobian_3D(CeedInt Q,CeedInt i,const CeedScalar (* dxdX_q)[3][CEED_Q_VLA],CeedScalar dXdx[3][3],CeedScalar * detJ_ptr)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 */
InvertMappingJacobian_2D(CeedInt Q,CeedInt i,const CeedScalar (* dxdX_q)[2][CEED_Q_VLA],CeedScalar dXdx[2][2],CeedScalar * detJ_ptr)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 */
NormalVectorFromdxdX_3D(CeedInt Q,CeedInt i,const CeedScalar (* dxdX_q)[3][CEED_Q_VLA],CeedScalar normal[3],CeedScalar * detJ_ptr)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 */
NormalVectorFromdxdX_2D(CeedInt Q,CeedInt i,const CeedScalar (* dxdX_q)[CEED_Q_VLA],CeedScalar normal[2],CeedScalar * detJ_ptr)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 */
InvertBoundaryMappingJacobian_3D(CeedInt Q,CeedInt i,const CeedScalar (* dxdX_q)[3][CEED_Q_VLA],CeedScalar dXdx[2][3])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