Operator.h 12.3 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20
21
22
23
24
25
26

/** \file Operator.h */

#ifndef AMDIS_OPERATOR_H
#define AMDIS_OPERATOR_H

#include <vector>
27
#include "AMDiS_fwd.h"
28
29
30
31
32
#include "FixVec.h"
#include "Flag.h"
#include "MatrixVector.h"
#include "ElInfo.h"
#include "AbstractFunction.h"
33
#include "SubAssembler.h"
34
35
36
37
#include "OperatorTerm.h"
#include "ZeroOrderTerm.h"
#include "FirstOrderTerm.h"
#include "SecondOrderTerm.h"
38
39
40

namespace AMDiS {

41
42
  using namespace std;

43
44
45
46
47
48
49
50
51
52
53
54
55
  /** \brief
   * An Operator holds all information needed to assemble the system matrix
   * and the right hand side. It consists of four OperatorTerm lists each storing
   * Terms of a specific order and type. You can define your own Operator by 
   * creating an empty Operator and than adding OperatorTerms to it.
   * By calling \ref getElementMatrix() or \ref getElementVector() one can 
   * initiate the assembling procedure. Therefor each Operator has its own
   * Assembler, especially created for this Operator, by the first call of
   * \ref getElementMatrix() or \ref getElementVector(). 
   */
  class Operator
  {
  public:
56
57
    /// Obsolete consructor. Calling this constructor yields an error. Will be removed
    /// in on of the next revisions.
58
    Operator(Flag operatorType,
59
60
	     const FiniteElemSpace *rowFeSpace,
	     const FiniteElemSpace *colFeSpace = NULL);
61

62
63
64
65
    /// Constructs an empty Operator of type operatorType for the given FiniteElemSpace.
    Operator(const FiniteElemSpace *rowFeSpace,
	     const FiniteElemSpace *colFeSpace = NULL);

66
    /// Destructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
67
    virtual ~Operator() {}
68

69
    /// Sets \ref optimized.
70
71
    inline void useOptimizedAssembler(bool opt) 
    {
72
      optimized = opt;
Thomas Witkowski's avatar
Thomas Witkowski committed
73
    }
74

75
    /// Returns \ref optimized.
76
77
    inline bool isOptimized() 
    {
78
      return optimized;
Thomas Witkowski's avatar
Thomas Witkowski committed
79
    }
80

81
    /// Adds a ZeroOrderTerm to the Operator
82
    template <typename T>
83
    void addZeroOrderTerm(T *term);
84

85
    /// Adds a FirstOrderTerm to the Operator
86
    template <typename T>
87
88
    void addFirstOrderTerm(T *term, FirstOrderType type = GRD_PHI);

89
    /// Adds a SecondOrderTerm to the Operator
90
    template <typename T>
91
    void addSecondOrderTerm(T *term);
92

93
    /// Adds a ZeroOrderTerm to the Operator
94
    void addTerm(ZeroOrderTerm *term);
95
96

    /// Adds a FirstOrderTerm to the Operator
97
98
    void addTerm(FirstOrderTerm *term,
      FirstOrderType type);
99
100

    /// Adds a SecondOrderTerm to the Operator
101
    void addTerm(SecondOrderTerm *term);
102

103
104
    /// Calculates the element matrix for this ElInfo and adds it multiplied by
    /// factor to userMat.
105
    virtual void getElementMatrix(const ElInfo *elInfo, 
106
				  ElementMatrix& userMat, 
107
108
				  double factor = 1.0);

109
110
    virtual void getElementMatrix(const ElInfo *rowElInfo,
				  const ElInfo *colElInfo,
111
112
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
113
114
				  bool rowColFeSpaceEqual,
				  ElementMatrix& userMat,				  
115
116
				  double factor = 1.0);

117
118
    /// Calculates the element vector for this ElInfo and adds it multiplied by
    /// factor to userVec.
119
    virtual void getElementVector(const ElInfo *elInfo, 
120
				  ElementVector& userVec, 
121
122
				  double factor = 1.0);

Thomas Witkowski's avatar
Thomas Witkowski committed
123
124
125
126
    virtual void getElementVector(const ElInfo *mainElInfo, 
				  const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
127
				  ElementVector& userVec,
Thomas Witkowski's avatar
Thomas Witkowski committed
128
				  double factor = 1.0);
129

Thomas Witkowski's avatar
Thomas Witkowski committed
130
131
    /// That function must be called after one assembling cycle has been finished.
    void finishAssembling();
132

133
134
    /// Returns \ref rowFeSpace.
    inline const FiniteElemSpace *getRowFeSpace() 
135
    { 
136
      return rowFeSpace; 
137
    }
138

139
140
    /// Returns \ref colFeSpace.
    inline const FiniteElemSpace *getColFeSpace() 
141
    { 
142
      return colFeSpace; 
143
    }
144

145
146
    /// Returns \ref auxFeSpaces.
    inline std::set<const FiniteElemSpace*>& getAuxFeSpaces()
147
    {
148
      return auxFeSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
149
150
    }

151
    /// Sets \ref uhOld.
152
153
    void setUhOld(const DOFVectorBase<double> *uhOld);

154
    /// Returns \ref uhOld.
155
156
    inline const DOFVectorBase<double> *getUhOld() 
    {
157
      return uhOld;
158
    }    
159

Thomas Witkowski's avatar
Thomas Witkowski committed
160
    /// Returns \ref assembler
161
    Assembler *getAssembler();
162

Thomas Witkowski's avatar
Thomas Witkowski committed
163
    /// Sets \ref assembler
164
    void setAssembler(Assembler *ass);
165

Thomas Witkowski's avatar
Thomas Witkowski committed
166
    /// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
167
168
    inline void setFillFlag(Flag f) 
    { 
169
      fillFlag = f; 
170
    }
171

Thomas Witkowski's avatar
Thomas Witkowski committed
172
    /// Sets \ref needDualTraverse.
173
174
    void setNeedDualTraverse(bool b) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
175
176
177
178
      needDualTraverse = b;
    }

