Operator.h 13.2 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
25
// ==                                                                        ==
// ============================================================================

/** \file Operator.h */

#ifndef AMDIS_OPERATOR_H
#define AMDIS_OPERATOR_H

#include <vector>
26
#include "AMDiS_fwd.h"
27
28
29
30
31
#include "FixVec.h"
#include "Flag.h"
#include "MatrixVector.h"
#include "ElInfo.h"
#include "AbstractFunction.h"
32
#include "OpenMP.h"
33
#include "SubAssembler.h"
34
35
36
37
#include "OperatorTerm.h"
#include "ZeroOrderTerm.h"
#include "FirstOrderTerm.h"
#include "SecondOrderTerm.h"
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

namespace AMDiS {

  /** \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:
54
55
    /// Obsolete consructor. Calling this constructor yields an error. Will be removed
    /// in on of the next revisions.
56
    Operator(Flag operatorType,
57
58
	     const FiniteElemSpace *rowFeSpace,
	     const FiniteElemSpace *colFeSpace = NULL);
59

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

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

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

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

79
    /// Adds a SecondOrderTerm to the Operator
80
81
    template <typename T>
    void addSecondOrderTerm(T *term);
82

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

    /// Adds a ZeroOrderTerm to the Operator
88
89
    template <typename T>
    void addZeroOrderTerm(T *term);
90
91
92
93
94
95
96


    /** \brief
     * Calculates the element matrix for this ElInfo and adds it multiplied by
     * factor to userMat.
     */
    virtual void getElementMatrix(const ElInfo *elInfo, 
97
				  ElementMatrix& userMat, 
98
99
				  double factor = 1.0);

100
101
    virtual void getElementMatrix(const ElInfo *rowElInfo,
				  const ElInfo *colElInfo,
102
103
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
104
105
				  bool rowColFeSpaceEqual,
				  ElementMatrix& userMat,				  
106
107
				  double factor = 1.0);

108
109
110
111
112
    /** \brief
     * Calculates the element vector for this ElInfo and adds it multiplied by
     * factor to userVec.
     */
    virtual void getElementVector(const ElInfo *elInfo, 
113
				  ElementVector& userVec, 
114
115
				  double factor = 1.0);

Thomas Witkowski's avatar
Thomas Witkowski committed
116
117
118
119
    virtual void getElementVector(const ElInfo *mainElInfo, 
				  const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
120
				  ElementVector& userVec,
Thomas Witkowski's avatar
Thomas Witkowski committed
121
				  double factor = 1.0);
122

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

126
127
    /// Returns \ref rowFeSpace.
    inline const FiniteElemSpace *getRowFeSpace() 
128
    { 
129
      return rowFeSpace; 
130
    }
131

132
133
    /// Returns \ref colFeSpace.
    inline const FiniteElemSpace *getColFeSpace() 
134
    { 
135
      return colFeSpace; 
136
    }
137

138
139
    /// Returns \ref auxFeSpaces.
    inline std::set<const FiniteElemSpace*>& getAuxFeSpaces()
140
    {
141
      return auxFeSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
142
143
    }

144
    /// Sets \ref uhOld.
145
146
    void setUhOld(const DOFVectorBase<double> *uhOld);

147
    /// Returns \ref uhOld.
148
149
    inline const DOFVectorBase<double> *getUhOld() 
    {
150
      return uhOld;
151
    }    
152

Thomas Witkowski's avatar
Thomas Witkowski committed
153
    /// Returns \ref assembler
154
    Assembler *getAssembler(int rank = 0);
155

Thomas Witkowski's avatar
Thomas Witkowski committed
156
    /// Sets \ref assembler
157
    void setAssembler(int rank, Assembler *ass);
158

Thomas Witkowski's avatar
Thomas Witkowski committed
159
    /// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
160
161
    inline void setFillFlag(Flag f) 
    { 
162
      fillFlag = f; 
163
    }
164

Thomas Witkowski's avatar
Thomas Witkowski committed
165
    /// Sets \ref needDualTraverse.
166
167
    void setNeedDualTraverse(bool b) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
168
169
170
171
      needDualTraverse = b;
    }

    /// Returns \ref fillFlag
172
173
    inline Flag getFillFlag() 
    { 
174
      return fillFlag; 
175
    }
176

Thomas Witkowski's avatar
Thomas Witkowski committed
177
    /// Returns \ref needDualTraverse
178
179
    bool getNeedDualTraverse() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
180
181
182
183
      return needDualTraverse;
    }

    /// Initialization of the needed SubAssemblers using the given quadratures. 
184
185
    void initAssembler(int rank,
		       Quadrature *quad2,
186
187
188
189
190
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


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

Thomas Witkowski's avatar
Thomas Witkowski committed
194
    /// Evaluation of all terms in \ref zeroOrder. 
195
196
    void evalZeroOrder(int nPoints,		       
		       const ElementVector& uhAtQP,
197
198
		       const WorldVector<double> *grdUhAtQP,
		       const WorldMatrix<double> *D2UhAtQP,
199
		       ElementVector& result,
200
201
		       double factor) const
    {
202
203
      int myRank = omp_get_thread_num();

204
      std::vector<OperatorTerm*>::const_iterator termIt;
205
206
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
207
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
208
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
209
    }
210

211
    /// Evaluation of all terms in \ref firstOrderGrdPsi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
212
    void evalFirstOrderGrdPsi(int nPoints,
213
			      const ElementVector& uhAtQP,
214
215
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
216
			      ElementVector& result,
217
218
			      double factor) const
    {
219
220
      int myRank = omp_get_thread_num();

221
      std::vector<OperatorTerm*>::const_iterator termIt;
222
223
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
224
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
225
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
226
    }
227

228
    /// Evaluation of all terms in \ref firstOrderGrdPhi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
229
    void evalFirstOrderGrdPhi(int nPoints,
230
			      const ElementVector& uhAtQP,
231
232
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
233
			      ElementVector& result,
234
235
			      double factor) const
    {
236
237
      int myRank = omp_get_thread_num();

238
      std::vector<OperatorTerm*>::const_iterator termIt;
239
240
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
241
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
242
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
243
    }
244

245
    /// Evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
246
    void evalSecondOrder(int nPoints,
247
			 const ElementVector& uhAtQP,
248
249
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
250
			 ElementVector& result,
251
252
			 double factor) const
    {
253
254
      int myRank = omp_get_thread_num();

255
      std::vector<OperatorTerm*>::const_iterator termIt;
256
257
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
258
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
259
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
260
    }
261

262
    /// Weak evaluation of all terms in \ref secondOrder. 
263
    void weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
Thomas Witkowski's avatar
Thomas Witkowski committed
264
			     std::vector<WorldVector<double> > &result) const;
