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 #include <ceed/types.h>
9
setup(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)10 CEED_QFUNCTION(setup)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
11 const CeedScalar *w = in[0], (*J)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[1];
12 CeedScalar(*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
13
14 // Quadrature point loop
15 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
16 // Qdata stored in Voigt convention
17 // J: 0 2 q_data: 0 2 adj(J): J11 -J01
18 // 1 3 2 1 -J10 J00
19 const CeedScalar J00 = J[0][0][i];
20 const CeedScalar J10 = J[0][1][i];
21 const CeedScalar J01 = J[1][0][i];
22 const CeedScalar J11 = J[1][1][i];
23 const CeedScalar qw = w[i] / (J00 * J11 - J10 * J01);
24 q_data[0][i] = qw * (J01 * J01 + J11 * J11);
25 q_data[1][i] = qw * (J00 * J00 + J10 * J10);
26 q_data[2][i] = -qw * (J00 * J01 + J10 * J11);
27 } // End of Quadrature Point Loop
28
29 return 0;
30 }
31
diff(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)32 CEED_QFUNCTION(diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
33 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*ug)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[1];
34 CeedScalar(*vg)[2][CEED_Q_VLA] = (CeedScalar(*)[2][CEED_Q_VLA])out[0];
35
36 const CeedInt dim = 2;
37 const CeedScalar num_comp = 2;
38 const CeedScalar scale[2][2] = {
39 {1.0, 2.0},
40 {3.0, 4.0},
41 };
42
43 // Quadrature point loop
44 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
45 // Read qdata (dXdxdXdxT symmetric matrix)
46 // Stored in Voigt convention
47 // 0 2
48 // 2 1
49 const CeedScalar dXdxdXdxT[2][2] = {
50 {q_data[0][i], q_data[2][i]},
51 {q_data[2][i], q_data[1][i]}
52 };
53
54 // Apply Poisson operator
55 // j = direction of vg
56 for (CeedInt j = 0; j < dim; j++) {
57 for (CeedInt k = 0; k < num_comp; k++) {
58 vg[j][k][i] = (ug[0][k][i] * dXdxdXdxT[0][j] * scale[0][j] + ug[1][k][i] * dXdxdXdxT[1][j] * scale[1][j]);
59 }
60 }
61 } // End of Quadrature Point Loop
62
63 return 0;
64 }
65