    /// Returns \ref fillFlag
179
180
    inline Flag getFillFlag() 
    { 
181
      return fillFlag; 
182
    }
183

Thomas Witkowski's avatar
Thomas Witkowski committed
184
    /// Returns \ref needDualTraverse
185
186
    bool getNeedDualTraverse() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
189
190
      return needDualTraverse;
    }

    /// Initialization of the needed SubAssemblers using the given quadratures. 
191
    void initAssembler(Quadrature *quad2,
192
193
194
195
196
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


Thomas Witkowski's avatar
Thomas Witkowski committed
197
    /// Calculates the needed quadrature degree for the given order. 
198
199
    int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);

Thomas Witkowski's avatar
Thomas Witkowski committed
200
    /// Evaluation of all terms in \ref zeroOrder. 
201
    void evalZeroOrder(int nPoints,		       
202
203
204
205
		       const mtl::dense_vector<double>& uhAtQP,
		       const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
		       const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
		       mtl::dense_vector<double>& result,
206
207
		       double factor) const
    {
208
209
      for (vector<OperatorTerm*>::const_iterator termIt = zeroOrder.begin(); 
	   termIt != zeroOrder.end(); ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
210
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
211
    }
212

213
    /// Evaluation of all terms in \ref firstOrderGrdPsi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
214
    void evalFirstOrderGrdPsi(int nPoints,
215
216
217
218
			      const mtl::dense_vector<double>& uhAtQP,
			      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			      mtl::dense_vector<double>& result,
219
			      double factor) const
220
221
222
    {      
      for (vector<OperatorTerm*>::const_iterator termIt = firstOrderGrdPsi.begin(); 
	   termIt != firstOrderGrdPsi.end(); ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
223
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
224
    }
225

226
    /// Evaluation of all terms in \ref firstOrderGrdPhi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
227
    void evalFirstOrderGrdPhi(int nPoints,
228
229
230
231
			      const mtl::dense_vector<double>& uhAtQP,
			      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			      mtl::dense_vector<double>& result,
232
233
			      double factor) const
    {
234
235
      for (vector<OperatorTerm*>::const_iterator termIt = firstOrderGrdPhi.begin(); 
	   termIt != firstOrderGrdPhi.end(); ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
236
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
237
    }
238

239
    /// Evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
240
    void evalSecondOrder(int nPoints,
241
242
243
244
			  const mtl::dense_vector<double>& uhAtQP,
			  const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			  const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			  mtl::dense_vector<double>& result,
245
			 double factor) const
246
247
248
    {      
      for (vector<OperatorTerm*>::const_iterator termIt = secondOrder.begin(); 
	   termIt != secondOrder.end(); ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
249
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
250
    }
251

252
    /// Weak evaluation of all terms in \ref secondOrder. 
253
254
    void weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
			     std::vector<WorldVector<double> > &result) const;
255
  
256
    /// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
257
258
259
260
261
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const
    {      
      for (vector<OperatorTerm*>::const_iterator termIt = secondOrder.begin(); 
	   termIt != secondOrder.end(); ++termIt)
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, LALt);
Thomas Witkowski's avatar
Thomas Witkowski committed
262
    }
263
  
Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
    /// Calls getLb() for each term in \ref firstOrderGrdPsi and adds the 
    /// results to Lb.
266
    void getLbGrdPsi(const ElInfo *elInfo, 
267
		     std::vector<mtl::dense_vector<double> >& Lb) const
268
269
270
271
    {      
      for (vector<OperatorTerm*>::const_iterator termIt = firstOrderGrdPsi.begin(); 
	   termIt != firstOrderGrdPsi.end(); ++termIt)
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
272
    }
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
    /// Calls getLb() for each term in \ref firstOrderGrdPhi and adds the 
    /// results to Lb.
276
    void getLbGrdPhi(const ElInfo *elInfo, 
277
		     std::vector<mtl::dense_vector<double> >& Lb) const
278
279
280
281
    {      
      for (vector<OperatorTerm*>::const_iterator termIt = firstOrderGrdPhi.begin(); 
	   termIt != firstOrderGrdPhi.end(); ++termIt)
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
282
    }
283

284
    /// Calls getC() for each term in \ref zeroOrder and adds the results to c.
285
    void getC(const ElInfo *elInfo, int nPoints, mtl::dense_vector<double>& c)
286
287
288
    {      
      for (vector<OperatorTerm*>::const_iterator termIt = zeroOrder.begin(); 
	   termIt != zeroOrder.end(); ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
289
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
Thomas Witkowski's avatar
Thomas Witkowski committed
290
    }
291

Thomas Witkowski's avatar
Thomas Witkowski committed
292
    /// Returns true, if there are second order terms. Returns false otherwise.
293
294
    inline bool secondOrderTerms() 
    {
295
      return secondOrder.size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
296
    }
297

Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
    /// Returns true, if there are first order terms (grdPsi). Returns 
    /// false otherwise.
300
301
    inline bool firstOrderTermsGrdPsi() 
    {
302
      return firstOrderGrdPsi.size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
303
    }
304

Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
    /// Returns true, if there are first order terms (grdPhi). Returns 
    /// false otherwise.
307
308
    inline bool firstOrderTermsGrdPhi() 
    {
309
      return firstOrderGrdPhi.size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
310
    }
311

312
    /// Returns true, if there are zero order terms. Returns false otherwise.
313
314
    inline bool zeroOrderTerms() 
    {
315
      return zeroOrder.size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
316
    }
317
318

  public:
319
    /// Constant type flag for matrix operators. Obsolete, will be removed!
320
321
    static const Flag MATRIX_OPERATOR;

322
    /// Constant type flag for vector operators. Obsolete, will be removed!
323
324
325
    static const Flag VECTOR_OPERATOR;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
326
    /// FiniteElemSpace for matrix rows and element vector
327
    const FiniteElemSpace *rowFeSpace;
328

329
330
    /// FiniteElemSpace for matrix columns. Can be equal to \rowFeSpace.
    const FiniteElemSpace *colFeSpace;
331

Thomas Witkowski's avatar
Thomas Witkowski committed
332
    /// List of aux fe spaces, e.g., if a term is multiplied with a DOF vector
333
    std::set<const FiniteElemSpace*> auxFeSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335

    /// Number of rows in the element matrix
336
337
    int nRow;

Thomas Witkowski's avatar
Thomas Witkowski committed
338
    /// Number of columns in the element matrix
339
340
    int nCol;

Thomas Witkowski's avatar
Thomas Witkowski committed
341
    /// Flag for mesh traversal
342
343
    Flag fillFlag;

Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
346
    /// If true, the operator needs a dual traverse over two meshes for assembling.
    bool needDualTraverse;

Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
349
    /// Calculates the element matrix and/or the element vector. It is
    /// created especially for this Operator, when \ref getElementMatrix()
    /// or \ref getElementVector is called for the first time.
350
    Assembler* assembler;
351

Thomas Witkowski's avatar
Thomas Witkowski committed
352
    /// List of all second order terms
353
    vector<OperatorTerm*> secondOrder;
354

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    /// List of all first order terms derived to psi
356
    vector<OperatorTerm*> firstOrderGrdPsi;
357

Thomas Witkowski's avatar
Thomas Witkowski committed
358
    /// List of all first order terms derived to phi
359
    vector<OperatorTerm*> firstOrderGrdPhi;
360

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    /// List of all zero order terms
362
    vector<OperatorTerm*> zeroOrder;
363

Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
366
    /// Pointer to the solution of the last timestep. Can be used for a more
    /// efficient assemblage of the element vector when the element matrix 
    /// was already computed.
367
368
    const DOFVectorBase<double> *uhOld;

Thomas Witkowski's avatar
Thomas Witkowski committed
369
    /// Spezifies whether optimized assemblers are used or not.
370
371
372
373
374
375
376
377
378
379
380
    bool optimized;

    friend class Assembler;
    friend class SubAssembler;
    friend class ZeroOrderAssembler;
    friend class FirstOrderAssembler;
    friend class SecondOrderAssembler;
  };

}

381
382
#include "Operator.hh"

383
#endif // AMDIS_OPERATOR_H