DOFVector.h 26.3 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

/** \file DOFVector.h */

#ifndef AMDIS_DOFVECTOR_H 
#define AMDIS_DOFVECTOR_H 
 
26
27
28
29
#include <vector> 
#include <memory> 
#include <list> 

30
#include "AMDiS_fwd.h"
31
32
33
34
35
36
37
38
39
40
41
#include "FixVec.h"
#include "Global.h" 
#include "Flag.h" 
#include "RCNeighbourList.h" 
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "DOFContainer.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
Thomas Witkowski's avatar
Thomas Witkowski committed
42
#include "BasisFunction.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
43
#include "FiniteElemSpace.h"
44
#include "SurfaceQuadrature.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
45
46
47
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/ParallelDofMapping.h"
#endif
48

49
50
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
51
52
  using namespace std;

53
54
55
56
  template<typename T> 
  class DOFVectorBase : public DOFIndexed<T>
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
57

58
59
    DOFVectorBase() 
      : feSpace(NULL),
60
	elementVector(3),
Thomas Witkowski's avatar
Thomas Witkowski committed
61
62
        boundaryManager(NULL),
        nBasFcts(0)
63
64
    {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
65
      dofMap = NULL;
66
67
#endif
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
68
    
Thomas Witkowski's avatar
Thomas Witkowski committed
69
    DOFVectorBase(const FiniteElemSpace *f, string n);
70

Thomas Witkowski's avatar
Thomas Witkowski committed
71
    virtual ~DOFVectorBase();
72

Thomas Witkowski's avatar
Thomas Witkowski committed
73
74
    /// For the given element, this function returns an array of all DOFs of 
    /// this DOFVector that are defined on this element.
75
76
    virtual void getLocalVector(const Element *el, 
				mtl::dense_vector<T>& localVec) const;
77

Thomas Witkowski's avatar
Thomas Witkowski committed
78
79
    /// Evaluates the DOF vector at a set of quadrature points defined on the 
    /// given element.
80
81
82
83
    void getVecAtQPs(const ElInfo *elInfo, 
		     const Quadrature *quad,
		     const FastQuadrature *quadFast,
		     mtl::dense_vector<T>& vecAtQPs) const;
84

85
86
87
88
89
    void getVecAtQPs(const ElInfo *smallElInfo, 
		     const ElInfo *largeElInfo,
		     const Quadrature *quad,
		     const FastQuadrature *quadFast,
		     mtl::dense_vector<T>& vecAtQPs) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
90

Thomas Witkowski's avatar
Thomas Witkowski committed
91
92
    /// Evaluates the gradient of a DOF vector at a set of quadrature points
    /// defined on the given element.
93
94
95
96
97
98
99
100
101
102
103
    void getGrdAtQPs( const ElInfo *elInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;

    void getGrdAtQPs( const ElInfo *smallElInfo,
		      const ElInfo *largeElInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
104
105
    /// Evaluates the comp'th component of the derivative of a DOF vector at a
    /// set of quadrature points defined on the given element.
106
107
108
109
110
111
112
113
114
115
116
117
118
    void getDerivativeAtQPs( const ElInfo *elInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      int comp,
		      mtl::dense_vector<T> &derivativeAtQPs) const;

    void getDerivativeAtQPs( const ElInfo *smallElInfo,
		      const ElInfo *largeElInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      int comp,
		      mtl::dense_vector<T> &derivativeAtQPs) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
119
120
    /// Evaluates the jacobian of a DOF vector at a set of quadrature points
    ///  defined on the given element.
121
122
123
124
    void getD2AtQPs(const ElInfo *elInfo,
		    const Quadrature *quad,
		    const FastQuadrature *quadFast,
		    mtl::dense_vector<typename D2Type<T>::type> &d2AtQPs) const;
125

Thomas Witkowski's avatar
Thomas Witkowski committed
126
    /// Returns the FE space the DOF vector is defined on.
127
    inline const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
128
    {
129
      return feSpace;
130
    }
131

Thomas Witkowski's avatar
Thomas Witkowski committed
132
133
    /// Assembles the element vector for the given ellement and adds the
    /// element matrix to the current DOF vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
134
135
136
137
138
139
140
141
142
    void assemble(T factor, ElInfo *elInfo,			    
		  const BoundaryType *bound, 
		  Operator *op = NULL);

    void assemble2(T factor, 
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
		   const BoundaryType *bound, 
		   Operator *op = NULL);
143
144

    void addElementVector(T sign,
145
			  const ElementVector& elVec, 
146
			  const BoundaryType *bound,
147
			  ElInfo *elInfo,
148
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
149

Thomas Witkowski's avatar
Thomas Witkowski committed
150
151
152
    /// That function must be called after the matrix assembling has been
    /// finished. This makes it possible to start some cleanup or matrix
    /// data compressing procedures.
Thomas Witkowski's avatar
Thomas Witkowski committed
153
    void finishAssembling();
154
155
156
157
158
159
160
161
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
162
    }
163

Thomas Witkowski's avatar
Thomas Witkowski committed
164
    inline vector<double*>::iterator getOperatorFactorBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
165
    {
166
      return operatorFactor.begin();
167
    }
168

Thomas Witkowski's avatar
Thomas Witkowski committed
169
    inline vector<double*>::iterator getOperatorFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
170
    {
171
      return operatorFactor.end();
172
    }
173

Thomas Witkowski's avatar
Thomas Witkowski committed
174
    inline vector<double*>::iterator getOperatorEstFactorBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
175
    {
176
      return operatorEstFactor.begin();
177
    }
178

Thomas Witkowski's avatar
Thomas Witkowski committed
179
    inline vector<double*>::iterator getOperatorEstFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
180
    {
181
      return operatorEstFactor.end();
182
    }
183
184


Thomas Witkowski's avatar
Thomas Witkowski committed
185
    inline vector<Operator*>::iterator getOperatorsBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
186
    {
187
      return operators.begin();
188
    }
189

Thomas Witkowski's avatar
Thomas Witkowski committed
190
    inline vector<Operator*>::iterator getOperatorsEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
191
    {
192
      return operators.end();
193
    }
194
195
196
197
198
199
200
201
202
203
204

    Flag getAssembleFlag();

    /** \brief
     * Evaluates \f[ u_h(x(\lambda)) = \sum_{i=0}^{m-1} vec[ind[i]] * 
     * \varphi^i(\lambda) \f] where \f$ \varphi^i \f$ is the i-th basis function,
     * \f$ x(\lambda) \f$ are the world coordinates of lambda and
     * \f$ m \f$ is the number of basis functions
     */
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

Thomas Witkowski's avatar
Thomas Witkowski committed
205
    inline vector<Operator*>& getOperators() 
Thomas Witkowski's avatar
Thomas Witkowski committed
206
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
207
      return operators; 
208
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
209

Thomas Witkowski's avatar
Thomas Witkowski committed
210
    inline vector<double*>& getOperatorFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
211
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
212
      return operatorFactor; 
213
    }
214

Thomas Witkowski's avatar
Thomas Witkowski committed
215
    inline vector<double*>& getOperatorEstFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
216
    { 
217
      return operatorEstFactor; 
218
    }
219

Thomas Witkowski's avatar
Thomas Witkowski committed
220
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
221
    inline string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
222
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
223
      return name; 
224
    } 
225

Thomas Witkowski's avatar
Thomas Witkowski committed
226
    inline void setName(string n)
227
    {
228
      name = n;
229
230
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
233
      return boundaryManager; 
234
    }
235

Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
238
      boundaryManager = bm;
239
    }
