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