OperatorTermBase.hpp 3.74 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

  public:
71
    GenericOperatorTerm(Term const& term, Traits traits = {})
72
      : term(term)
73
      , traits(traits)
74
75
    {}

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

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

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

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

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

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

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

120
121
  private:
    Term term;
122
    Traits traits;
123

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

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

130
} // end namespace AMDiS