240

Thomas Witkowski's avatar
Thomas Witkowski committed
241
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
242
    void setDofMap(FeSpaceDofMap& m)
Thomas Witkowski's avatar
Thomas Witkowski committed
243
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
244
      dofMap = &m;
Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
247
248
    }

    inline bool isRankDof(DegreeOfFreedom dof)
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
249
      TEST_EXIT_DBG(dofMap)("No rank dofs set!\n");
250

251
      return dofMap->isRankDof(dof);
Thomas Witkowski's avatar
Thomas Witkowski committed
252
    }
253
254
255
256
257
258
259
260
261
262
263

    inline void setDirichletDofValue(DegreeOfFreedom dof,
				     T value)
    {
      dirichletDofValues[dof] = value;
    }

    map<DegreeOfFreedom, T>& getDirichletValues()
    {
      return dirichletDofValues;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
#endif

266
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
267
    ///
268
269
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
270
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
271
    string name;
272

Thomas Witkowski's avatar
Thomas Witkowski committed
273
    ///
274
    ElementVector elementVector;
275

Thomas Witkowski's avatar
Thomas Witkowski committed
276
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
277
    vector<Operator*> operators;
278

Thomas Witkowski's avatar
Thomas Witkowski committed
279
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
280
    vector<double*> operatorFactor;
281

Thomas Witkowski's avatar
Thomas Witkowski committed
282
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
283
    vector<double*> operatorEstFactor;
284

Thomas Witkowski's avatar
Thomas Witkowski committed
285
    ///
286
287
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
288
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
289
    int nBasFcts;
290

Thomas Witkowski's avatar
Thomas Witkowski committed
291
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
292
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
295
    FeSpaceDofMap *dofMap;
296

297
298
    map<DegreeOfFreedom, T> dirichletDofValues;
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
299
  };
