SubAssembler.h 7.57 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
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ==  http://www.amdis-fem.org                                              ==
// ==                                                                        ==
// ============================================================================
//
// 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.



/** \file SubAssembler.h */

23
24
25
26
27
#ifndef AMDIS_SUB_ASSEMBLER_H
#define AMDIS_SUB_ASSEMBLER_H

#include <map>
#include <vector>
28
#include "AMDiS_fwd.h"
29
30
#include "FixVec.h"
#include "Flag.h"
31
#include <boost/any.hpp>
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

namespace AMDiS {

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Base class for SecondOrderAssembler, FirstOrderAssembler, 
   * ZeroOrderAssembler. The task of a SubAssembler is to assemble a list of 
   * terms of a special order and add their contributions to a DOFMatrix or a 
   * DOFVector. An Assembler can consist of up to four SubAssemblers: one 
   * SecondOrderAssembler for second order terms, one ZeroOrderAssembler for 
   * terms of order zero, and two FirstOrderAssemblers. One for terms with
   * derivatives of the basis functions connected to to row DOFs and one for 
   * those connected to the column DOFs.
   */
  class SubAssembler
  {
  public:
    /** \brief
     * Creates a SubAssembler belonging to assembler for the terms of given 
     * order of Operator op. If order is equal to one, type spezifies what kind 
     * of FirstOrderType are to assemble. During construction of a SubAssembler
     * the needs and properties of the terms are considered.
     */
    SubAssembler(Operator *op,
		 Assembler *assembler,
		 Quadrature *quadrat,
		 int order, 
		 bool optimized,
		 FirstOrderType type = GRD_PHI);

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

67
68
    /// Calculates the element matrix for elInfo and adds it to mat. Memory for
    /// mat must be provided by the caller.
69
    virtual void calculateElementMatrix(const ElInfo *elInfo,
70
					ElementMatrix& mat) = 0;
71

72
73
    /// Calculates the element vector for elInfo and adds it to vec. Memory for
    /// vec must be provided by the caller.
74
    virtual void calculateElementVector(const ElInfo *elInfo,
75
					ElementVector& vec) = 0;
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
    /// Returns \ref terms
78
79
    inline std::vector<OperatorTerm*> *getTerms() 
    { 
80
      return &terms; 
Thomas Witkowski's avatar
Thomas Witkowski committed
81
    }
82

Thomas Witkowski's avatar
Thomas Witkowski committed
83
    /// Returns \ref quadrature.
84
85
    inline Quadrature *getQuadrature() 
    {
86
      return quadrature;
Thomas Witkowski's avatar
Thomas Witkowski committed
87
    }
88

Thomas Witkowski's avatar
Thomas Witkowski committed
89
    /// Sets \ref quadrature to q.
90
91
    inline void setQuadrature(Quadrature* q) 
    {
92
      quadrature = q;
Thomas Witkowski's avatar
Thomas Witkowski committed
93
    }
94
95
96
97
98
99
100
101
102
  
    /** \brief
     * Returns a vector with the world coordinates of the quadrature points
     * of \ref quadrature on the element of elInfo. 
     * Used by \ref OperatorTerm::initElement().
     */
    WorldVector<double>* getCoordsAtQPs(const ElInfo* elInfo, 
					Quadrature *quad = NULL);

103
104
    /// DOFVector dv evaluated at quadrature points.
    /// Used by \ref OperatorTerm::initElement().
105
    template<typename T>
106
107
108
109
    void getVectorAtQPs(DOFVectorBase<T>* dv, 
			const ElInfo* elInfo,
			Quadrature *quad,
			mtl::dense_vector<T>& vecAtQPs);
Thomas Witkowski's avatar
Thomas Witkowski committed
110

Thomas Witkowski's avatar
Thomas Witkowski committed
111
    ///
112
    template<typename T>
113
114
115
116
117
    void getVectorAtQPs(DOFVectorBase<T>* dv, 
			const ElInfo* smallElInfo,
			const ElInfo* largeElInfo,
			Quadrature *quad,
			mtl::dense_vector<T>& vecAtQPs);
118
    
119
120
    /// Gradients of DOFVector dv evaluated at quadrature points.
    /// Used by \ref OperatorTerm::initElement().
121
122
123
124
125
    template<typename T>
    void getGradientsAtQPs( DOFVectorBase<T>* dv,
			    const ElInfo* elInfo,
			    Quadrature *quad,
			    mtl::dense_vector<typename GradientType<T>::type>& grdAtQPs);
126

127
128
129
130
131
132
133
    template<typename T>
    void getGradientsAtQPs( DOFVectorBase<T>* dv,
			    const ElInfo* smallElInfo,
			    const ElInfo* largeElInfo,
			    Quadrature *quad,
			    mtl::dense_vector<typename GradientType<T>::type>& grdAtQPs);

134
135
136
137
    /// The comp'th component of the derivative of DOFVector dv evaluated at
    /// quadrature points. Used by \ref OperatorTerm::initElement().
    /// Attention: not caching at the moment! Using cache if gradients for read 
    /// but not for write.
138
139
140
141
142
143
144
145
146
147
148
149
150
151
    template<typename T>
    void getDerivativeAtQPs(DOFVectorBase<T>* dv,
			    const ElInfo* elInfo,
			    Quadrature *quad,
			    int comp,
			    mtl::dense_vector<T>& grdAtQPs);

