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