Operator.h 12.6 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 "OpenMP.h"
34
#include "SubAssembler.h"
35
36
37
38
#include "OperatorTerm.h"
#include "ZeroOrderTerm.h"
#include "FirstOrderTerm.h"
#include "SecondOrderTerm.h"
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

272
      std::vector<OperatorTerm*>::const_iterator termIt;
273
274
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
275
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
276
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
Thomas Witkowski's avatar
Thomas Witkowski committed
277
    }
278
  
279
280
281
282
    /// 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
283
    {
284
285
      int myRank = omp_get_thread_num();

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

293
294
295
296
    /// 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
297
    {
298
299
      int myRank = omp_get_thread_num();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

}

409
410
#include "Operator.hh"

411
#endif // AMDIS_OPERATOR_H