    template<typename T>
    void getDerivativeAtQPs(DOFVectorBase<T>* dv,
			    const ElInfo* smallElInfo,
			    const ElInfo* largeElInfo,
			    Quadrature *quad,
			    int comp,
			    mtl::dense_vector<T>& grdAtQPs);
Thomas Witkowski's avatar
Thomas Witkowski committed
152

153
154
155
    /// Called once for each ElInfo when \ref calculateElementMatrix() or
    /// \ref calculateElementVector() is called for the first time for this
    /// Element.
Thomas Witkowski's avatar
Thomas Witkowski committed
156
157
    virtual void initElement(const ElInfo *smallElInfo,
			     const ElInfo *largeElInfo = NULL,
158
159
			     Quadrature *quad = NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
160
    /// Returns \ref psiFast.
161
162
    const FastQuadrature *getPsiFast() const 
    { 
163
      return psiFast; 
Thomas Witkowski's avatar
Thomas Witkowski committed
164
    }
165

Thomas Witkowski's avatar
Thomas Witkowski committed
166
    // Returns \ref phiFast.
167
168
    const FastQuadrature *getPhiFast() const 
    { 
169
      return phiFast; 
Thomas Witkowski's avatar
Thomas Witkowski committed
170
    }
171

172
173
174
175
176
177
    /// Returns \ref name.
    std::string getName() const
    {
      return name;
    }

178
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
179
    /// Updates \ref psiFast and \ref phiFast.
180
181
182
183
184
    FastQuadrature *updateFastQuadrature(FastQuadrature *quadFast,
					 const BasisFunction *psi,
					 Flag updateFlag);
  
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
185
    /// Problem dimension
186
187
    int dim;

188
189
190
191
192
193
    /// Row FiniteElemSpace.
    const FiniteElemSpace *rowFeSpace;

    /// Column FiniteElemSpace.
    const FiniteElemSpace *colFeSpace;

194
195
    /// Number of rows of the element matrix and length of the element
    /// vector. Is equal to the number of row basis functions
196
197
    int nRow;

198
199
    /// Number of columns of the element matrix. Is equal to the number
    /// of column basis functions
200
201
    int nCol;

202
203
    /// Used for \ref getVectorAtQPs() and \ref getGradientsAtQPs().
    struct ValuesAtQPs {
Thomas Witkowski's avatar
Thomas Witkowski committed
204
      ValuesAtQPs()
205
	: valid(false)
Thomas Witkowski's avatar
Thomas Witkowski committed
206
      {}
207
208
      ValuesAtQPs(boost::any values_, bool valid_)
	: values(values_), valid(valid_)
Thomas Witkowski's avatar
Thomas Witkowski committed
209
      {}
210
211
      
      boost::any values; // used for mtl::dense_vector<T>
212
213
214
      bool valid;
    };

215
216
217
    std::map<const void*, ValuesAtQPs* > cachedValuesAtQPs;
    std::map<const void*, ValuesAtQPs* > cachedGradientsAtQPs;
    
218
219
    /// Set and updated by \ref initElement() for each ElInfo. 
    /// coordsAtQPs[i] points to the coordinates of the i-th quadrature point.
220
221
    WorldVector<double> *coordsAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
222
    /// Used for \ref getCoordsAtQPs().
223
224
    bool coordsValid;

225
226
    /// Used for \ref getCoordsAtQP(). Stores the number of allocated 
    /// WorldVectors.
227
228
    int coordsNumAllocated;

Thomas Witkowski's avatar
Thomas Witkowski committed
229
    /// Quadrature object to be used for assembling.
230
231
    Quadrature *quadrature;

Thomas Witkowski's avatar
Thomas Witkowski committed
232
    /// FastQuadrature for row basis functions
233
234
    FastQuadrature *psiFast;

Thomas Witkowski's avatar
Thomas Witkowski committed
235
    /// FastQuadrature for column basis functions
236
237
    FastQuadrature *phiFast;

Thomas Witkowski's avatar
Thomas Witkowski committed
238
    /// Flag that specifies whether the element matrix is symmetric.
239
240
    bool symmetric;

Thomas Witkowski's avatar
Thomas Witkowski committed
241
    /// List of all terms with a contribution to this SubAssembler
242
    std::vector<OperatorTerm*> terms;
243

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    ///
245
246
    bool opt;

Thomas Witkowski's avatar
Thomas Witkowski committed
247
    ///
248
249
    bool firstCall;

250
251
252
    /// Name of the assembler. Is used to print information about used assembler.
    std::string name;

253
254
255
256
257
    friend class Assembler;
  };

}

258
259
#include "SubAssembler.hh"

260
#endif // AMDSIS_SUB_ASSEMBLER_H