Operator.hpp 4.83 KB
Newer Older
1
2
#pragma once

3
4
#include <limits>

5
#include "utility/FunctorVector.hpp"
6
#include "Mapper.hpp"
7
8
9

namespace Dec
{
10
11
  /**
    * \defgroup operators Operators
12
13
    * \brief Predefined terms representing (differential) operators, to assemble a linear system.
    *
14
15
16
    * @{
    **/

17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
    /*
    * Operators are used to simplify the assembling of linear systems. Therefore, standard
    * differential operators are provided that can be chained to describe a partial differential
    * equation in simple terms.
    *
    * Basic usage:
    * ```
    * DOFMatrix<double> mat;
    * OperatorXYZ<GridView> op(gridview);
    *
    * op.build(mat); // assemble operator into matrix.
    * ```
    *
    * The method `build()` creates a new inserter to insert coefficients of one specific
    * operator. When combining several operators, use the `assemble()` method. This expects
    * an inserter instead of a matrix.
    * ```
    * op.init(mat);
    * { // local scope
    *   auto ins = mat.inserter();
    *   op.assemble(ins);
    * }
    * op.finish(mat);
    **/


Praetorius, Simon's avatar
Praetorius, Simon committed
43
44
  class OperatorAccessor
  {
45
    template <class Model, class GridView, int N, int M> friend class Operator;
Praetorius, Simon's avatar
Praetorius, Simon committed
46
47
48
49

    template <class Derived>
    static auto nzrows(Derived const& derived)
    {
50
      return derived.nzrows_impl();
Praetorius, Simon's avatar
Praetorius, Simon committed
51
52
    }

53
54
    template <class Derived, class Inserter>
    static void assemble(Derived const& derived, Inserter& A, float_type factor)
Praetorius, Simon's avatar
Praetorius, Simon committed
55
    {
56
      derived.assemble_impl(A, factor);
Praetorius, Simon's avatar
Praetorius, Simon committed
57
58
    }

59
60
    template <class Derived, class Inserter, class Entity>
    static void assembleRow(Derived const& derived, Inserter& A, Entity const& e, float_type factor)
Praetorius, Simon's avatar
Praetorius, Simon committed
61
    {
62
      derived.assembleRow_impl(A, e, factor);
Praetorius, Simon's avatar
Praetorius, Simon committed
63
64
    }
  };
65

Praetorius, Simon's avatar
Praetorius, Simon committed
66

67
68
  /// \brief Base class for all operators, using CRT-Pattern.
  /**
69
   * Usage: Specify class Derived<GridView> : Operator<Derived, GridView, N, M>.
70
71
72
   * The parameter N and M specify the grid dimension of the
   * rows/ column entities.
   **/
73
  template <class Model, class GridView, int N, int M>
74
75
  class Operator
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
76
77
    friend class OperatorAccessor;

78
79
80
    template <class,class,class,int,int> friend class ChainedOperator;
    template <class,class,int,int> friend class ScaleOperator;

81
82
    static_assert( 0 <= N && N <= GridView::dimension, "N out of range [0,dim]" );
    static_assert( 0 <= M && M <= GridView::dimension, "M out of range [0,dim]" );
83

84
  public:
85
    static constexpr int dimension = GridView::dimension;
86

87
    /// Constructor, stores a reference to the `grid`.
88
89
    Operator(GridView const& gv)
      : gv_(gv)
90
91
    {}

92

93
94
  public: // Methods:

95
    /// Returns reference to the grid
96
    GridView const& gridView() const
97
    {
98
      return gv_;
99
100
101
    }

    /// Returns the number of rows of the matrix
102
103
    std::size_t rows() const
    {
104
      return gv_.size(dimension - N);
105
106
    }

107
    /// Returns the number of columns of the matrix
108
109
    std::size_t cols() const
    {
110
      return gv_.size(dimension - M);
111
112
    }

113
    /// Returns the number of nonzeros in the major dimension
Praetorius, Simon's avatar
Praetorius, Simon committed
114
115
116
117
118
    auto nzrows() const
    {
      return OperatorAccessor::nzrows(derived());
    }

119
120
121
    template <class Matrix>
    void init(Matrix& A) const
    {
122
      A.resize(rows(), cols());
123
124
125
    }

    template <class Matrix>
126
    void finish(Matrix& /*A*/) const
127
    {}
128

129
130
    template <class Inserter>
    void assemble(Inserter& ins, float_type factor = 1) const
131
    {
132
      OperatorAccessor::assemble(derived(), ins, factor);
133
134
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
135
136
137
138
    template <class Matrix>
    void build(Matrix& A, float_type factor = 1) const
    {
      init(A);
139
      {
140
        auto ins = A.inserter();
141
142
        assemble(ins, factor);
      }
Praetorius, Simon's avatar
Praetorius, Simon committed
143
144
145
      finish(A);
    }

146

Praetorius, Simon's avatar
Praetorius, Simon committed
147
  protected:
Praetorius, Simon's avatar
Praetorius, Simon committed
148

149
150
    template <class Inserter, class Entity>
    void assembleRow(Inserter& A, Entity const& e, float_type factor) const
151
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
152
      OperatorAccessor::assembleRow(derived(), A, e, factor);
153
154
155
    }


Praetorius, Simon's avatar
Praetorius, Simon committed
156
157
158
    template <class T>
    T regularize(T const& x) const
    {
159
      return max(10*std::numeric_limits<T>::epsilon(), x);
Praetorius, Simon's avatar
Praetorius, Simon committed
160
161
162
    }


Praetorius, Simon's avatar
Praetorius, Simon committed
163
  protected: // Implementation details. Should be overridden in subclasses
164
165

    // default implementation: 1 columns per row
166
    auto nzrows_impl() const
167
    {
168
      return make_vector([](std::size_t) { return 20; });
169
170
171
    }

    // default implementation: traverse entities of dim N and assemble row
172
173
    template <class Inserter>
    void assemble_impl(Inserter& A, float_type factor) const
174
    {
175
      for (auto const& e : entities(gv_, Dim_<N>))
Praetorius, Simon's avatar
Praetorius, Simon committed
176
        assembleRow(A, e, factor);
177
178
179
    }

    // default implementation: do nothing
180
181
    template <class Inserter, class Entity>
    void assembleRow_impl(Inserter& /*A*/, Entity const& /*e*/, float_type /*factor*/) const
182
183
184
185
186
    {
      /* do nothing */
    }


Praetorius, Simon's avatar
Praetorius, Simon committed
187
  public:
188

Praetorius, Simon's avatar
Praetorius, Simon committed
189
    Model& derived()
190
191
192
193
    {
      return static_cast<Model&>(*this);
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
194
    Model const& derived() const
195
196
197
    {
      return static_cast<Model const&>(*this);
    }
198

199

200
201
  protected: // Member variables:

202
    GridView gv_;
203

204
205
  };

206
207
  /** @} */

208
} // end namespace Dec