ZeroOrderTest.hpp 2.26 KB
Newer Older
1
2
3
4
#pragma once

#include <type_traits>

5
#include <amdis/GridFunctionOperator.hpp>
6
#include <amdis/LocalBasisCache.hpp>
7
#include <amdis/common/ValueCategory.hpp>
8
9
10

namespace AMDiS
{
11
12
13
14
15
  /**
   * \addtogroup operators
   * @{
   **/

16
17
18
19
20
21
  namespace tag
  {
    struct test {};
  }


22
  /// zero-order vector-operator \f$ (c\, \psi) \f$
23
24
25
26
  template <class LocalContext, class GridFct>
  class GridFunctionOperator<tag::test, LocalContext, GridFct>
      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test, LocalContext, GridFct>,
                                        LocalContext, GridFct>
27
  {
28
    using Self = GridFunctionOperator;
29
    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
30

31
    static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
32
33

  public:
34
35
    GridFunctionOperator(tag::test, GridFct const& expr)
      : Super(expr, 0)
36
37
    {}

38
39
40
41
    template <class Context, class Node, class ElementVector>
    void getElementVector(Context const& context,
                          Node const& node,
                          ElementVector& elementVector)
42
    {
43
      static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only");
44
45
46

      auto const& localFE = node.finiteElement();

47
      using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
48
      LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
49
50

      auto const& quad = this->getQuadratureRule(context.type(), node);
51
52
      auto const& shapeValuesCache = cache.values(context,quad);
      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
53
        // Position of the current quadrature point in the reference element
54
        decltype(auto) local = context.local(quad[iq].position());
55
56

        // The multiplicative factor in the integral transformation formula
57
58
        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position())
                                                      * quad[iq].weight();
59

60
        auto const& shapeValues = shapeValuesCache[iq];
61
        for (std::size_t i = 0; i < localFE.size(); ++i) {
62
63
          const auto local_i = node.localIndex(i);
          elementVector[local_i] += factor * shapeValues[i];
64
65
66
67
68
        }
      }
    }
  };

69
70
  /** @} **/

71
} // end namespace AMDiS