300

301

Thomas Witkowski's avatar
Thomas Witkowski committed
302
  /// Specifies which operation should be done after coarsening
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
  typedef enum{
    NO_OPERATION = 0,   
    COARSE_RESTRICT = 1,
    COARSE_INTERPOL = 2 
  } CoarsenOperation;
 

  /** \ingroup DOFAdministration 
   * \brief
   * The DOFs described above are just integers that can be used as indices into 
   * vectors and matrices. During refinement and coarsening of the mesh, the 
   * number of used DOFs, the meaning of one integer index, and even the total 
   * range of DOFs change. To be able to handle these changes automatically for 
   * all vectors, which are indexed by the DOFs, special data structures are 
   * used which contain such vector data. Lists of these structures are kept in 
   * DOFAdmin, so that all vectors in the lists can be resized together with the
   * range of DOFs. During refinement and coarsening of elements, values can be
   * interpolated automatically to new DOFs, and restricted from old DOFs.
   */
  template<typename T> 
  class DOFVector : public DOFVectorBase<T>, public Serializable
  {  
  public:
    /** \ingroup DOFAdministration
     * \brief
     * Enables the access of DOFVector<T>::Iterator. Alias for DOFIterator<T>
     */
    class Iterator : public DOFIterator<T> {
    public:
      Iterator(DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(c, type)
334
      {}
335
336
337

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
338
      {}
339
340
341
342
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
343
344
345
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
346

Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
      DOFVector<T> *create() 
      { 
349
	return new DOFVector<T>(feSpace, ""); 
350
      }
351

Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
      void free(DOFVector<T> *vec) 
      { 
354
	delete vec; 
355
      }
356
357
358
359
360
361

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
362
    /// Empty constructor. No initialization!
363
364
365
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
366
    {}
367

Thomas Witkowski's avatar
Thomas Witkowski committed
368
    /// Constructs a DOFVector with name n belonging to FiniteElemSpace f
Thomas Witkowski's avatar
Thomas Witkowski committed
369
    DOFVector(const FiniteElemSpace* f, string n); 
370

Thomas Witkowski's avatar
Thomas Witkowski committed
371
    /// Initialization.
Thomas Witkowski's avatar
Thomas Witkowski committed
372
    void init(const FiniteElemSpace* f, string n);
373

Thomas Witkowski's avatar
Thomas Witkowski committed
374
    /// Copy Constructor
Thomas Witkowski's avatar
Thomas Witkowski committed
375
376
    DOFVector(const DOFVector& rhs) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
377
      *this = rhs;   
378
379
380
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
381
    }
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
    /// Destructor
384
385
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
386
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
387
    typename vector<T>::iterator begin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
388
    { 
389
      return vec.begin(); 
390
    }
391

Thomas Witkowski's avatar
Thomas Witkowski committed
392
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
393
    typename vector<T>::iterator end() 
Thomas Witkowski's avatar
Thomas Witkowski committed
394
    { 
395
      return vec.end(); 
396
    }
397
  
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
    /// Used by DOFAdmin to compress this DOFVector. Implementation of
    /// \ref DOFIndexedBase::compress()
400
    virtual void compressDOFIndexed(int first, int last,
Thomas Witkowski's avatar
Thomas Witkowski committed
401
				    vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
402
403

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
404
405
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
406
      coarsenOperation = op; 
407
    }
408

Thomas Witkowski's avatar
Thomas Witkowski committed
409
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
410
411
    inline CoarsenOperation getCoarsenOperation() 
    { 
412
      return coarsenOperation; 
413
    }
414

Thomas Witkowski's avatar
Thomas Witkowski committed
415
    /// Restriction after coarsening. Implemented for DOFVector<double>
416
    inline void coarseRestrict(RCNeighbourList&, int) {}
417

418
    /// Interpolation after refinement. Implemented for DOFVector<double>
419
    inline void refineInterpol(RCNeighbourList&, int) {}
420

Thomas Witkowski's avatar
Thomas Witkowski committed
421
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
422
    vector<T>& getVector() 
Thomas Witkowski's avatar
Thomas Witkowski committed
423
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
424
      return vec;
425
    }
426

Thomas Witkowski's avatar
Thomas Witkowski committed
427
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
428
429
    inline int getSize() const 
    { 
430
      return vec.size();
431
    } 
432

Thomas Witkowski's avatar
Thomas Witkowski committed
433
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
434
435
    inline int getUsedSize() const 
    { 
436
      return this->feSpace->getAdmin()->getUsedSize(); 
437
    }
438

Thomas Witkowski's avatar
Thomas Witkowski committed
439
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
    inline void resize(int n) 
    { 
442
      FUNCNAME("DOFVector<T>::resize()");
443
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
444
      vec.resize(n);
445
    } 
446

Thomas Witkowski's avatar
Thomas Witkowski committed
447
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
448
449
    inline void resize(int n, T init) 
    { 
450
      FUNCNAME("DOFVector<T>::resize()");
451
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
452
      vec.resize(n, init);
453
    } 
454

Thomas Witkowski's avatar
Thomas Witkowski committed
455
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
    inline const T& operator[](DegreeOfFreedom i) const 
    {
458
      FUNCNAME("DOFVector<T>::operator[]");
459
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
460
	("Illegal vector index %d.\n", i);
461
      return vec[i];
462
    } 
463

Thomas Witkowski's avatar
Thomas Witkowski committed
464
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
465
466
    inline T& operator[](DegreeOfFreedom i) 
    {
467
      FUNCNAME("DOFVector<T>::operator[]");
468
469

      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size())) 
