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