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"