xref: /libCEED/examples/mfem/bp1.h (revision 5aed82e4fa97acf4ba24a7f10a35f5303a6798e0)
1*5aed82e4SJeremy L Thompson // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors.
23d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
34d537eeaSYohann //
43d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause
54d537eeaSYohann //
63d8e8822SJeremy L Thompson // This file is part of CEED:  http://github.com/ceed
74d537eeaSYohann 
8a57ca787Svaleriabarra #ifndef bp1_h
9a57ca787Svaleriabarra #define bp1_h
10a57ca787Svaleriabarra 
11c9c2c079SJeremy L Thompson #include <ceed.h>
12c9c2c079SJeremy L Thompson 
134d537eeaSYohann /// A structure used to pass additional data to f_build_mass
142b730f8bSJeremy L Thompson struct BuildContext {
152b730f8bSJeremy L Thompson   CeedInt dim, space_dim;
162b730f8bSJeremy L Thompson };
174d537eeaSYohann 
184d537eeaSYohann /// libCEED Q-function for building quadrature data for a mass operator
192b730f8bSJeremy L Thompson CEED_QFUNCTION(f_build_mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
204d537eeaSYohann   // in[0] is Jacobians with shape [dim, nc=dim, Q]
214d537eeaSYohann   // in[1] is quadrature weights, size (Q)
224d537eeaSYohann   BuildContext     *bc = (BuildContext *)ctx;
23a2fa7910SValeria Barra   const CeedScalar *J = in[0], *w = in[1];
24a2fa7910SValeria Barra   CeedScalar       *qdata = out[0];
25ee07ded2SValeria Barra 
264d537eeaSYohann   switch (bc->dim + 10 * bc->space_dim) {
274d537eeaSYohann     case 11:
28ee07ded2SValeria Barra       // Quadrature Point Loop
292b730f8bSJeremy L Thompson       CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { qdata[i] = J[i] * w[i]; }
304d537eeaSYohann       break;
314d537eeaSYohann     case 22:
32ee07ded2SValeria Barra       // Quadrature Point Loop
332b730f8bSJeremy L Thompson       CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
344d537eeaSYohann         // 0 2
354d537eeaSYohann         // 1 3
36a2fa7910SValeria Barra         qdata[i] = (J[i + Q * 0] * J[i + Q * 3] - J[i + Q * 1] * J[i + Q * 2]) * w[i];
374d537eeaSYohann       }
384d537eeaSYohann       break;
394d537eeaSYohann     case 33:
40ee07ded2SValeria Barra       // Quadrature Point Loop
412b730f8bSJeremy L Thompson       CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
424d537eeaSYohann         // 0 3 6
434d537eeaSYohann         // 1 4 7
444d537eeaSYohann         // 2 5 8
45a2fa7910SValeria Barra         qdata[i] = (J[i + Q * 0] * (J[i + Q * 4] * J[i + Q * 8] - J[i + Q * 5] * J[i + Q * 7]) -
464d537eeaSYohann                     J[i + Q * 1] * (J[i + Q * 3] * J[i + Q * 8] - J[i + Q * 5] * J[i + Q * 6]) +
472b730f8bSJeremy L Thompson                     J[i + Q * 2] * (J[i + Q * 3] * J[i + Q * 7] - J[i + Q * 4] * J[i + Q * 6])) *
482b730f8bSJeremy L Thompson                    w[i];
494d537eeaSYohann       }
504d537eeaSYohann       break;
514d537eeaSYohann   }
524d537eeaSYohann   return 0;
534d537eeaSYohann }
544d537eeaSYohann 
554d537eeaSYohann /// libCEED Q-function for applying a mass operator
562b730f8bSJeremy L Thompson CEED_QFUNCTION(f_apply_mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
57a2fa7910SValeria Barra   const CeedScalar *u = in[0], *qdata = in[1];
584d537eeaSYohann   CeedScalar       *v = out[0];
59ee07ded2SValeria Barra 
60ee07ded2SValeria Barra   // Quadrature Point Loop
612b730f8bSJeremy L Thompson   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { v[i] = qdata[i] * u[i]; }
624d537eeaSYohann   return 0;
634d537eeaSYohann }
64a57ca787Svaleriabarra 
65a57ca787Svaleriabarra #endif  // bp1_h
66