OperatorTermBase.hpp 3.72 KB
Newer Older
1
2
#pragma once

3
4
5
6
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

7
8
#include <vector>
#include <type_traits>
9

10
#include <dune/geometry/quadraturerules.hh>
11

12
13
#include <dune/amdis/OperatorEvaluation.hpp>
#include <dune/amdis/common/ValueCategory.hpp>
14

15
namespace AMDiS
16
{
17

18
  /// Abstract base class of all operator terms
19
20
21
22
23
24
  template <class MeshView>
  class OperatorTerm
  {
  protected:
    using Codim0  = typename MeshView::template Codim<0>;
    using Element = typename Codim0::Entity;
25

26
27
    static constexpr int dim = MeshView::dimension;
    static constexpr int dow = MeshView::dimensionworld;
28

29
30
    using QuadratureRule = Dune::QuadratureRule<double, dim>;
    using PointList = std::vector<Dune::QuadraturePoint<double, dim>>;
31

32
  public:
33
    virtual void init(Element const& element,
34
                      PointList const& points) = 0;
35
36
37

    virtual double evalZot(size_t iq,
        Dune::FieldVector<double,1> const& test,
38
        Dune::FieldVector<double,1> const trial = 1.0) const = 0;
39
40
41

    virtual double evalFot1(size_t iq,
        Dune::FieldVector<double,1> const& test,
42
        Dune::FieldVector<double,dow> const& grad_trial) const = 0;
43
44

    virtual double evalFot2(size_t iq,
45
        Dune::FieldVector<double,dow> const& grad_test,
46
        Dune::FieldVector<double,1> const trial = 1.0) const = 0;
47
48

    virtual double evalSot(size_t iq,
49
50
        Dune::FieldVector<double,dow> const& grad_test,
        Dune::FieldVector<double,dow> const& grad_trial) const = 0;
51
52

    virtual int getDegree(Dune::GeometryType const& t) const = 0;
53
  };
54
55
56



57
58
  /// Base class for all operators based on expressions
  template <class MeshView, class Term, class Traits = tag::none>
59
  class GenericOperatorTerm
60
61
      : public OperatorTerm<MeshView>
      , public OperatorEvaluation
62
  {
63
64
65
    using Super   = OperatorTerm<MeshView>;
    using Element = typename Super::Element;
    using PointList = typename Super::PointList;
66

67
68
    static constexpr int dim = MeshView::dimension;
    static constexpr int dow = MeshView::dimensionworld;
69
70
71
72
73
74

  public:
    GenericOperatorTerm(Term const& term)
      : term(term)
    {}

75
    virtual void init(Element const& element,
76
77
78
                      PointList const& points) override
    {
      term.init(element, points);
79

80
81
82
83
84
      // cache term evaluation
      values.resize(points.size());
      for (size_t iq = 0; iq < points.size(); ++iq)
        values[iq] = term[iq];
    }
85
86
87

    virtual double evalZot(size_t iq,
        Dune::FieldVector<double,1> const& test,
88
89
90
91
        Dune::FieldVector<double,1> const trial = 1.0) const override
    {
      return this->evalZotImpl(_cat{}, _traits{}, values[iq], test, trial);
    }
92
93
94

    virtual double evalFot1(size_t iq,
        Dune::FieldVector<double,1> const& test,
95
        Dune::FieldVector<double,dow> const& grad_trial) const override
96
97
98
    {
      return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_trial, test);
    }
99
100

    virtual double evalFot2(size_t iq,
101
        Dune::FieldVector<double,dow> const& grad_test,
102
103
104
105
        Dune::FieldVector<double,1> const trial = 1.0) const override
    {
      return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_test, trial);
    }
106
107

    virtual double evalSot(size_t iq,
108
109
        Dune::FieldVector<double,dow> const& grad_test,
        Dune::FieldVector<double,dow> const& grad_trial) const override
110
111
112
    {
      return this->evalSotImpl(_cat{}, _traits{}, values[iq], grad_test, grad_trial);
    }
113
114

    virtual int getDegree(Dune::GeometryType const& t) const override
115
    {
116
      return term.getDegree(t);
117
    }
118

119
120
  private:
    Term term;
121

122
123
124
    using value_type = std::decay_t< decltype( std::declval<Term>()[std::declval<size_t>()] ) >;
    using _cat    = ValueCategory_t<value_type>;
    using _traits = Traits;
125

126
127
    std::vector<value_type> values;
  };
128

129
} // end namespace AMDiS