OperatorTermBase.hpp 3.6 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
    static constexpr int dim = Element::dimension;
27

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

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

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

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

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

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

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



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

66
67
68
69
70
71
72
    static constexpr int dim = Element::dimension;

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

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

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

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

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

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

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

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

117
118
  private:
    Term term;
119

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

124
125
    std::vector<value_type> values;
  };
126

127
} // end namespace AMDiS