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

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

6
7
8
9
10
#include <dune/amdis/OperatorTermBase.hpp>
#include <dune/amdis/terms/TermGenerator.hpp>

#include <dune/amdis/common/TypeDefs.hpp>
#include <dune/amdis/utility/GetDegree.hpp>
11
12
13

namespace AMDiS
{
14
15
  enum FirstOrderType
  {
16
17
    GRD_PSI,
    GRD_PHI
18
19
20
  };
    
    
21
22
23
24
25
  template <class MeshView>
  class Operator
  {
    using Self = Operator;
    using OperatorTermType = OperatorTerm<MeshView>;
26
27
28
29
30

    using ElementVector = Impl::ElementVector;
    using ElementMatrix = Impl::ElementMatrix;
    
    using IdType = typename MeshView::Grid::LocalIdSet::IdType;
31
    
32
33
34
  public:    
    /// Add coefficients for zero-order operator < c * u, v >.
    /// The coefficient \p c must be a scalar expression.
35
    template <class Term>
36
    Self& addZOT(Term const& c)
37
    {
38
      return addZOTImpl(toTerm(c));
39
40
    }
    
41
    template <class... Args>
42
    static shared_ptr<Self> zot(Args&&... args)
43
    {
44
      auto op = make_shared<Self>();
45
46
47
48
      op->addZOT(std::forward<Args>(args)...);
      return op;
    }
    
49
    
50
51
52
53
    /// Add coefficients for first-order operator < psi, term * grad(phi) > or 
    /// < grad(psi), term * phi >, where the first corresponds to 
    /// \p firstOrderType == GRD_PHI and the second one to \p firstOrderType == GRD_PSI.
    /// The coefficient \p b must be a vector expression.
54
    template <class Term>
55
    Self& addFOT(Term const& b, FirstOrderType firstOrderType)
56
    {
57
      return addFOTImpl(toTerm(b), firstOrderType);
58
59
    }
    
60
61
62
63
    /// Add coefficients for first-order operator < psi, b * d_i(phi) > or 
    /// < d_i(psi), b * phi >, where the first corresponds to 
    /// \p firstOrderType == GRD_PHI and the second one to \p firstOrderType == GRD_PSI.
    /// The coefficient \p b must be a scalar expression.
64
    template <class Term, size_t I>
65
    Self& addFOT(Term const& b, const index_<I> _i, FirstOrderType firstOrderType)
66
    {
67
      return addFOTImpl(toTerm(b), _i, firstOrderType);
68
69
    }
    
70
    template <class... Args>
71
    static shared_ptr<Self> fot(Args&&... args)
72
    {
73
      auto op = make_shared<Self>();
74
75
76
77
      op->addFOT(std::forward<Args>(args)...);
      return op;
    }
    
78
    
79
80
    /// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
    /// where \p A can be a matrix or a scalar expression. 
81
    template <class Term>
82
    Self& addSOT(Term const& A)
83
    {
84
      return addSOTImpl(toTerm(A));
85
86
    }
    
87
88
89
90
    /// 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.
91
    template <class Term, size_t I, size_t J>
92
    Self& addSOT(Term const& A, const index_<I> _i, const index_<J> _j)
93
    {
94
      return addSOTImpl(toTerm(A), _i, _j);
95
96
    }
    
97
    template <class... Args>
98
    static shared_ptr<Self> sot(Args&&... args)
99
    {
100
      auto op = make_shared<Self>();
101
102
103
      op->addSOT(std::forward<Args>(args)...);
      return op;
    }
104
    
105
106
107
108
    
    /// Calls \ref zot(), \ref for() or \ref sot(), depending on template 
    /// parameter \p Order.
    template <size_t Order, class... Args>
109
    static shared_ptr<Self> create(Args&&... args)
110
    {
111
      return create(index_<Order>{}, std::forward<Args>(args)...);
112
113
114
    }
    
    
115
116
117
118
119
    /// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
    template <class RowFeSpace, class ColFeSpace>
    void init(RowFeSpace const& rowFeSpace,
              ColFeSpace const& colFeSpace);
    
120

121
122
123
    /// 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
124
    int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);
125
126
      
    
127
    template <class RowView, class ColView>
128
129
130
131
132
    bool getElementMatrix(RowView const& rowView,
                          ColView const& colView,
                          ElementMatrix& elementMatrix,
                          double* factor = NULL);
    
133
    
134
    template <class RowView>
135
136
137
    bool getElementVector(RowView const& rowView,
                          ElementVector& elementVector,
                          double* factor = NULL);
138
139
    
  protected:
140
    template <class RowView, class ColView>
141
142
    void assembleZeroOrderTerms(RowView const& rowView,
				ColView const& colView,
143
				ElementMatrix& elementMatrix);
144
    
145
    template <class RowView>
146
    void assembleZeroOrderTerms(RowView const& rowView,
147
				ElementVector& elementVector);
148
    
149
    template <class RowView, class ColView>
150
151
    void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
				       ColView const& colView,
152
				       ElementMatrix& elementMatrix);
153
    
154
    template <class RowView, class ColView>
155
156
    void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
				       ColView const& colView,
157
				       ElementMatrix& elementMatrix);
158
    
159
    template <class RowView>
160
    void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
161
                                       ElementVector& elementVector);
162
    
163
    template <class RowView, class ColView>
164
165
    void assembleSecondOrderTerms(RowView const& rowView,
				  ColView const& colView,
166
				  ElementMatrix& elementMatrix);
167
    
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
  private:
    template <class Term>
    Self& addZOTImpl(Term const& term);
    
    template <class Term>
    Self& addFOTImpl(Term const& term, FirstOrderType firstOrderType);
    
    template <class Term, size_t I>
    Self& addFOTImpl(Term const& term, const index_<I>, 
                     FirstOrderType firstOrderType);
    
    template <class Term>
    Self& addSOTImpl(Term const& term);
    
    template <class Term, size_t I, size_t J>
    Self& addSOTImpl(Term const& term, const index_<I>, const index_<J>);
    
185
    
186
    template <class... Args>
187
    static shared_ptr<Self> create(index_<0>, Args&&... args)
188
189
190
191
192
    {
      return zot(std::forward<Args>(args)...);
    }
    
    template <class... Args>
193
    static shared_ptr<Self> create(index_<1>, Args&&... args)
194
195
196
197
198
    {
      return fot(std::forward<Args>(args)...);
    }
    
    template <class... Args>
199
    static shared_ptr<Self> create(index_<2>, Args&&... args)
200
201
202
203
    {
      return sot(std::forward<Args>(args)...);
    }
      
204
205
206
  private:
    template <class LocalView>
    IdType getElementId(LocalView const& localView);
207
      
208
  private:
209
210
    /// List of all second order terms
    std::list<OperatorTermType*> secondOrder;
211

212
213
    /// List of all first order terms derived to psi
    std::list<OperatorTermType*> firstOrderGrdPsi;
214

215
216
    /// List of all first order terms derived to phi
    std::list<OperatorTermType*> firstOrderGrdPhi;
217

218
219
220
221
222
    /// List of all zero order terms
    std::list<OperatorTermType*> zeroOrder;
    
    int psiDegree = 1;
    int phiDegree = 1;
223
224
225
226
227
228
    
    IdType lastMatrixId = 0;
    IdType lastVectorId = 0;
    
    ElementMatrix cachedElementMatrix;
    ElementVector cachedElementVector;
229
230
231
232
233
  };
  
} // end namespace AMDiS

#include "Operator.inc.hpp"