470
 	("Illegal vector index %d.\n", i); 
471

472
      return vec[i];
473
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
474
 
Thomas Witkowski's avatar
Thomas Witkowski committed
475
    /// Calculates Integral of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
    double Int(Quadrature* q = NULL) const
    {
      return Int(-1, q);
    }

    /** \brief
     * Calculates Integral of this DOFVector. 
     *
     * \param[in]  meshlevel  Then mesh level on which the integral should be
     *                        calculated. Usually, only -1 for the leaf level
     *                        makes sence, but other values can be used, e.g.,
     *                        for debugging.
     * \param[in]  q          Quadrature object. If not specified, the function
     *                        creates a new quadrature object.
     */
    double Int(int meshLevel, Quadrature* q = NULL) const;

493

Thomas Witkowski's avatar
Thomas Witkowski committed
494
495
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
496
497
498
499
500
501
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
502

Thomas Witkowski's avatar
Thomas Witkowski committed
503
504
505
    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
506
507
508
509
510
511
512
    double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = NULL) const
    {
      FUNCNAME("DOFVector::IntOnBoundaryNormal())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
513
    /// Calculates L1 norm of this DOFVector
514
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
515
 
Thomas Witkowski's avatar
Thomas Witkowski committed
516
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
517
518
    inline double L2Norm(Quadrature* q = NULL) const 
    {
519
      return sqrt(L2NormSquare());
520
    }
521

Thomas Witkowski's avatar
Thomas Witkowski committed
522
    /// Calculates square of L2 norm of this DOFVector
523
524
    double L2NormSquare(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
525
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
526
527
    inline double H1Norm(Quadrature* q = NULL) const 
    {
528
      return sqrt(H1NormSquare());
529
    }
530

Thomas Witkowski's avatar
Thomas Witkowski committed
531
    /// Calculates square of H1 norm of this DOFVector
532
533
    double H1NormSquare(Quadrature* q = NULL) const;  

Thomas Witkowski's avatar
Thomas Witkowski committed
534
    /// Calculates euclidian norm of this DOFVector
535
536
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
537
    /// Returns square of the euclidian norm.
538
539
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
540
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
541
542
    inline double l2norm() const 
    { 
543
      return nrm2();
544
    }
545

Thomas Witkowski's avatar
Thomas Witkowski committed
546
    /// Calculates the absolute sum of this DOFVector
547
548
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
549
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
550
551
    inline double l1norm() const 
    { 
552
      return asum();
553
    } 
554

Thomas Witkowski's avatar
Thomas Witkowski committed
555
    /// Calculates doublewell of this DOFVector
556
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
557
 
Thomas Witkowski's avatar
Thomas Witkowski committed
558
    /// Calculates the sum of this DOFVector
559
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
560
 
Thomas Witkowski's avatar
Thomas Witkowski committed
561
    /// Sets \ref vec[i] = val, i=0 , ... , size
562
563
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
564
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
565
566
    inline DOFVector<T>& operator=(T d) 
    {
567
568
      set(d); 
      return *this;
569
    } 
570

Thomas Witkowski's avatar
Thomas Witkowski committed
571
    /// Assignment operator between two vectors
572
    DOFVector<T>& operator=(const DOFVector<T>&); 
573

Thomas Witkowski's avatar
Thomas Witkowski committed
574
    /// vec[i] = v.vec[i]
575
576
    void copy(const DOFVector<T>& v); 

577
    /// Returns minimum of DOFVector.
578
579
    T min() const; 

580
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
581
582
    T max() const;

583
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
584
    T absMax() const;
585

586
587
588
589
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
590
591
    void print() const; 

592
    ///
593
594
    int calcMemoryUsage() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
595
596
    /// Computes the coefficients of the interpolant of the function fct and
    /// stores these in the DOFVector
597
598
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

599
    void interpol(DOFVector<T> *v, double factor = 1.0);
600

601

Thomas Witkowski's avatar
Thomas Witkowski committed
602
603
604
605
606
607
    /// Eval DOFVector at given point p. If oldElInfo != NULL the search for 
    /// the element, where p is inside, starts from oldElInfo. implemented for:
    /// double, WorldVector< double >
    inline const T& evalAtPoint(WorldVector<double> &p, 
				ElInfo *oldElInfo = NULL, 
				T* value = NULL) const 
608
609
610
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
Praetorius, Simon's avatar
Praetorius, Simon committed
611
      return *value;
612
    }
613

Thomas Witkowski's avatar
Thomas Witkowski committed
614
615
616
617
618
619
620
    /// Determine the DegreeOfFreedom that has coords with minimal euclidean 
    /// distance to WorldVector p. return true if DOF is found, and false 
    /// otherwise.
    const bool getDofIdxAtPoint(WorldVector<double> &p, 
				DegreeOfFreedom &idx, 
				ElInfo *oldElInfo = NULL, 
				bool useOldElInfo = false) const;
621

622
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
623
    void serialize(ostream &out) 
Thomas Witkowski's avatar
Thomas Witkowski committed
624
    {
625
626
627
      unsigned int size = vec.size();
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
628
    }
629

630
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
631
    void deserialize(istream &in) 
Thomas Witkowski's avatar
Thomas Witkowski committed
632
    {
633
634
635
636
      unsigned int size;
      in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int));
      vec.resize(size);
      in.read(reinterpret_cast<char*>(&(vec[0])), size * sizeof(T));
637
    }
