OperatorTermBase.hpp 5.79 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/amdis/QuadratureRules.hpp>
11

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

16
namespace AMDiS
17
{
18
  /// Abstract base class of all operator terms
19
  template <class LocalContext>
20
21
22
  class OperatorTerm
  {
  protected:
23
24
    static constexpr int dim = LocalContext::dimension;
    static constexpr int dow = LocalContext::Geometry::coorddimension;
25

26
    using QuadratureRule = QuadratureRuleFactory_t<LocalContext,double,dim>;
27

28
  public:
29
30
    // initialize operator-term on the current element
    virtual void init(LocalContext const& element, QuadratureRule const& points) = 0;
31

32
33
    // evaluate operator-term at quadrature point as zero-order term,
    // with no gradients involved
34
    virtual double evalZot(std::size_t iq,
35
        Dune::FieldVector<double,1> const& test,
36
        Dune::FieldVector<double,1> const trial = 1.0) const = 0;
37

38
39
    // evaluate operator-term at quadrature point as first-order term,
    // with gradients on the trial-function
40
    virtual double evalFot1(std::size_t iq,
41
        Dune::FieldVector<double,1> const& test,
42
        Dune::FieldVector<double,dow> const& grad_trial) const = 0;
43

44
45
    // evaluate operator-term at quadrature point as first-order term,
    // with gradients on the test-function
46
    virtual double evalFot2(std::size_t iq,
47
        Dune::FieldVector<double,dow> const& grad_test,
48
        Dune::FieldVector<double,1> const trial = 1.0) const = 0;
49

50
51
    // evaluate operator-term at quadrature point as second-order term,
    // with gradients on the trial and test-function
52
    virtual double evalSot(std::size_t iq,
53
54
        Dune::FieldVector<double,dow> const& grad_test,
        Dune::FieldVector<double,dow> const& grad_trial) const = 0;
55

56
57
    // return the polynomial degree, necessary to intergrate the operator-term
    // on an element accurately
58
    virtual int getDegree(Dune::GeometryType const& t) const = 0;
59
  };
60
61


62
63
64
65
66
67
68
69
70
71
72
73
  /// \brief Base class for all operator-terms based on expressions
  /**
   * Represents an operator-term with a coefficient the scales the classical operator, e.g.
   * - zero order terms `<expr * u, v>`
   * - first order terms (GRD_PHI) `<(vecexpr * grad(u)), v>`, `<expr * d_i(u), v>` (a)
   * - first order terms (GRD_PSI) `<vecexpr * u, grad(v)>`, `<expr * u, d_i(v)>` (a)
   * - second order terms `<expr * grad(u), grad(v)>`, `<(matexpr *grad(u)), grad(v)>`, `<expr * d_i(u), d_j(v)>` (b)
   **/
  // LocalContext ... The type of the element where to evaluate, either entity (codim=0), or intersection (codim=1)
  // Expression ..... The type of the expression to evaluate
  // Traits ......... Specialization for the operator evaluation [tag:none | VectorComponent (a) | MatrixComponent (b)]
  template <class LocalContext, class Expression, class Traits = tag::none>
74
  class GenericOperatorTerm
75
      : public OperatorTerm<LocalContext>
76
  {
77
    using Super = OperatorTerm<LocalContext>;
78
    using QuadratureRule = typename Super::QuadratureRule;
79

80
    static constexpr int dow = Super::dow;
81
82

  public:
83
84
85
    /// Constructor, stores a copy of the expression `expr` and of the traits `traits`.
    GenericOperatorTerm(Expression const& expr, Traits traits = {})
      : expr(expr)
86
      , traits(traits)
87
88
    {}

89
90
91
92
93
94
95
96
97
98
    /// \brief Initialize the expression on the current (inside-)element, in the quadrature points.
    /**
     * This initialization does not only bind the operator-term to the inside-element,
     * but also calculates values in all quadrature points. These pre-calculated values
     * can be accessed using the `operator[]` method of the expression.
     *
     * Note: For entities the expression is bound to the `element` directly, but
     * for intersections, the expression is bound to inside-entity of the element.
     **/
    virtual void init(LocalContext const& element, QuadratureRule const& points) override
99
    {
100
      expr.init(get_entity(element), points);
101
    }
102

103
    /// Calculates `expr[iq] * test * trial`
104
    virtual double evalZot(std::size_t iq,
105
        Dune::FieldVector<double,1> const& test,
106
107
        Dune::FieldVector<double,1> const trial = 1.0) const override
    {
108
      return Evaluate::zot(_cat{}, traits, expr[iq], test, trial);
109
    }
110

111
    /// Calculates `expr[iq] * test * grad(trial)`, maybe partial derivative only.
112
    virtual double evalFot1(std::size_t iq,
113
        Dune::FieldVector<double,1> const& test,
114
        Dune::FieldVector<double,dow> const& grad_trial) const override
115
    {
116
      return Evaluate::fot(_cat{}, traits, expr[iq], grad_trial, test);
117
    }
118

119
    /// Calculates `expr[iq] * grad(test) * trial`, maybe partial derivative only.
120
    virtual double evalFot2(std::size_t iq,
121
        Dune::FieldVector<double,dow> const& grad_test,
122
123
        Dune::FieldVector<double,1> const trial = 1.0) const override
    {
124
      return Evaluate::fot(_cat{}, traits, expr[iq], grad_test, trial);
125
    }
126

127
    /// Calculates `expr[iq] * grad(test) * grad(trial)`, maybe partial derivatives only.
128
    virtual double evalSot(std::size_t iq,
129
130
        Dune::FieldVector<double,dow> const& grad_test,
        Dune::FieldVector<double,dow> const& grad_trial) const override
131
    {
132
      return Evaluate::sot(_cat{}, traits, expr[iq], grad_test, grad_trial);
133
    }
134

135
    /// Returns the polynomial degree of the expression
136
    virtual int getDegree(Dune::GeometryType const& t) const override
137
    {
138
      return expr.getDegree(t);
139
    }
140

141
  private:
142
143
144
145
    /// A local expression, implementing the Concept \ref Concepts::LocalExpression.
    Expression expr;

    /// One of [tag::none, VectorComponent, MatrixComponent]
146
    Traits traits;
147

148
    using value_type = std::decay_t< decltype( std::declval<Expression>()[std::size_t(0)] ) >;
Praetorius, Simon's avatar
Praetorius, Simon committed
149
    using _cat = ValueCategory_t<value_type>;
150
  };
151

152
} // end namespace AMDiS