Operator.h 14 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
    void weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
Thomas Witkowski's avatar
Thomas Witkowski committed
283
			     std::vector<WorldVector<double> > &result) const;
284
  
285
    /// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
Thomas Witkowski's avatar
Thomas Witkowski committed
286
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
287
    {
288
289
      int myRank = omp_get_thread_num();

290
      std::vector<OperatorTerm*>::const_iterator termIt;
291
292
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
293
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
294
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
Thomas Witkowski's avatar
Thomas Witkowski committed
295
    }
296
  
297
298
299
300
    /// 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
301
    {
302
303
      int myRank = omp_get_thread_num();

304
      std::vector<OperatorTerm*>::const_iterator termIt;
305
306
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
307
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
308
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
309
    }
310

311
312
313
314
    /// 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
315
    {
316
317
      int myRank = omp_get_thread_num();

318
      std::vector<OperatorTerm*>::const_iterator termIt;
319
320
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
321
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
322
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
323
    }
324

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

330
      std::vector<OperatorTerm*>::const_iterator termIt;
331
332
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
333
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
334
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
Thomas Witkowski's avatar
Thomas Witkowski committed
335
    }
336

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

343
    /// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
344
345
    inline bool firstOrderTermsGrdPsi() 
    {
346
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
347
    }
348

349
    /// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
350
351
    inline bool firstOrderTermsGrdPhi() 
    {
352
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
353
    }
354

355
    /// Returns true, if there are zero order terms. Returns false otherwise.
356
357
    inline bool zeroOrderTerms() 
    {
358
      return zeroOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
359
    }
360
361

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
362
    /// Constant type flag for matrix operators
363
364
    static const Flag MATRIX_OPERATOR;

Thomas Witkowski's avatar
Thomas Witkowski committed
365
    /// Constant type flag for vector operators
366
367
368
    static const Flag VECTOR_OPERATOR;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
369
    /// FiniteElemSpace for matrix rows and element vector
370
371
    const FiniteElemSpace *rowFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
372
    /// FiniteElemSpace for matrix columns. Can be equal to \rowFESpace.
373
374
    const FiniteElemSpace *colFESpace;

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

    /// Number of rows in the element matrix
379
380
    int nRow;

Thomas Witkowski's avatar
Thomas Witkowski committed
381
    /// Number of columns in the element matrix
382
383
    int nCol;

Thomas Witkowski's avatar
Thomas Witkowski committed
384
    /// Type of this Operator.
385
386
    Flag type;

Thomas Witkowski's avatar
Thomas Witkowski committed
387
    /// Flag for mesh traversal
388
389
    Flag fillFlag;

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

393
394
395
396
397
    /** \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.
     */
398
    std::vector<Assembler*> assembler;
399

Thomas Witkowski's avatar
Thomas Witkowski committed
400
    /// List of all second order terms
401
    std::vector< std::vector<OperatorTerm*> > secondOrder;
402

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
409
    /// List of all zero order terms
410
    std::vector< std::vector<OperatorTerm*> > zeroOrder;
411
412
413
414
415
416
417
418
419

    /** \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
420
    /// Spezifies whether optimized assemblers are used or not.
421
422
423
424
425
426
427
428
429
430
431
    bool optimized;

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

}

432
433
#include "Operator.hh"

434
#endif // AMDIS_OPERATOR_H