Operator.hpp 6.99 KB
Newer Older
1
2
3
4
5
#pragma once

#include <list>

#include "OperatorTerm.hpp"
6
#include "TermGenerator.hpp"
7
8
9

namespace AMDiS
{
10
11
  enum FirstOrderType
  {
12
13
    GRD_PSI,
    GRD_PHI
14
15
16
  };
    
    
17
18
19
20
21
22
  template <class MeshView>
  class Operator
  {
    using Self = Operator;
    using OperatorTermType = OperatorTerm<MeshView>;
    
23
24
25
  public:    
    /// Add coefficients for zero-order operator < c * u, v >.
    /// The coefficient \p c must be a scalar expression.
26
    template <class Term>
27
    Self& addZOT(Term const& c)
28
    {
29
      return addZOTImpl(toTerm(c));
30
31
    }
    
32
33
34
35
36
37
38
39
    template <class... Args>
    static std::shared_ptr<Self> zot(Args&&... args)
    {
      auto op = std::make_shared<Self>();
      op->addZOT(std::forward<Args>(args)...);
      return op;
    }
    
40
    
41
42
43
44
    /// 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.
45
    template <class Term>
46
    Self& addFOT(Term const& b, FirstOrderType firstOrderType)
47
    {
48
      return addFOTImpl(toTerm(b), firstOrderType);
49
50
    }
    
51
52
53
54
    /// 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.
55
    template <class Term, size_t I>
56
    Self& addFOT(Term const& b, const index_<I> _i, FirstOrderType firstOrderType)
57
    {
58
      return addFOTImpl(toTerm(b), _i, firstOrderType);
59
60
    }
    
61
62
63
64
65
66
67
68
    template <class... Args>
    static std::shared_ptr<Self> fot(Args&&... args)
    {
      auto op = std::make_shared<Self>();
      op->addFOT(std::forward<Args>(args)...);
      return op;
    }
    
69
    
70
71
    /// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
    /// where \p A can be a matrix or a scalar expression. 
72
    template <class Term>
73
    Self& addSOT(Term const& A)
74
    {
75
      return addSOTImpl(toTerm(A));
76
77
    }
    
78
79
80
81
    /// 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.
82
    template <class Term, size_t I, size_t J>
83
    Self& addSOT(Term const& A, const index_<I> _i, const index_<J> _j)
84
    {
85
      return addSOTImpl(toTerm(A), _i, _j);
86
87
    }
    
88
89
90
91
92
93
94
    template <class... Args>
    static std::shared_ptr<Self> sot(Args&&... args)
    {
      auto op = std::make_shared<Self>();
      op->addSOT(std::forward<Args>(args)...);
      return op;
    }
95
    
96
97
98
99
100
    
    /// Calls \ref zot(), \ref for() or \ref sot(), depending on template 
    /// parameter \p Order.
    template <size_t Order, class... Args>
    static std::shared_ptr<Self> create(Args&&... args)
101
    {
102
      return create(index_<Order>{}, std::forward<Args>(args)...);
103
104
105
    }
    
    
106
107
108
109
110
    /// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
    template <class RowFeSpace, class ColFeSpace>
    void init(RowFeSpace const& rowFeSpace,
              ColFeSpace const& colFeSpace);
    
111

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

205
206
    /// List of all first order terms derived to psi
    std::list<OperatorTermType*> firstOrderGrdPsi;
207

208
209
    /// List of all first order terms derived to phi
    std::list<OperatorTermType*> firstOrderGrdPhi;
210

211
212
213
214
215
    /// List of all zero order terms
    std::list<OperatorTermType*> zeroOrder;
    
    int psiDegree = 1;
    int phiDegree = 1;
216
217
218
219
220
  };
  
} // end namespace AMDiS

#include "Operator.inc.hpp"