Operator.hpp 4.01 KB
Newer Older
1
2
3
4
5
6
7
8
#pragma once

#include <list>

#include "OperatorTerm.hpp"

namespace AMDiS
{
9
10
  enum FirstOrderType
  {
11
12
    GRD_PSI,
    GRD_PHI
13
14
15
  };
    
    
16
17
18
19
20
21
22
  template <class MeshView>
  class Operator
  {
    using Self = Operator;
    using OperatorTermType = OperatorTerm<MeshView>;
    
  public:
23
24
25
26
    template <class RowFeSpace, class ColFeSpace>
    void init(RowFeSpace const& rowFeSpace,
	      ColFeSpace const& colFeSpace);
      
27
    template <class RowView, class ColView, class ElementMatrix>
28
    bool getElementMatrix(RowView const& rowView,
29
			  ColView const& colView,
30
31
			  ElementMatrix& elementMatrix,
                          double* factor = NULL);
32
33
    
    template <class RowView, class ElementVector>
34
35
36
37
    bool getElementVector(RowView const& rowView,
			  ElementVector& elementVector,
                          double* factor = NULL);
    
38
    
39
40
    template <class Term>
    Self& addZOT(Term const& term)
41
    {
42
43
      zeroOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
      return *this;
44
45
    }
    
46
    
47
48
49
    template <class Term>
    Self& addFOT(Term const& term, FirstOrderType firstOrderType)
    {
50
51
52
53
54
55
      using OpTerm = GenericOperatorTerm<MeshView, Term>;
      if (firstOrderType == GRD_PHI)
        firstOrderGrdPhi.push_back(new OpTerm(term));
      else
        firstOrderGrdPsi.push_back(new OpTerm(term));
      return *this;
56
57
58
59
60
    }
    
    template <class Term, size_t I>
    Self& addFOT(Term const& term, const index_<I>, FirstOrderType firstOrderType)
    {
61
62
63
64
65
66
67
      using OpTerm = GenericOperatorTerm<MeshView, Term, VectorComponent<I>>;
      
      if (firstOrderType == GRD_PHI)
        firstOrderGrdPhi.push_back(new OpTerm(term));
      else
        firstOrderGrdPsi.push_back(new OpTerm(term));
      return *this;
68
69
    }
    
70
    
71
72
    template <class Term>
    Self& addSOT(Term const& term)
73
    {
74
75
      secondOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
      return *this;
76
77
    }
    
78
79
80
    template <class Term, size_t I, size_t J>
    Self& addSOT(Term const& term, const index_<I>, const index_<J>)
    {
81
82
83
      using OpTerm = GenericOperatorTerm<MeshView, Term, MatrixComponent<I,J>>;
      secondOrder.push_back(new OpTerm(term));
      return *this;
84
85
    }
    
86
87

    /// Calculates the needed quadrature degree for the given order.
88
    int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);
89
90
91
92
93
94
    
    
  protected:
    template <class RowView, class ColView, class ElementMatrix>
    void assembleZeroOrderTerms(RowView const& rowView,
				ColView const& colView,
95
96
				ElementMatrix& elementMatrix,
                                double factor);
97
    
98
99
    template <class RowView, class ElementVector>
    void assembleZeroOrderTerms(RowView const& rowView,
100
101
				ElementVector& elementVector,
                                double factor);
102
    
103
104
105
    template <class RowView, class ColView, class ElementMatrix>
    void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
				       ColView const& colView,
106
107
				       ElementMatrix& elementMatrix,
                                       double factor);
108
109
110
111
    
    template <class RowView, class ColView, class ElementMatrix>
    void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
				       ColView const& colView,
112
113
				       ElementMatrix& elementMatrix,
                                       double factor);
114
    
115
116
117
    template <class RowView, class ColView, class ElementMatrix>
    void assembleSecondOrderTerms(RowView const& rowView,
				  ColView const& colView,
118
119
				  ElementMatrix& elementMatrix,
                                  double factor);
120
121
    
  private:
122
123
    /// List of all second order terms
    std::list<OperatorTermType*> secondOrder;
124

125
126
    /// List of all first order terms derived to psi
    std::list<OperatorTermType*> firstOrderGrdPsi;
127

128
129
    /// List of all first order terms derived to phi
    std::list<OperatorTermType*> firstOrderGrdPhi;
130

131
132
133
134
135
    /// List of all zero order terms
    std::list<OperatorTermType*> zeroOrder;
    
    int psiDegree = 1;
    int phiDegree = 1;
136
137
138
139
140
  };
  
} // end namespace AMDiS

#include "Operator.inc.hpp"