xref: /libCEED/examples/mfem/bp1.h (revision 3d8e882215d238700cdceb37404f76ca7fa24eaa)
1*3d8e8822SJeremy L Thompson // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
2*3d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
34d537eeaSYohann //
4*3d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause
54d537eeaSYohann //
6*3d8e8822SJeremy L Thompson // This file is part of CEED:  http://github.com/ceed
74d537eeaSYohann 
8a57ca787Svaleriabarra #ifndef bp1_h
9a57ca787Svaleriabarra #define bp1_h
10a57ca787Svaleriabarra 
114d537eeaSYohann /// A structure used to pass additional data to f_build_mass
124d537eeaSYohann struct BuildContext { CeedInt dim, space_dim; };
134d537eeaSYohann 
144d537eeaSYohann /// libCEED Q-function for building quadrature data for a mass operator
154d537eeaSYohann CEED_QFUNCTION(f_build_mass)(void *ctx, const CeedInt Q,
164d537eeaSYohann                              const CeedScalar *const *in, CeedScalar *const *out) {
174d537eeaSYohann   // in[0] is Jacobians with shape [dim, nc=dim, Q]
184d537eeaSYohann   // in[1] is quadrature weights, size (Q)
194d537eeaSYohann   BuildContext *bc = (BuildContext *)ctx;
20a2fa7910SValeria Barra   const CeedScalar *J = in[0], *w = in[1];
21a2fa7910SValeria Barra   CeedScalar *qdata = out[0];
22ee07ded2SValeria Barra 
234d537eeaSYohann   switch (bc->dim + 10*bc->space_dim) {
244d537eeaSYohann   case 11:
25ee07ded2SValeria Barra     // Quadrature Point Loop
26ee07ded2SValeria Barra     CeedPragmaSIMD
274d537eeaSYohann     for (CeedInt i=0; i<Q; i++) {
28a2fa7910SValeria Barra       qdata[i] = J[i] * w[i];
294d537eeaSYohann     }
304d537eeaSYohann     break;
314d537eeaSYohann   case 22:
32ee07ded2SValeria Barra     // Quadrature Point Loop
33ee07ded2SValeria Barra     CeedPragmaSIMD
344d537eeaSYohann     for (CeedInt i=0; i<Q; i++) {
354d537eeaSYohann       // 0 2
364d537eeaSYohann       // 1 3
37a2fa7910SValeria Barra       qdata[i] = (J[i+Q*0]*J[i+Q*3] - J[i+Q*1]*J[i+Q*2]) * w[i];
384d537eeaSYohann     }
394d537eeaSYohann     break;
404d537eeaSYohann   case 33:
41ee07ded2SValeria Barra     // Quadrature Point Loop
42ee07ded2SValeria Barra     CeedPragmaSIMD
434d537eeaSYohann     for (CeedInt i=0; i<Q; i++) {
444d537eeaSYohann       // 0 3 6
454d537eeaSYohann       // 1 4 7
464d537eeaSYohann       // 2 5 8
47a2fa7910SValeria Barra       qdata[i] = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) -
484d537eeaSYohann                   J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) +
49a2fa7910SValeria Barra                   J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6])) * w[i];
504d537eeaSYohann     }
514d537eeaSYohann     break;
524d537eeaSYohann   }
534d537eeaSYohann   return 0;
544d537eeaSYohann }
554d537eeaSYohann 
564d537eeaSYohann /// libCEED Q-function for applying a mass operator
574d537eeaSYohann CEED_QFUNCTION(f_apply_mass)(void *ctx, const CeedInt Q,
584d537eeaSYohann                              const CeedScalar *const *in, CeedScalar *const *out) {
59a2fa7910SValeria Barra   const CeedScalar *u = in[0], *qdata = in[1];
604d537eeaSYohann   CeedScalar *v = out[0];
61ee07ded2SValeria Barra 
62ee07ded2SValeria Barra   // Quadrature Point Loop
63ee07ded2SValeria Barra   CeedPragmaSIMD
644d537eeaSYohann   for (CeedInt i=0; i<Q; i++) {
65a2fa7910SValeria Barra     v[i] = qdata[i] * u[i];
664d537eeaSYohann   }
674d537eeaSYohann   return 0;
684d537eeaSYohann }
69a57ca787Svaleriabarra 
70a57ca787Svaleriabarra #endif // bp1_h
71