638

639
640
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
641
642
643

    WorldVector<DOFVector<T>*> *getGradient(WorldVector<DOFVector<T>*> *grad) const;

644
645
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
646
647

  protected: 
648

649
    /// Data container
Thomas Witkowski's avatar
Thomas Witkowski committed
650
    vector<T> vec; 
651
 
652
    /// Specifies what operation should be performed after coarsening
653
654
655
    CoarsenOperation coarsenOperation;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
656

657
  template<>
658
659
660
661
662
663
664
665
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

  template<>
  double DOFVector<WorldVector<double> >::IntOnBoundaryNormal(
    BoundaryType boundaryType, Quadrature* q) const;

  template<>
Praetorius, Simon's avatar
Praetorius, Simon committed
666
667
  const double& DOFVector<double>::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
668
669

  template<>
Praetorius, Simon's avatar
Praetorius, Simon committed
670
671
  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
672

673
  template<>
674
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
675
676

  template<>
677
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
678

Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
681
    return v.min();
682
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
683

Thomas Witkowski's avatar
Thomas Witkowski committed
684
685
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
686
    return v.max();
687
  }
688
689
690
691
692
693
694
695
696
697


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Thomas Witkowski's avatar
Thomas Witkowski committed
698
699
700
    /// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
    /// as DOFContainer at DOFAdmin
    DOFVectorDOF(const FiniteElemSpace* feSpace_, string name_)
701
702
703
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
704
    }
705
  
706
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
707
708
    ~DOFVectorDOF() 
    {
709
      feSpace->getAdmin()->removeDOFContainer(this);
710
    }
711

Thomas Witkowski's avatar
Thomas Witkowski committed
712
713
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
716
      return DOFVector<DegreeOfFreedom>::operator[](i);
