xref: /libCEED/tests/t522-operator.h (revision df9b1082f3fea590280bb79ace83362d3424df46)
1 // Copyright (c) 2017-2022, 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 CEED_QFUNCTION(setup)(void *ctx, const CeedInt Q,
9                       const CeedScalar *const *in,
10                       CeedScalar *const *out) {
11   const CeedScalar *qw = in[0], *J = in[1];
12   CeedScalar *qd = out[0];
13 
14   for (CeedInt i=0; i<Q; i++) {
15     // Qdata stored in Voigt convention
16     // J: 0 2   qd: 0 2   adj(J):  J22 -J12
17     //    1 3       2 1           -J21  J11
18     const CeedScalar J11 = J[i+Q*0];
19     const CeedScalar J21 = J[i+Q*1];
20     const CeedScalar J12 = J[i+Q*2];
21     const CeedScalar J22 = J[i+Q*3];
22     const CeedScalar w = qw[i] / (J11*J22 - J21*J12);
23     qd[i+Q*0] =   w * (J12*J12 + J22*J22);
24     qd[i+Q*1] =   w * (J11*J11 + J21*J21);
25     qd[i+Q*2] = - w * (J11*J12 + J21*J22);
26   } // End of Quadrature Point Loop
27   return 0;
28 }
29 
30 CEED_QFUNCTION(diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in,
31                      CeedScalar *const *out) {
32   const CeedScalar *qd = in[0], *ug = in[1];
33   CeedScalar *vg = out[0];
34 
35   for (CeedInt i=0; i<Q; i++) {
36     // Read spatial derivatives of u
37     const CeedScalar du[2]        =  {ug[i+Q*0],
38                                       ug[i+Q*1]
39                                      };
40 
41     // Read qdata (dXdxdXdxT symmetric matrix)
42     // Stored in Voigt convention
43     // 0 2
44     // 2 1
45     // *INDENT-OFF*
46     const CeedScalar dXdxdXdxT[2][2] = {{qd[i+0*Q],
47                                          qd[i+2*Q]},
48                                         {qd[i+2*Q],
49                                          qd[i+1*Q]}
50                                        };
51     // *INDENT-ON*
52     // j = direction of vg
53     for (int j=0; j<2; j++)
54       vg[i+j*Q] = (du[0] * dXdxdXdxT[0][j] +
55                    du[1] * dXdxdXdxT[1][j]);
56   }
57   return 0;
58 }
59