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