Operator.hpp 5.63 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
    /// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
    template <class RowFeSpace, class ColFeSpace>
    void init(RowFeSpace const& rowFeSpace,
              ColFeSpace const& colFeSpace);
    
77

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

152
153
    /// List of all first order terms derived to psi
    std::list<OperatorTermType*> firstOrderGrdPsi;
154

155
156
    /// List of all first order terms derived to phi
    std::list<OperatorTermType*> firstOrderGrdPhi;
157

158
159
160
161
162
    /// List of all zero order terms
    std::list<OperatorTermType*> zeroOrder;
    
    int psiDegree = 1;
    int phiDegree = 1;
163
164
165
166
167
  };
  
} // end namespace AMDiS

#include "Operator.inc.hpp"