Operator.hpp 7.75 KB
Newer Older
1
2
3
#pragma once

#include <list>
4
#include <vector>
5

6
7
#include <dune/geometry/type.hh>

8
#include <dune/amdis/OperatorTermBase.hpp>
9
10
11
#include <dune/amdis/ZeroOrderAssembler.hpp>
#include <dune/amdis/FirstOrderAssembler.hpp>
#include <dune/amdis/SecondOrderAssembler.hpp>
12
13
14
15
#include <dune/amdis/terms/TermGenerator.hpp>

#include <dune/amdis/common/TypeDefs.hpp>
#include <dune/amdis/utility/GetDegree.hpp>
16
17
18

namespace AMDiS
{
19
20
21

  // the LocalContext is either an Codim=0-EntityType or an IntersectionType
  template <class MeshView, class LocalContext>
22
23
24
  class Operator
  {
    using Self = Operator;
25
    using OperatorTermType = OperatorTerm<LocalContext>;
26
27
28

    using ElementVector = Impl::ElementVector;
    using ElementMatrix = Impl::ElementMatrix;
29

30
    using IdType = typename MeshView::Grid::LocalIdSet::IdType;
31

32
33
34
35
    template <class, class> friend class ZeroOrderAssembler;
    template <class, class, FirstOrderType> friend class FirstOrderAssembler;
    template <class, class> friend class SecondOrderAssembler;

36
  public:
37
38
39
40
41
42
43
    /// \brief Add coefficients for zero-order operator < C(phi), psi >.
    /**
     * If `c` is a scalar expression, then C(phi) := c * phi,
     * otherwise if `c` is derived from \ref ZeroOrderOperatorTerm, it is evaluated
     * using its `evalZot` method at quadrature points and represents the local
     * evaluation of the scalar product `c = C(phi(lambda))*v(lambda)` with lambda a local coordinate.
     **/
44
    template <class Term>
45
    Self& addZOT(Term const& c)
46
    {
47
      return addZOTImpl(toTerm(c));
48
    }
49

50
    template <class... Args>
51
    static shared_ptr<Self> zot(Args&&... args)
52
    {
53
      auto op = make_shared<Self>();
54
55
56
      op->addZOT(std::forward<Args>(args)...);
      return op;
    }
57
58


59
60
    /// Add coefficients for first-order operator < psi, b * grad(phi) > or
    /// < grad(psi), b * phi >, where the first corresponds to
61
    /// \p type == GRD_PHI and the second one to \p type == GRD_PSI.
62
    /// The coefficient \p b must be a vector expression.
63
    template <class Term>
64
    Self& addFOT(Term const& b, FirstOrderType type)
65
    {
66
      return addFOTImpl(toTerm(b), type);
67
    }
68
69
70

    /// Add coefficients for first-order operator < psi, b * d_i(phi) > or
    /// < d_i(psi), b * phi >, where the first corresponds to
71
    /// \p type == GRD_PHI and the second one to \p type == GRD_PSI.
72
    /// The coefficient \p b must be a scalar expression.
73
    template <class Term>
74
    Self& addFOT(Term const& b, std::size_t i, FirstOrderType type)
75
    {
76
      return addFOTImpl(toTerm(b), i, type);
77
    }
78

79
    template <class... Args>
80
    static shared_ptr<Self> fot(Args&&... args)
81
    {
82
      auto op = make_shared<Self>();
83
84
85
      op->addFOT(std::forward<Args>(args)...);
      return op;
    }
86
87


88
    /// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
89
    /// where \p A can be a matrix or a scalar expression.
90
    template <class Term>
91
    Self& addSOT(Term const& A)
92
    {
93
      return addSOTImpl(toTerm(A));
94
    }
95

96
97
98
99
    /// Add coefficients for second-order operator < d_i(psi), A * d_j(phi) >,
    /// where \p A can be a matrix or a scalar expression. The first index \p _i
    /// corresponds to the derivative component of the test-function and the
    /// second index \p _j to the derivative component of the trial function.
100
101
    template <class Term>
    Self& addSOT(Term const& A, std::size_t i, std::size_t j)
102
    {
103
      return addSOTImpl(toTerm(A), i, j);
104
    }
105

106
    template <class... Args>
107
    static shared_ptr<Self> sot(Args&&... args)
108
    {
109
      auto op = make_shared<Self>();
110
111
112
      op->addSOT(std::forward<Args>(args)...);
      return op;
    }
113
114


115
  public: // initialize and assemble operator on element
116

117
118
119
120
    /// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
    template <class RowFeSpace, class ColFeSpace>
    void init(RowFeSpace const& rowFeSpace,
              ColFeSpace const& colFeSpace);
121

122

123
124
125
    /// Calculates the needed quadrature degree for the given term-order \p order.
    /// This takes into account the polynomial degree of trial and test functions
    /// and the polynomial degree of the coefficient functions
126
    template <class Geometry>
127
128
    int getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order,
                            FirstOrderType firstOrderType = GRD_PHI);
129
130


131
    template <class RowView, class ColView>
132
    bool getElementMatrix(LocalContext const& element,
133
                          RowView const& rowView,
134
135
136
                          ColView const& colView,
                          ElementMatrix& elementMatrix,
                          double* factor = NULL);
137

138
    template <class RowView>
139
    bool getElementVector(LocalContext const& element,
140
                          RowView const& rowView,
141
142
                          ElementVector& elementVector,
                          double* factor = NULL);
143

144
145
146

