xref: /libCEED/examples/mfem/bp1.h (revision 3e551a327d6c97f9de071b988b42ffdb7bed19a7)
1 // Copyright (c) 2017-2024, 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 
10 /// A structure used to pass additional data to f_build_mass
11 struct BuildContext {
12   CeedInt dim, space_dim;
13 };
14 
15 /// libCEED Q-function for building quadrature data for a mass operator
16 CEED_QFUNCTION(f_build_mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
17   // in[0] is Jacobians with shape [dim, nc=dim, Q]
18   // in[1] is quadrature weights, size (Q)
19   BuildContext     *bc = (BuildContext *)ctx;
20   const CeedScalar *J = in[0], *w = in[1];
21   CeedScalar       *qdata = out[0];
22 
23   switch (bc->dim + 10 * bc->space_dim) {
24     case 11:
25       // Quadrature Point Loop
26       CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { qdata[i] = J[i] * w[i]; }
27       break;
28     case 22:
29       // Quadrature Point Loop
30       CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
31         // 0 2
32         // 1 3
33         qdata[i] = (J[i + Q * 0] * J[i + Q * 3] - J[i + Q * 1] * J[i + Q * 2]) * w[i];
34       }
35       break;
36     case 33:
37       // Quadrature Point Loop
38       CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
39         // 0 3 6
40         // 1 4 7
41         // 2 5 8
42         qdata[i] = (J[i + Q * 0] * (J[i + Q * 4] * J[i + Q * 8] - J[i + Q * 5] * J[i + Q * 7]) -
43                     J[i + Q * 1] * (J[i + Q * 3] * J[i + Q * 8] - J[i + Q * 5] * J[i + Q * 6]) +
44                     J[i + Q * 2] * (J[i + Q * 3] * J[i + Q * 7] - J[i + Q * 4] * J[i + Q * 6])) *
45                    w[i];
46       }
47       break;
48   }
49   return 0;
50 }
51 
52 /// libCEED Q-function for applying a mass operator
53 CEED_QFUNCTION(f_apply_mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
54   const CeedScalar *u = in[0], *qdata = in[1];
55   CeedScalar       *v = out[0];
56 
57   // Quadrature Point Loop
58   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { v[i] = qdata[i] * u[i]; }
59   return 0;
60 }
61