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 (2D) for HONEE
6 #include <ceed/types.h>
7 #include "setupgeo_helpers.h"
8 #include "utils.h"
9
10 // *****************************************************************************
11 // This QFunction sets up the geometric factors required for integration and coordinate transformations
12 //
13 // Reference (parent) coordinates: X
14 // Physical (current) coordinates: x
15 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation)
16 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j}
17 //
18 // All quadrature data is stored in 10 field vector of quadrature data.
19 //
20 // We require the determinant of the Jacobian to properly compute integrals of the form: int( v u )
21 //
22 // Determinant of Jacobian:
23 // detJ = J11*J22 - J21*J12
24 // Jij = Jacobian entry ij
25 //
26 // Stored: w detJ
27 // in q_data[0]
28 //
29 // We require the transpose of the inverse of the Jacobian to properly compute integrals of the form: int( gradv u )
30 //
31 // Inverse of Jacobian:
32 // dXdx_i,j = Aij / detJ
33 // Aij = Adjugate ij
34 //
35 // Stored: Aij / detJ
36 // in q_data[1:4] as
37 // (detJ^-1) * [A11 A12]
38 // [A21 A22]
39 // *****************************************************************************
Setup2d(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)40 CEED_QFUNCTION(Setup2d)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
41 const CeedScalar(*J)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[0];
42 const CeedScalar(*w) = in[1];
43 CeedScalar(*q_data) = out[0];
44
45 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
46 CeedScalar dXdx[2][2], detJ;
47 InvertMappingJacobian_2D(Q, i, J, dXdx, &detJ);
48 const CeedScalar wdetJ = w[i] * detJ;
49
50 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data);
51 StoredValuesPack(Q, i, 1, 4, (const CeedScalar *)dXdx, q_data);
52 }
53 return 0;
54 }
55
56 // *****************************************************************************
57 // This QFunction sets up the geometric factor required for integration when reference coordinates are in 1D and the physical coordinates are in 2D
58 //
59 // Reference (parent) 1D coordinates: X
60 // Physical (current) 2D coordinates: x
61 // Change of coordinate vector:
62 // J1 = dx_1/dX
63 // J2 = dx_2/dX
64 //
65 // detJb is the magnitude of (J1,J2)
66 //
67 // All quadrature data is stored in 3 field vector of quadrature data.
68 //
69 // We require the determinant of the Jacobian to properly compute integrals of the form: int( u v )
70 //
71 // Stored: w detJb
72 // in q_data_sur[0]
73 //
74 // Normal vector is given by the cross product of (J1,J2)/detJ and ẑ
75 //
76 // Stored: (J1,J2,0) x (0,0,1) / detJb
77 // in q_data_sur[1:2] as
78 // (detJb^-1) * [ J2 ]
79 // [-J1 ]
80 // *****************************************************************************
SetupBoundary2d(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)81 CEED_QFUNCTION(SetupBoundary2d)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
82 const CeedScalar(*J)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
83 const CeedScalar(*w) = in[1];
84 CeedScalar(*q_data_sur) = out[0];
85
86 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
87 CeedScalar normal[2], detJb;
88 NormalVectorFromdxdX_2D(Q, i, J, normal, &detJb);
89 const CeedScalar wdetJ = w[i] * detJb;
90
91 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur);
92 StoredValuesPack(Q, i, 1, 2, normal, q_data_sur);
93 }
94 return 0;
95 }
96
97 // *****************************************************************************
98 // This QFunction sets up the geometric factor required for integration when reference coordinates are 2D and the physical coordinates are in 3D
99 //
100 // In otherwords, when a 2D topology element is embedded in a 3D physical space.
101 //
102 // Reference (parent) 2D coordinates: X
103 // Physical (current) 3D coordinates: x
104 // Change of coordinate matrix:
105 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2]
106 // Inverse change of coordinate matrix:
107 // dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3]
108 //
109 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j}
110 //
111 // detJb is the magnitude of (J1,J2,J3)
112 //
113 // dXdx is calculated via Moore–Penrose inverse:
114 //
115 // dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX
116 // = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k
117 //
118 // All quadrature data is stored in 10 field vector of quadrature data.
119 //
120 // We require the determinant of the Jacobian to properly compute integrals of
121 // the form: int( u v )
122 //
123 // Stored: w detJb
124 // in q_data_sur[0]
125 //
126 // Normal vector = (J1,J2,J3) / detJb
127 //
128 // Stored: (J1,J2,J3) / detJb
129 //
130 // Stored: dXdx_{i,j}
131 // in q_data_sur[1:6] as
132 // [dXdx_11 dXdx_12 dXdx_13]
133 // [dXdx_21 dXdx_22 dXdx_23]
134 // *****************************************************************************
Setup2D_3Dcoords(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)135 CEED_QFUNCTION(Setup2D_3Dcoords)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
136 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0];
137 const CeedScalar(*w) = in[1];
138 CeedScalar(*q_data_sur) = out[0];
139
140 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
141 CeedScalar detJb, normal[3], dXdx[2][3];
142
143 NormalVectorFromdxdX_3D(Q, i, J, normal, &detJb);
144 InvertBoundaryMappingJacobian_3D(Q, i, J, dXdx);
145 const CeedScalar wdetJ = w[i] * detJb;
146
147 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur);
148 StoredValuesPack(Q, i, 1, 6, (const CeedScalar *)dXdx, q_data_sur);
149 }
150 return 0;
151 }
152
153 /**
154 @brief Compute geometric factors for integration, gradient transformations, and coordinate transformations on element faces.
155
156 Reference (parent) 1D coordinates are given by `X` and physical (current) 2D coordinates are given by `x`.
157 The change of coordinate matrix is given by`dxdX_{i,j} = dx_i/dX_j (indicial notation) [2 * 1]`.
158
159 @param[in] ctx QFunction context, unused
160 @param[in] Q Number of quadrature points
161 @param[in] in Input arrays
162 - 0 - Jacobian of cell coordinates
163 - 1 - Jacobian of face coordinates
164 - 2 - quadrature weights
165 @param[out] out Output array
166 - 0 - qdata, `w detNb`, `dXdx`, and `N`
167
168 @return An error code: 0 - success, otherwise - failure
169 **/
Setup2DBoundaryGradient(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)170 CEED_QFUNCTION(Setup2DBoundaryGradient)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
171 const CeedScalar(*J_cell)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[0];
172 const CeedScalar(*J_face)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
173 const CeedScalar(*w) = in[2];
174 CeedScalar(*q_data_sur) = out[0];
175
176 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
177 CeedScalar detJ_face, normal[2], dXdx[2][2];
178
179 NormalVectorFromdxdX_2D(Q, i, J_face, normal, &detJ_face);
180 const CeedScalar wdetJ = w[i] * detJ_face;
181 InvertMappingJacobian_2D(Q, i, J_cell, dXdx, NULL);
182
183 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur);
184 StoredValuesPack(Q, i, 1, 4, (CeedScalar *)dXdx, q_data_sur);
185 StoredValuesPack(Q, i, 5, 2, normal, q_data_sur);
186 }
187 return CEED_ERROR_SUCCESS;
188 }
189