OperatorTerm.h 8.44 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  TU Dresden                                                            ==
// ==                                                                        ==
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
// ==                                                                        ==
// ============================================================================

/** \file OperatorTerm.h */

#ifndef AMDIS_OPERATORTERM_H
#define AMDIS_OPERATORTERM_H

#include <set>
#include "AMDiS_fwd.h"
#include "SubAssembler.h"

namespace AMDiS {

  /** 
   * \ingroup Assembler
   * 
   * \brief
   * Base class for ZeroOrderTerm, FirstOrderTerm and SecondOrderTerm. 
   * OperatorTerms are the building blocks of an Operator. Each OperatorTerm
   * has its properties which are regarded, when constructing 
   * an Assembler for the corresponding Operator.
   */
  class OperatorTerm
  {
  public:
    /** \brief
     * Constructs an OperatorTerm with initially no properties.
     * degree_ is used to determine the degree of the needed quadrature
     * for the assemblage.  
     */
    OperatorTerm(int deg) 
      : properties(0), 
	degree(deg),
	dimOfWorld(Global::getGeo(WORLD)),
	bOne(-1)
    {}

    /// Destructor.
    virtual ~OperatorTerm() {}

    /** \brief
     * Virtual method. It's called by SubAssembler::initElement() for
     * each OperatorTerm belonging to this SubAssembler. E.g., vectors
     * and coordinates at quadrature points can be calculated here.
     */
    virtual void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL) 
    {}

    virtual void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
			     SubAssembler*, Quadrature *quad = NULL)
    {}


    /// Returs \auxFeSpaces, the list of all aux fe spaces the operator makes use off.
    inline std::set<const FiniteElemSpace*>& getAuxFeSpaces() 
    {
      return auxFeSpaces;
    }

    /// Specifies whether the matrix of the term is symmetric
    void setSymmetric(bool symm);

    /// Returns true, if the term is piecewise constant, returns false otherwise
    inline bool isPWConst() 
    { 
      return (degree == 0);
    }

    /// Returns true, if the term has a symmetric matrix, returns false otherwise.
    bool isSymmetric();

    /// Returns \ref degree.
    inline int getDegree() 
    { 
      return degree; 
    }

    /// Sets one component of the b vector to be one. See \ref bOne.
    void setB(int b)
    {
      bOne = b;
    }

    /// Evaluation of the OperatorTerm at all quadrature points.
    virtual void eval(int nPoints,
		      const double *uhAtQP,
		      const WorldVector<double> *grdUhAtQP,
		      const WorldMatrix<double> *D2UhAtQP,
		      double *result,
		      double factor) = 0;

    /** \brief
     * Determines the value of a dof vector at the quadrature points of a given 
     * element. It is used by all VecAtQP like operator terms.
     */
    double *getVectorAtQPs(DOFVectorBase<double>* vec,
			   const ElInfo* elInfo, 
			   SubAssembler* subAssembler,
			   Quadrature *quad);		   

    /** \brief
     * Determines the value of a dof vector at the quadrature points of a given 
     * element. This function is used, if an operator is assembled on two different
     * meshes using the dual traverse. The element, i.e. the small or the large one,
     * is choosen, which corresponds to the mesh the dof vector is defined on.
     */
    double *getVectorAtQPs(DOFVectorBase<double>* vec,
			   const ElInfo* smallElInfo, 
			   const ElInfo* largeElInfo, 
			   SubAssembler* subAssembler,
			   Quadrature *quad);	

    ///
    WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec,
					   const ElInfo* elInfo,
					   SubAssembler* subAssembler,
					   Quadrature *quad);

    ///
    WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec,
					   const ElInfo* smallElInfo, 
					   const ElInfo* largeElInfo,
					   SubAssembler* subAssembler,
					   Quadrature *quad);

  protected:
    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$.
    void lalt(const DimVec<WorldVector<double> >& Lambda,
	      const WorldMatrix<double>& matrix,
	      DimMat<double>& LALt,
	      bool symm,
	      double factor) const;

    /** \brief
     * Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for \f$ A \f$
     * the matrix having a ONE in the position \f$ (K,L) \f$
     * and ZEROS in all other positions.
     */
    static void lalt_kl(const DimVec<WorldVector<double> >& Lambda,
			int k, int l,
			DimMat<double>& LALt,
			double factor);

    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity.
    inline void l1lt(const DimVec<WorldVector<double> >& Lambda, 
		     DimMat<double>& LALt,
		     double factor) const
    {
      const int dim = LALt.getNumRows();
      
      for (int i = 0; i < dim; i++) {
	double val = 0.0;
	for (int k = 0; k < dimOfWorld; k++)
	  val += Lambda[i][k] * Lambda[i][k];
	val *= factor;
	LALt[i][i] += val;
	for (int j = i + 1; j < dim; j++) {
	  val = 0.0;
	  for (int k = 0; k < dimOfWorld; k++)
	    val += Lambda[i][k] * Lambda[j][k];
	  val *= factor;
	  LALt[i][j] += val;
	  LALt[j][i] += val;
	}    
      }      
    }

    /// Evaluation of \f$ \Lambda \cdot b\f$.
    inline void lb(const DimVec<WorldVector<double> >& Lambda,
		   const WorldVector<double>& b,
		   DimVec<double>& Lb,
		   double factor) const
    {
      const int dim = Lambda.getSize();

      for (int i = 0; i < dim; i++) {
	double val = 0.0;
	
	for (int j = 0; j < dimOfWorld; j++)
	  val += Lambda[i][j] * b[j];
	
	Lb[i] += val * factor;
      }    
    }

    /// Evaluation of \f$ \Lambda \cdot b\f$.
    inline void lb_one(const DimVec<WorldVector<double> >& Lambda,
		       DimVec<double>& Lb,
		       double factor) const
    {
      const int dim = Lambda.getSize();

      for (int i = 0; i < dim; i++)
	Lb[i] += Lambda[i][bOne] * factor;
    }

    /** \brief
     * Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in 
     * each component.
     */
    inline void l1(const DimVec<WorldVector<double> >& Lambda,
		   DimVec<double>& Lb,
		   double factor) const
    {
      const int dim = Lambda.getSize();

      for (int i = 0; i < dim; i++) {
	double val = 0.0;
      
	for (int j = 0; j < dimOfWorld; j++)
	  val += Lambda[i][j];

	Lb[i] += val * factor;
      }    
    }

  protected:
    /// Stores the properties of this OperatorTerm
    Flag properties;

    /// Polynomial degree of the term. Used to detemine the degree of the quadrature.
    int degree;

    /// Stores the dimension of the world.
    int dimOfWorld;

    /// List off all fe spaces, the operator term makes use off.
    std::set<const FiniteElemSpace*> auxFeSpaces;

    /// Pointer to the Operator this OperatorTerm belongs to.
    Operator* operat;

    /** \brief
     * In many cases, the vector b in the evaluation \f$ \Lambda \cdot b\f$ has zeros
     * in all components expect one that is set to one. Using the function \ref lb is
     * then unnecessary time consuming. Instead, this variable defines the component
     * of the vector b to be one. The function \ref lb_one is used if this variable is
     * not -1.
     */
    int bOne;

    /// Flag for piecewise constant terms
    static const Flag PW_CONST;

    /// Flag for symmetric terms
    static const Flag SYMMETRIC;

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

}

#endif