717
    }
718

Thomas Witkowski's avatar
Thomas Witkowski committed
719
720
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
721
722
723
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

724
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
725
726
    int getSize() const 
    {
727
      return DOFVector<DegreeOfFreedom>::getSize();
728
    }
729

730
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
731
732
    void resize(int size) 
    {
733
734
735
736
737
738
739
740
741
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
742

743
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
  double norm(DOFVector<T> *vec) 
  {
746
    return vec->nrm2();
747
  }
748
749

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
  double L2Norm(DOFVector<T> *vec) 
  {
752
    return vec->L2Norm();
753
  }
754
755

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
756
757
  double H1Norm(DOFVector<T> *vec) 
  {
758
    return vec->H1Norm();
759
  }
760
761

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
762
  void print(DOFVector<T> *vec) 
763
  {
764
    vec->print();
765
  }
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797

  // point wise multiplication
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y);

  // multiplication with scalar
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal);

  // scalar product
  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y);

  // addition
  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y);

  // subtraction
  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y);

  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d);

  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v);

  template<typename T>
  const DOFVector<T>& operator+(const DOFVector<T>&v1 , const DOFVector<T>& v2);

  // y = a*x + y
  template<typename T>
798
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
799
800
801
802
803
804
805
806
807
808
809
810
811

  // matrix vector product
  template<typename T>
  void mv(MatrixTranspose transpose, 
	  const DOFMatrix &a, 
	  const DOFVector<T> &x,
	  DOFVector<T> &result,
	  bool add = false); 

  template<typename T>
  void xpay(double a,const DOFVector<T>& x,DOFVector<T>& y);

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
812
813
  inline void scal(T a, DOFVector<T>& y) 
  {
814
815
    y *= a;
  }
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831

  template<typename T>
  inline const DOFVector<T>& mult(double scal,
				  const DOFVector<T>& v,
				  DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v1,
				 const DOFVector<T>& v2,
				 DOFVector<T>& result);

832
833
834
835
836
837
838
839
840
  template<typename T>
  inline const DOFVector<T>& mod(const DOFVector<T>& v,
				 DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& Tanh(const DOFVector<T>& v,
				  DOFVector<T>& result);


841
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
842
843
  inline void set(DOFVector<T>& vec, T d) 
  {
844
    vec.set(d);
845
  }
846
847

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
848
849
  inline void setValue(DOFVector<T>& vec, T d) 
  {
850
    vec.set(d);
851
  }
852
853

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
854
855
  inline int size(DOFVector<T> *vec) 
  {
856
    return vec->getUsedSize();
857
  }
858

859
860
861
862
863
864
  template<typename T>
  inline int size(const DOFVector<T>& vec) 
  {
    return vec.getUsedSize();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
865
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
866
  inline void checkFeSpace(const FiniteElemSpace* feSpace, const vector<T>& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
867
868
869
870
871
872
873
874
  {
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       feSpace->getAdmin()->getUsedSize());
  }

875
876
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
877
878
  
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
879
880
  vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     vector<DOFVector<double>*> *res);
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905

  
  /// Computes the integral: \f$ \int f(\vec{x})\,\text{d}\vec{x}\f$
  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct);
  
  template<typename TOut>
  TOut integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
  {
    return integrate_Coords(feSpace, fct);
  }
  
  /// Computes the integral: \f$ \int f(v(\vec{x}))\,\text{d}\vec{x}\f$
  template<typename TOut, typename T>
  TOut integrate_Vec(const DOFVector<T> &vec,
		  AbstractFunction<TOut, T> *fct);
  
  template<typename TOut, typename T>
  TOut integrate(const DOFVector<T> &vec,
		  AbstractFunction<TOut, T> *fct)
  {
    return integrate_Vec(vec, fct);
  }
906
  
907
  /** \brief
908
909
   * Computes the integral of a function that includes two different DOFVectors:
   * \f$ \int f(v(\vec{x}), w(\vec{x}))\,\text{d}\vec{x}\f$
Thomas Witkowski's avatar
Thomas Witkowski committed
910
911
   * This function works also for the case that the DOFVectors are defined on 
   * two different meshes.
912
   */
913
914
915
916
917
918
919
920
921
922
923
924
925
  template<typename TOut, typename T1, typename T2>
  TOut integrate_Vec2(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct);

  template<typename TOut, typename T1, typename T2>
  TOut integrate(const DOFVector<T1> &vec1,
		   const DOFVector<T2