Operator.h 14.3 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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

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.
   * An Operator can by of type MATRIX_OPERATOR, if it's used to assemble the
   * system matrix on the left hand side, or it can be of type VECTOR_OPERATOR,
   * if it's used to assemble the right hand side vector. If an Operator gives
   * contributions to both sides of the system it is a MATRIX_OPERATOR and a
   * VECTOR_OPERATOR in one instance. This allows to efficiently reuse element 
   * matrices once calculated.
   * 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:
    /** \brief
     * Constructs an empty Operator of type operatorType for the given 
     * FiniteElemSpace. 
     * The type is given by a Flag that can contain the values MATRIX_OPERATOR,
     * VECTOR_OPERATOR, or MATRIX_OPERATOR | VECTOR_OPERATOR. This type specifies 
     * whether the Operator is used on the left hand side, the right hand side,
     * or on both sides of the system. 
     */
    Operator(Flag operatorType,
Thomas Witkowski's avatar
Thomas Witkowski committed
69
70
	     const FiniteElemSpace *rowFESpace,
	     const FiniteElemSpace *colFESpace = NULL);
71

72
    /// Destructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
73
    virtual ~Operator() {}
74

75
    /// Sets \ref optimized.
76
77
    inline void useOptimizedAssembler(bool opt) 
    {
78
      optimized = opt;
Thomas Witkowski's avatar
Thomas Witkowski committed
79
    }
80

81
    /// Returns \ref optimized.
82
83
    inline bool isOptimized() 
    {
84
      return optimized;
Thomas Witkowski's avatar
Thomas Witkowski committed
85
    }
86

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

91
    /// Adds a FirstOrderTerm to the Operator
92
    template <typename T>
93
94
95
    void addFirstOrderTerm(T *term, FirstOrderType type = GRD_PHI);

    /// Adds a ZeroOrderTerm to the Operator
96
97
    template <typename T>
    void addZeroOrderTerm(T *term);
98
99
100
101
102
103
104


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

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

115
116
117
118
119
    /** \brief
     * Calculates the element vector for this ElInfo and adds it multiplied by
     * factor to userVec.
     */
    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

Thomas Witkowski's avatar
Thomas Witkowski committed
133
    /// Returns \ref rowFESpace.
134
135
    inline const FiniteElemSpace *getRowFESpace() 
    { 
136
      return rowFESpace; 
137
    }
138

Thomas Witkowski's avatar
Thomas Witkowski committed
139
    /// Returns \ref colFESpace.
140
141
    inline const FiniteElemSpace *getColFESpace() 
    { 
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(int rank = 0);
162

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

Thomas Witkowski's avatar
Thomas Witkowski committed
166
    /// Returns whether this is a matrix operator.
167
168
    inline const bool isMatrixOperator() 
    {
169
      return type.isSet(MATRIX_OPERATOR);
170
    }
171

Thomas Witkowski's avatar
Thomas Witkowski committed
172
    /// Returns whether this is a vector operator
173
174
    inline const bool isVectorOperator() 
    {
175
      return type.isSet(VECTOR_OPERATOR);
176
    }
177

Thomas Witkowski's avatar
Thomas Witkowski committed
178
    /// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
179
180
    inline void setFillFlag(Flag f) 
    { 
181
      fillFlag = f; 
182
    }
183

Thomas Witkowski's avatar
Thomas Witkowski committed
184
    /// Sets \ref needDualTraverse.
185
186
    void setNeedDualTraverse(bool b) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
189
190
      needDualTraverse = b;
    }

    /// Returns \ref fillFlag
191
192
    inline Flag getFillFlag() 
    { 
193
      return fillFlag; 
194
    }
195

Thomas Witkowski's avatar
Thomas Witkowski committed
196
    /// Returns \ref needDualTraverse
197
198
    bool getNeedDualTraverse() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
201
202
      return needDualTraverse;
    }

    /// Initialization of the needed SubAssemblers using the given quadratures. 