  private: // implementation details

147
148
    template <class Term>
    Self& addZOTImpl(Term const& term);
149

150
151
    template <class Term>
    Self& addFOTImpl(Term const& term, FirstOrderType firstOrderType);
152

153
154
    template <class Term>
    Self& addFOTImpl(Term const& term, std::size_t i, FirstOrderType firstOrderType);
155

156
157
    template <class Term>
    Self& addSOTImpl(Term const& term);
158

159
160
    template <class Term>
    Self& addSOTImpl(Term const& term, std::size_t i, std::size_t j);
161
162


163
  protected: // sum up constribution from all operator-terms
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193

    template <class... Args>
    void evalZot(double& result, int iq, double factor, Args const&... args)
    {
      for (auto* operatorTerm : zeroOrder)
        result += operatorTerm->evalZot(iq, args...) * factor;
    }

    template <class... Args>
    void evalFot1(double& result, int iq, double factor, Args const&... args)
    {
      for (auto* operatorTerm : firstOrderGrdPhi)
        result += operatorTerm->evalFot1(iq, args...) * factor;
    }

    template <class... Args>
    void evalFot2(double& result, int iq, double factor, Args const&... args)
    {
      for (auto* operatorTerm : firstOrderGrdPsi)
        result += operatorTerm->evalFot2(iq, args...) * factor;
    }

    template <class... Args>
    void evalSot(double& result, int iq, double factor, Args const&... args)
    {
      for (auto* operatorTerm : secondOrder)
        result += operatorTerm->evalSot(iq, args...) * factor;
    }


194
195
  protected: // initialize opoerator-terms

196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
    template <class E, class Quadrature>
    void initZot(E const& element, Quadrature const& quad)
    {
      for (auto* operatorTerm : zeroOrder)
        operatorTerm->init(element, quad);
    }

    template <class E, class Quadrature>
    void initFot1(E const& element, Quadrature const& quad)
    {
      for (auto* operatorTerm : firstOrderGrdPhi)
        operatorTerm->init(element, quad);
    }

    template <class E, class Quadrature>
    void initFot2(E const& element, Quadrature const& quad)
    {
      for (auto* operatorTerm : firstOrderGrdPsi)
        operatorTerm->init(element, quad);
    }

    template <class E, class Quadrature>
    void initSot(E const& element, Quadrature const& quad)
    {
      for (auto* operatorTerm : secondOrder)
        operatorTerm->init(element, quad);
    }

224

225
  private:
226

227
228
    template <class LocalView>
    IdType getElementId(LocalView const& localView);
229

230

231
  private:
232

233
234
    /// List of all second order terms
    std::list<OperatorTermType*> secondOrder;
235

236
237
    /// List of all first order terms derived to psi
    std::list<OperatorTermType*> firstOrderGrdPsi;
238

239
240
    /// List of all first order terms derived to phi
    std::list<OperatorTermType*> firstOrderGrdPhi;
241

242
243
    /// List of all zero order terms
    std::list<OperatorTermType*> zeroOrder;
244

245
246
    int psiDegree = 1;
    int phiDegree = 1;
247

248
249
    IdType lastMatrixId = 0;
    IdType lastVectorId = 0;
250

251
252
    ElementMatrix cachedElementMatrix;
    ElementVector cachedElementVector;
253
  };
254

255
256
257
} // end namespace AMDiS

#include "Operator.inc.hpp"