Operator.hpp 5.73 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
    /// 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.
37
    template <class Term>
38
    Self& addFOT(Term const& b, FirstOrderType firstOrderType)
39
    {
40
      return addFOTImpl(toTerm(b), firstOrderType);
41
42
    }
    
43
44
45
46
    /// 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.
47
    template <class Term, size_t I>
48
    Self& addFOT(Term const& b, const index_<I> _i, FirstOrderType firstOrderType)
49
    {
50
      return addFOTImpl(toTerm(b), _i, firstOrderType);
51
52
    }
    
53
    
54
55
    /// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
    /// where \p A can be a matrix or a scalar expression. 
56
    template <class Term>
57
    Self& addSOT(Term const& A)
58
    {
59
      return addSOTImpl(toTerm(A));
60
61
    }
    
62
63
64
65
    /// 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.
66
    template <class Term, size_t I, size_t J>
67
    Self& addSOT(Term const& A, const index_<I> _i, const index_<J> _j)
68
    {
69
      return addSOTImpl(toTerm(A), _i, _j);
70
71
    }
    
72
73
74
75
76
77
78
    
    static std::shared_ptr<Self> create()
    {
      return std::make_shared<Self>();
    }
    
    
79
80
81
82
83
    /// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
    template <class RowFeSpace, class ColFeSpace>
    void init(RowFeSpace const& rowFeSpace,
              ColFeSpace const& colFeSpace);
    
84

85
86
87
    /// 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
88
    int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);
89
90
91
92
93
94
95
96
      
    
    template <class RowView, class ColView, class ElementMatrix>
    bool getElementMatrix(RowView const& rowView,
                          ColView const& colView,
                          ElementMatrix& elementMatrix,
                          double* factor = NULL);
    
97
    
98
99
100
101
    template <class RowView, class ElementVector>
    bool getElementVector(RowView const& rowView,
                          ElementVector& elementVector,
                          double* factor = NULL);
102
103
104
105
106
    
  protected:
    template <class RowView, class ColView, class ElementMatrix>
    void assembleZeroOrderTerms(RowView const& rowView,
				ColView const& colView,
107
108
				ElementMatrix& elementMatrix,
                                double factor);
109
    
110
111
    template <class RowView, class ElementVector>
    void assembleZeroOrderTerms(RowView const& rowView,
112
113
				ElementVector& elementVector,
                                double factor);
114
    
115
116
117
    template <class RowView, class ColView, class ElementMatrix>
    void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
				       ColView const& colView,
118
119
				       ElementMatrix& elementMatrix,
                                       double factor);
120
121
122
123
    
    template <class RowView, class ColView, class ElementMatrix>
    void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
				       ColView const& colView,
124
125
				       ElementMatrix& elementMatrix,
                                       double factor);
126
    
127
128
129
130
131
    template <class RowView, class ElementVector>
    void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
                                       ElementVector& elementVector,
                                       double factor);
    
132
133
134
    template <class RowView, class ColView, class ElementMatrix>
    void assembleSecondOrderTerms(RowView const& rowView,
				  ColView const& colView,
135
136
				  ElementMatrix& elementMatrix,
                                  double factor);
137
    
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
  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>);
    
155
  private:
156
157
    /// List of all second order terms
    std::list<OperatorTermType*> secondOrder;
158

159
160
    /// List of all first order terms derived to psi
    std::list<OperatorTermType*> firstOrderGrdPsi;
161

162
163
    /// List of all first order terms derived to phi
    std::list<OperatorTermType*> firstOrderGrdPhi;
164

165
166
167
168
169
    /// List of all zero order terms
    std::list<OperatorTermType*> zeroOrder;
    
    int psiDegree = 1;
    int phiDegree = 1;
170
171
172
173
174
  };
  
} // end namespace AMDiS

#include "Operator.inc.hpp"