203
204
    void initAssembler(int rank,
		       Quadrature *quad2,
205
206
207
208
209
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


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

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

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

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

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

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

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

264
    /// Evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
265
266
    void evalSecondOrder(int nPoints,
			 const double *uhAtQP,
267
268
269
270
271
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
    {
272
273
      int myRank = omp_get_thread_num();

274
      std::vector<OperatorTerm*>::const_iterator termIt;
275
276
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
277
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
278
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    }
280

281
    /// Weak evaluation of all terms in \ref secondOrder. 
282
283
    void weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
			     std::vector<WorldVector<double> > &result) const
284
    {
285
286
      int myRank = omp_get_thread_num();

287
      for (std::vector<OperatorTerm*>::const_iterator termIt = secondOrder[myRank].begin(); 
288
	   termIt != secondOrder[myRank].end(); 
289
	   ++termIt)
290
	static_cast<SecondOrderTerm*>(*termIt)->weakEval(grdUhAtQP, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
291
    }
292
  
293
    /// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
Thomas Witkowski's avatar
Thomas Witkowski committed
294
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
295
    {
296
297
      int myRank = omp_get_thread_num();

298
      std::vector<OperatorTerm*>::const_iterator termIt;
299
300
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
301
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
302
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
Thomas Witkowski's avatar
Thomas Witkowski committed
303
    }
304
  
305
306
307
308
    /// 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
309
    {
310
311
      int myRank = omp_get_thread_num();

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

319
320
321
322
    /// 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
323
    {
324
325
      int myRank = omp_get_thread_num();

326
      std::vector<OperatorTerm*>::const_iterator termIt;
327
328
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
329
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
330
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
331
    }
332

333
    /// Calls getC() for each term in \ref zeroOrder and adds the results to c.
Thomas Witkowski's avatar
Thomas Witkowski committed
334
    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c)
335
    {
336
337
      int myRank = omp_get_thread_num();

338
      std::vector<OperatorTerm*>::const_iterator termIt;
339
340
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
341
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
342
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
Thomas Witkowski's avatar
Thomas Witkowski committed
343
    }
344

Thomas Witkowski's avatar
Thomas Witkowski committed
345
    /// Returns true, if there are second order terms. Returns false otherwise.
346
347
    inline bool secondOrderTerms() 
    {
348
      return secondOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
349
    }
350

351
    /// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
352
353
    inline bool firstOrderTermsGrdPsi() 
    {
354
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
355
    }
356

357
    /// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
358
359
    inline bool firstOrderTermsGrdPhi() 
    {
360
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
361
    }
362

363
    /// Returns true, if there are zero order terms. Returns false otherwise.
364
365
    inline bool zeroOrderTerms() 
    {
366
      return zeroOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
367
    }
368
369

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
370
    /// Constant type flag for matrix operators
371
372
    static const Flag MATRIX_OPERATOR;

Thomas Witkowski's avatar
Thomas Witkowski committed
373
    /// Constant type flag for vector operators
374
375
376
    static const Flag VECTOR_OPERATOR;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
377
    /// FiniteElemSpace for matrix rows and element vector
378
379
    const FiniteElemSpace *rowFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
380
    /// FiniteElemSpace for matrix columns. Can be equal to \rowFESpace.
381
382
    const FiniteElemSpace *colFESpace;

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

    /// Number of rows in the element matrix
387
388
    int nRow;

Thomas Witkowski's avatar
Thomas Witkowski committed
389
    /// Number of columns in the element matrix
390
391
    int nCol;

Thomas Witkowski's avatar
Thomas Witkowski committed
392
    /// Type of this Operator.
393
394
    Flag type;

Thomas Witkowski's avatar
Thomas Witkowski committed
395
    /// Flag for mesh traversal
396
397
    Flag fillFlag;

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

401
402
403
404
405
    /** \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.
     */
406
    std::vector<Assembler*> assembler;
407

Thomas Witkowski's avatar
Thomas Witkowski committed
408
    /// List of all second order terms
409
    std::vector< std::vector<OperatorTerm*> > secondOrder;
410

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
417
    /// List of all zero order terms
418
    std::vector< std::vector<OperatorTerm*> > zeroOrder;
419
420
421
422
423
424
425
426
427

    /** \brief
     * Pointer to the solution of the last timestep. Can be used if the 
     * Operator is MATRIX_OPERATOR and VECTOR_OPERATOR for a more efficient
     * assemblage of the element vector when the element matrix was already
     * computed.
     */
    const DOFVectorBase<double> *uhOld;

Thomas Witkowski's avatar
Thomas Witkowski committed
428
    /// Spezifies whether optimized assemblers are used or not.
429
430
431
432
433
434
435
436
437
438
439
    bool optimized;

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

}

440
441
#include "Operator.hh"

442
#endif // AMDIS_OPERATOR_H