xref: /libCEED/examples/mfem/bp1.h (revision c9c2c07970382857cc7b4a28d359710237b91a3e)
13d8e8822SJeremy L Thompson // Copyright (c) 2017-2022, 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 
11*c9c2c079SJeremy L Thompson #include <ceed.h>
12*c9c2c079SJeremy L Thompson 
134d537eeaSYohann /// A structure used to pass additional data to f_build_mass
144d537eeaSYohann struct BuildContext { CeedInt dim, space_dim; };
154d537eeaSYohann 
164d537eeaSYohann /// libCEED Q-function for building quadrature data for a mass operator
174d537eeaSYohann CEED_QFUNCTION(f_build_mass)(void *ctx, const CeedInt Q,
184d537eeaSYohann                              const CeedScalar *const *in, CeedScalar *const *out) {
194d537eeaSYohann   // in[0] is Jacobians with shape [dim, nc=dim, Q]
204d537eeaSYohann   // in[1] is quadrature weights, size (Q)
214d537eeaSYohann   BuildContext *bc = (BuildContext *)ctx;
22a2fa7910SValeria Barra   const CeedScalar *J = in[0], *w = in[1];
23a2fa7910SValeria Barra   CeedScalar *qdata = out[0];
24ee07ded2SValeria Barra 
254d537eeaSYohann   switch (bc->dim + 10*bc->space_dim) {
264d537eeaSYohann   case 11:
27ee07ded2SValeria Barra     // Quadrature Point Loop
28ee07ded2SValeria Barra     CeedPragmaSIMD
294d537eeaSYohann     for (CeedInt i=0; i<Q; i++) {
30a2fa7910SValeria Barra       qdata[i] = J[i] * w[i];
314d537eeaSYohann     }
324d537eeaSYohann     break;
334d537eeaSYohann   case 22:
34ee07ded2SValeria Barra     // Quadrature Point Loop
35ee07ded2SValeria Barra     CeedPragmaSIMD
364d537eeaSYohann     for (CeedInt i=0; i<Q; i++) {
374d537eeaSYohann       // 0 2
384d537eeaSYohann       // 1 3
39a2fa7910SValeria Barra       qdata[i] = (J[i+Q*0]*J[i+Q*3] - J[i+Q*1]*J[i+Q*2]) * w[i];
404d537eeaSYohann     }
414d537eeaSYohann     break;
424d537eeaSYohann   case 33:
43ee07ded2SValeria Barra     // Quadrature Point Loop
44ee07ded2SValeria Barra     CeedPragmaSIMD
454d537eeaSYohann     for (CeedInt i=0; i<Q; i++) {
464d537eeaSYohann       // 0 3 6
474d537eeaSYohann       // 1 4 7
484d537eeaSYohann       // 2 5 8
49a2fa7910SValeria Barra       qdata[i] = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) -
504d537eeaSYohann                   J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) +
51a2fa7910SValeria Barra                   J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6])) * w[i];
524d537eeaSYohann     }
534d537eeaSYohann     break;
544d537eeaSYohann   }
554d537eeaSYohann   return 0;
564d537eeaSYohann }
574d537eeaSYohann 
584d537eeaSYohann /// libCEED Q-function for applying a mass operator
594d537eeaSYohann CEED_QFUNCTION(f_apply_mass)(void *ctx, const CeedInt Q,
604d537eeaSYohann                              const CeedScalar *const *in, CeedScalar *const *out) {
61a2fa7910SValeria Barra   const CeedScalar *u = in[0], *qdata = in[1];
624d537eeaSYohann   CeedScalar *v = out[0];
63ee07ded2SValeria Barra 
64ee07ded2SValeria Barra   // Quadrature Point Loop
65ee07ded2SValeria Barra   CeedPragmaSIMD
664d537eeaSYohann   for (CeedInt i=0; i<Q; i++) {
67a2fa7910SValeria Barra     v[i] = qdata[i] * u[i];
684d537eeaSYohann   }
694d537eeaSYohann   return 0;
704d537eeaSYohann }
71a57ca787Svaleriabarra 
72a57ca787Svaleriabarra #endif // bp1_h
73