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 #include <ceed/types.h>
7
8 #include "setupgeo_helpers.h"
9 #include "utils.h"
10
11 // *****************************************************************************
12 // This QFunction sets up the geometric factors required for integration and coordinate transformations
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 // All quadrature data is stored in 10 field vector of quadrature data.
20 //
21 // We require the determinant of the Jacobian to properly compute integrals of the form: int( v u )
22 //
23 // Determinant of Jacobian:
24 // detJ = J11*A11 + J21*A12 + J31*A13
25 // Jij = Jacobian entry ij
26 // Aij = Adjugate ij
27 //
28 // Stored: w detJ
29 // in q_data[0]
30 //
31 // We require the transpose of the inverse of the Jacobian to properly compute integrals of the form: int( gradv u )
32 //
33 // Inverse of Jacobian:
34 // dXdx_i,j = Aij / detJ
35 //
36 // Stored: Aij / detJ
37 // in q_data[1:9] as
38 // (detJ^-1) * [A11 A12 A13]
39 // [A21 A22 A23]
40 // [A31 A32 A33]
41 // *****************************************************************************
Setup(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)42 CEED_QFUNCTION(Setup)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
43 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0];
44 const CeedScalar(*w) = in[1];
45 CeedScalar(*q_data) = out[0];
46
47 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
48 CeedScalar detJ, dXdx[3][3];
49 InvertMappingJacobian_3D(Q, i, J, dXdx, &detJ);
50 const CeedScalar wdetJ = w[i] * detJ;
51
52 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data);
53 StoredValuesPack(Q, i, 1, 9, (const CeedScalar *)dXdx, q_data);
54 }
55 return 0;
56 }
57
58 // *****************************************************************************
59 // This QFunction sets up the geometric factor required for integration when reference coordinates are in 2D and the physical coordinates are in 3D
60 //
61 // Reference (parent) 2D coordinates: X
62 // Physical (current) 3D coordinates: x
63 // Change of coordinate matrix:
64 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2]
65 // Inverse change of coordinate matrix:
66 // dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3]
67 //
68 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j}
69 //
70 // detJb is the magnitude of (J1,J2,J3)
71 //
72 // dXdx is calculated via Moore–Penrose inverse:
73 //
74 // dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX
75 // = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k
76 //
77 // All quadrature data is stored in 10 field vector of quadrature data.
78 //
79 // We require the determinant of the Jacobian to properly compute integrals of
80 // the form: int( u v )
81 //
82 // Stored: w detJb
83 // in q_data_sur[0]
84 //
85 // Normal vector = (J1,J2,J3) / detJb
86 //
87 // - TODO Could possibly remove normal vector, as it could be calculated in the Qfunction from dXdx
88 // See https://github.com/CEED/libCEED/pull/868#discussion_r871979484
89 // Stored: (J1,J2,J3) / detJb
90 // in q_data_sur[1:3] as
91 // (detJb^-1) * [ J1 ]
92 // [ J2 ]
93 // [ J3 ]
94 //
95 // Stored: dXdx_{i,j}
96 // in q_data_sur[4:9] as
97 // [dXdx_11 dXdx_12 dXdx_13]
98 // [dXdx_21 dXdx_22 dXdx_23]
99 // *****************************************************************************
SetupBoundary(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)100 CEED_QFUNCTION(SetupBoundary)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
101 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0];
102 const CeedScalar(*w) = in[1];
103 CeedScalar(*q_data_sur) = out[0];
104
105 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
106 CeedScalar detJb, normal[3], dXdx[2][3];
107
108 NormalVectorFromdxdX_3D(Q, i, J, normal, &detJb);
109 InvertBoundaryMappingJacobian_3D(Q, i, J, dXdx);
110 const CeedScalar wdetJ = w[i] * detJb;
111
112 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur);
113 StoredValuesPack(Q, i, 1, 3, normal, q_data_sur);
114 StoredValuesPack(Q, i, 4, 6, (const CeedScalar *)dXdx, q_data_sur);
115 }
116 return 0;
117 }
118
119 /**
120 @brief Compute geometric factors for integration, gradient transformations, and coordinate transformations on element faces.
121
122 Reference (parent) 2D coordinates are given by `X` and physical (current) 3D coordinates are given by `x`.
123 The change of coordinate matrix is given by`dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2]`.
124
125 `(N_1, N_2, N_3)` is given by the cross product of the columns of `dxdX_{i,j}`.
126
127 `detNb` is the magnitude of `(N_1, N_2, N_3)`.
128
129 @param[in] ctx QFunction context, unused
130 @param[in] Q Number of quadrature points
131 @param[in] in Input arrays
132 - 0 - Jacobian of cell coordinates
133 - 1 - Jacobian of face coordinates
134 - 2 - quadrature weights
135 @param[out] out Output array
136 - 0 - qdata, `w detNb`, `dXdx`, and `N`
137
138 @return An error code: 0 - success, otherwise - failure
139 **/
SetupBoundaryGradient(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)140 CEED_QFUNCTION(SetupBoundaryGradient)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
141 const CeedScalar(*J_cell)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0];
142 const CeedScalar(*J_face)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1];
143 const CeedScalar(*w) = in[2];
144 CeedScalar(*q_data_sur) = out[0];
145
146 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
147 CeedScalar detJ_face, normal[3], dXdx[3][3];
148
149 NormalVectorFromdxdX_3D(Q, i, J_face, normal, &detJ_face);
150 const CeedScalar wdetJ = w[i] * detJ_face;
151 InvertMappingJacobian_3D(Q, i, J_cell, dXdx, NULL);
152
153 StoredValuesPack(Q, i, 0, 1, &wdetJ, q_data_sur);
154 StoredValuesPack(Q, i, 1, 9, (CeedScalar *)dXdx, q_data_sur);
155 StoredValuesPack(Q, i, 10, 3, normal, q_data_sur);
156 }
157 return CEED_ERROR_SUCCESS;
158 }
159