xref: /libCEED/examples/petsc/qfunctions/bps/bp1.h (revision 3d8e882215d238700cdceb37404f76ca7fa24eaa)
1 // Copyright (c) 2017-2022, 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 /// @file
9 /// libCEED QFunctions for mass operator example using PETSc
10 
11 #ifndef bp1_h
12 #define bp1_h
13 
14 #include <math.h>
15 
16 // -----------------------------------------------------------------------------
17 // This QFunction sets up the geometric factors required to apply the
18 //   mass operator
19 //
20 // The quadrature data is stored in the array q_data.
21 //
22 // We require the determinant of the Jacobian to properly compute integrals of
23 //   the form: int( u v )
24 //
25 // Qdata: det_J * w
26 //
27 // -----------------------------------------------------------------------------
28 CEED_QFUNCTION(SetupMassGeo)(void *ctx, const CeedInt Q,
29                              const CeedScalar *const *in,
30                              CeedScalar *const *out) {
31   const CeedScalar *J = in[1], *w = in[2]; // Note: *X = in[0]
32   CeedScalar *q_data = out[0];
33 
34   // Quadrature Point Loop
35   CeedPragmaSIMD
36   for (CeedInt i=0; i<Q; i++) {
37     const CeedScalar det_J = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) -
38                              J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) +
39                              J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6]));
40     q_data[i] = det_J * w[i];
41   } // End of Quadrature Point Loop
42   return 0;
43 }
44 
45 // -----------------------------------------------------------------------------
46 // This QFunction sets up the rhs and true solution for the problem
47 // -----------------------------------------------------------------------------
48 CEED_QFUNCTION(SetupMassRhs)(void *ctx, const CeedInt Q,
49                              const CeedScalar *const *in,
50                              CeedScalar *const *out) {
51   const CeedScalar *x = in[0], *w = in[1];
52   CeedScalar *true_soln = out[0], *rhs = out[1];
53 
54   // Quadrature Point Loop
55   CeedPragmaSIMD
56   for (CeedInt i=0; i<Q; i++) {
57     true_soln[i] = sqrt(x[i]*x[i] + x[i+Q]*x[i+Q] + x[i+2*Q]*x[i+2*Q]);
58     rhs[i] = w[i] * true_soln[i];
59   } // End of Quadrature Point Loop
60   return 0;
61 }
62 
63 // -----------------------------------------------------------------------------
64 // This QFunction applies the mass operator for a scalar field.
65 //
66 // Inputs:
67 //   u     - Input vector at quadrature points
68 //   q_data - Geometric factors
69 //
70 // Output:
71 //   v     - Output vector (test functions) at quadrature points
72 //
73 // -----------------------------------------------------------------------------
74 CEED_QFUNCTION(Mass)(void *ctx, const CeedInt Q,
75                      const CeedScalar *const *in, CeedScalar *const *out) {
76   const CeedScalar *u = in[0], *q_data = in[1];
77   CeedScalar *v = out[0];
78 
79   // Quadrature Point Loop
80   CeedPragmaSIMD
81   for (CeedInt i=0; i<Q; i++)
82     v[i] = q_data[i] * u[i];
83 
84   return 0;
85 }
86 // -----------------------------------------------------------------------------
87 
88 #endif // bp1_h
89