265
  
266
    /// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
Thomas Witkowski's avatar
Thomas Witkowski committed
267
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
268
    {
269
270
      int myRank = omp_get_thread_num();

271
      std::vector<OperatorTerm*>::const_iterator termIt;
272
273
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
274
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
275
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
Thomas Witkowski's avatar
Thomas Witkowski committed
276
    }
277
  
278
279
280
281
    /// Calls getLb() for each term in \ref firstOrderGrdPsi and adds the results to Lb.
    void getLbGrdPsi(const ElInfo *elInfo, 
		     int nPoints, 
		     VectorOfFixVecs<DimVec<double> >& Lb) const
282
    {
283
284
      int myRank = omp_get_thread_num();

285
      std::vector<OperatorTerm*>::const_iterator termIt;
286
287
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
288
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
289
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
290
    }
291

292
293
294
295
    /// Calls getLb() for each term in \ref firstOrderGrdPhi and adds the results to Lb.
    void getLbGrdPhi(const ElInfo *elInfo, 
		     int nPoints, 
		     VectorOfFixVecs<DimVec<double> >& Lb) const
296
    {
297
298
      int myRank = omp_get_thread_num();

299
      std::vector<OperatorTerm*>::const_iterator termIt;
300
301
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
302
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
303
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
304
    }
305

306
    /// Calls getC() for each term in \ref zeroOrder and adds the results to c.
307
    void getC(const ElInfo *elInfo, int nPoints, ElementVector& c)
308
    {
309
310
      int myRank = omp_get_thread_num();

311
      std::vector<OperatorTerm*>::const_iterator termIt;
312
313
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
314
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
315
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
Thomas Witkowski's avatar
Thomas Witkowski committed
316
    }
317

Thomas Witkowski's avatar
Thomas Witkowski committed
318
    /// Returns true, if there are second order terms. Returns false otherwise.
319
320
    inline bool secondOrderTerms() 
    {
321
      return secondOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
322
    }
323

324
    /// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
325
326
    inline bool firstOrderTermsGrdPsi() 
    {
327
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
328
    }
329

330
    /// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
331
332
    inline bool firstOrderTermsGrdPhi() 
    {
333
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
334
    }
335

336
    /// Returns true, if there are zero order terms. Returns false otherwise.
337
338
    inline bool zeroOrderTerms() 
    {
339
      return zeroOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
340
    }
341
342

  public:
343
    /// Constant type flag for matrix operators. Obsolete, will be removed!
344
345
    static const Flag MATRIX_OPERATOR;

346
    /// Constant type flag for vector operators. Obsolete, will be removed!
347
348
349
    static const Flag VECTOR_OPERATOR;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
350
    /// FiniteElemSpace for matrix rows and element vector
351
    const FiniteElemSpace *rowFeSpace;
352

353
354
    /// FiniteElemSpace for matrix columns. Can be equal to \rowFeSpace.
    const FiniteElemSpace *colFeSpace;
355

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

    /// Number of rows in the element matrix
360
361
    int nRow;

Thomas Witkowski's avatar
Thomas Witkowski committed
362
    /// Number of columns in the element matrix
363
364
    int nCol;

Thomas Witkowski's avatar
Thomas Witkowski committed
365
    /// Flag for mesh traversal
366
367
    Flag fillFlag;

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

371
372
373
374
375
    /** \brief
     * 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.
     */
376
    std::vector<Assembler*> assembler;
377

Thomas Witkowski's avatar
Thomas Witkowski committed
378
    /// List of all second order terms
379
    std::vector< std::vector<OperatorTerm*> > secondOrder;
380

Thomas Witkowski's avatar
Thomas Witkowski committed
381
    /// List of all first order terms derived to psi
382
    std::vector< std::vector<OperatorTerm*> > firstOrderGrdPsi;
383

Thomas Witkowski's avatar
Thomas Witkowski committed
384
    /// List of all first order terms derived to phi
385
    std::vector< std::vector<OperatorTerm*> > firstOrderGrdPhi;
386

Thomas Witkowski's avatar
Thomas Witkowski committed
387
    /// List of all zero order terms
388
    std::vector< std::vector<OperatorTerm*> > zeroOrder;
389
390

    /** \brief
391
392
     * 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.
393
394
395
     */
    const DOFVectorBase<double> *uhOld;

Thomas Witkowski's avatar
Thomas Witkowski committed
396
    /// Spezifies whether optimized assemblers are used or not.
397
398
399
400
401
402
403
404
405
406
407
    bool optimized;

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

}

408
409
#include "Operator.hh"

410
#endif // AMDIS_OPERATOR_H