xref: /libCEED/examples/mfem/bp1.h (revision c172e0a6730c614967cdf24760cc19b8004cf29b)
1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
3 // All Rights reserved. See files LICENSE and NOTICE for details.
4 //
5 // This file is part of CEED, a collection of benchmarks, miniapps, software
6 // libraries and APIs for efficient high-order finite element and spectral
7 // element discretizations for exascale applications. For more information and
8 // source code availability see http://github.com/ceed.
9 //
10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11 // a collaborative effort of two U.S. Department of Energy organizations (Office
12 // of Science and the National Nuclear Security Administration) responsible for
13 // the planning and preparation of a capable exascale ecosystem, including
14 // software, applications, hardware, advanced system engineering and early
15 // testbed platforms, in support of the nation's exascale computing imperative.
16 
17 /// A structure used to pass additional data to f_build_mass
18 struct BuildContext { CeedInt dim, space_dim; };
19 
20 /// libCEED Q-function for building quadrature data for a mass operator
21 CEED_QFUNCTION(f_build_mass)(void *ctx, const CeedInt Q,
22                              const CeedScalar *const *in, CeedScalar *const *out) {
23   // in[0] is Jacobians with shape [dim, nc=dim, Q]
24   // in[1] is quadrature weights, size (Q)
25   BuildContext *bc = (BuildContext *)ctx;
26   const CeedScalar *J = in[0], *qw = in[1];
27   CeedScalar *rho = out[0];
28   switch (bc->dim + 10*bc->space_dim) {
29   case 11:
30     for (CeedInt i=0; i<Q; i++) {
31       rho[i] = J[i] * qw[i];
32     }
33     break;
34   case 22:
35     for (CeedInt i=0; i<Q; i++) {
36       // 0 2
37       // 1 3
38       rho[i] = (J[i+Q*0]*J[i+Q*3] - J[i+Q*1]*J[i+Q*2]) * qw[i];
39     }
40     break;
41   case 33:
42     for (CeedInt i=0; i<Q; i++) {
43       // 0 3 6
44       // 1 4 7
45       // 2 5 8
46       rho[i] = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) -
47                 J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) +
48                 J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6])) * qw[i];
49     }
50     break;
51   }
52   return 0;
53 }
54 
55 /// libCEED Q-function for applying a mass operator
56 CEED_QFUNCTION(f_apply_mass)(void *ctx, const CeedInt Q,
57                              const CeedScalar *const *in, CeedScalar *const *out) {
58   const CeedScalar *u = in[0], *w = in[1];
59   CeedScalar *v = out[0];
60   for (CeedInt i=0; i<Q; i++) {
61     v[i] = w[i] * u[i];
62   }
63   return 0;
64 }
65