ZeroOrderTestvec.hpp 2.07 KB
Newer Older
1
2
3
4
#pragma once

#include <type_traits>

5
6
#include <dune/amdis/GridFunctionOperator.hpp>
#include <dune/amdis/common/ValueCategory.hpp>
7
8
9
10
11
12
13
14
15
16

namespace AMDiS
{
  namespace tag
  {
    struct testvec {};
  }


  // zero-order vector-operator <b * Psi>
17
18
19
  template <class GridFct>
  class GridFunctionOperator<tag::testvec, GridFct>
      : public GridFunctionOperatorBase<GridFct>
20
  {
21
    using Super = GridFunctionOperatorBase<GridFct>;
22

23
    static_assert( Category::Vector<typename GridFct::Range>, "Expression must be of vector type." );
24
25

  public:
26
27
    GridFunctionOperator(tag::testvec, GridFct const& expr, int degree)
      : Super(expr, 0, degree)
28
29
30
31
32
33
34
35
36
    {}

    template <class Context, class QuadratureRule,
              class ElementVector, class Node>
    void calculateElementVector(Context const& context,
                                QuadratureRule const& quad,
                                ElementVector& elementVector,
                                Node const& node)
    {
37
38
      static_assert(Node::isPower,
        "Operator can be applied to Power-Nodes only.");
39
40
41
42
43
44
45
46
47
48

      static const std::size_t CHILDREN = Node::CHILDREN;

      auto const& localFE = node.child(0).finiteElement();

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
        decltype(auto) local = context.position(quad[iq].position());

        // The multiplicative factor in the integral transformation formula
49
50
        double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
        const auto exprValue = Super::coefficient(local);
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

        localFE.localBasis().evaluateFunction(local, shapeValues_);

        for (std::size_t i = 0; i < localFE.size(); ++i) {
          auto value = exprValue * (factor * shapeValues_[i]);
          for (std::size_t k = 0; k < CHILDREN; ++k) {
            const int local_ki = node.child(k).localIndex(i);
            elementVector[local_ki] += value[k];
          }
        }
      }
    }

  private:
    std::vector<Dune::FieldVector<double,1>> shapeValues_;
  };

} // end namespace AMDiS