xref: /libCEED/examples/ceed/ex1-volume.h (revision 7b6fac6b22b60e6240f91980eb9d87fbfb688f61)
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 
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(build_mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
17   struct BuildContext *build_data = (struct BuildContext *)ctx;
18 
19   // in[0] is Jacobians with shape [dim, dim, Q]
20   // in[1] is quadrature weights with shape [1, Q]
21   const CeedScalar *w      = in[1];
22   CeedScalar       *q_data = out[0];
23 
24   switch (build_data->dim + 10 * build_data->space_dim) {
25     case 11: {
26       const CeedScalar(*J)[1][CEED_Q_VLA] = (const CeedScalar(*)[1][CEED_Q_VLA])in[0];
27 
28       // Quadrature Point Loop
29       CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { q_data[i] = J[0][0][i] * w[i]; }  // End of Quadrature Point Loop
30     } break;
31     case 22: {
32       const CeedScalar(*J)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[0];
33 
34       // Quadrature Point Loop
35       CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
36         q_data[i] = (J[0][0][i] * J[1][1][i] - J[0][1][i] * J[1][0][i]) * w[i];
37       }  // End of Quadrature Point Loop
38     } break;
39     case 33: {
40       const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0];
41 
42       // Quadrature Point Loop
43       CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
44         q_data[i] =
45             (J[0][0][i] * (J[1][1][i] * J[2][2][i] - J[1][2][i] * J[2][1][i]) - J[0][1][i] * (J[1][0][i] * J[2][2][i] - J[1][2][i] * J[2][0][i]) +
46              J[0][2][i] * (J[1][0][i] * J[2][1][i] - J[1][1][i] * J[2][0][i])) *
47             w[i];
48       }  // End of Quadrature Point Loop
49     } break;
50   }
51   return CEED_ERROR_SUCCESS;
52 }
53 
54 /// libCEED Q-function for applying a mass operator
55 CEED_QFUNCTION(apply_mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
56   // in[0], out[0] are solution variables with shape [1, Q]
57   // in[1] is quadrature data with shape [1, Q]
58   const CeedScalar *u = in[0], *q_data = in[1];
59   CeedScalar       *v = out[0];
60 
61   // Quadrature Point Loop
62   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { v[i] = q_data[i] * u[i]; }  // End of Quadrature Point Loop
63   return CEED_ERROR_SUCCESS;
64 }
65