DOFVector.h 26.1 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
    }
#endif

255
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
256
    ///
257
258
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
259
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
260
    string name;
261

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    ///
263
    ElementVector elementVector;
264

Thomas Witkowski's avatar
Thomas Witkowski committed
265
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
266
    vector<Operator*> operators;
267

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    vector<double*> operatorFactor;
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
272
    vector<double*> operatorEstFactor;
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
    ///
275
276
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
277
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    int nBasFcts;
279

Thomas Witkowski's avatar
Thomas Witkowski committed
280
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
282
283

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
284
    FeSpaceDofMap *dofMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
285
#endif
286

Thomas Witkowski's avatar
Thomas Witkowski committed
287
  };
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
  /// Specifies which operation should be done after coarsening
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
  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)
321
      {}
322
323
324

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
325
      {}
326
327
328
329
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
330
331
332
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
333

Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
      DOFVector<T> *create() 
      { 
336
	return new DOFVector<T>(feSpace, ""); 
337
      }
338

Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
      void free(DOFVector<T> *vec) 
      { 
341
	delete vec; 
342
      }
343
344
345
346
347
348

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
349
    /// Empty constructor. No initialization!
350
351
352
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
353
    {}
354

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

Thomas Witkowski's avatar
Thomas Witkowski committed
358
    /// Initialization.
Thomas Witkowski's avatar
Thomas Witkowski committed
359
    void init(const FiniteElemSpace* f, string n);
360

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    /// Copy Constructor
Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
    DOFVector(const DOFVector& rhs) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
364
      *this = rhs;   
365
366
367
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
368
    }
369

Thomas Witkowski's avatar
Thomas Witkowski committed
370
    /// Destructor
371
372
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
373
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
374
    typename vector<T>::iterator begin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
375
    { 
376
      return vec.begin(); 
377
    }
378

Thomas Witkowski's avatar
Thomas Witkowski committed
379
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
380
    typename vector<T>::iterator end() 
Thomas Witkowski's avatar
Thomas Witkowski committed
381
    { 
382
      return vec.end(); 
383
    }
384
  
Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
    /// Used by DOFAdmin to compress this DOFVector. Implementation of
    /// \ref DOFIndexedBase::compress()
387
    virtual void compressDOFIndexed(int first, int last,
Thomas Witkowski's avatar
Thomas Witkowski committed
388
				    vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
393
      coarsenOperation = op; 
394
    }
395

Thomas Witkowski's avatar
Thomas Witkowski committed
396
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
397
398
    inline CoarsenOperation getCoarsenOperation() 
    { 
399
      return coarsenOperation; 
400
    }
401

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

405
    /// Interpolation after refinement. Implemented for DOFVector<double>
406
    inline void refineInterpol(RCNeighbourList&, int) {}
407

Thomas Witkowski's avatar
Thomas Witkowski committed
408
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
409
    vector<T>& getVector() 
Thomas Witkowski's avatar
Thomas Witkowski committed
410
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
411
      return vec;
412
    }
413

Thomas Witkowski's avatar
Thomas Witkowski committed
414
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
    inline int getSize() const 
    { 
417
      return vec.size();
418
    } 
419

Thomas Witkowski's avatar
Thomas Witkowski committed
420
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
    inline int getUsedSize() const 
    { 
423
      return this->feSpace->getAdmin()->getUsedSize(); 
424
    }
425

Thomas Witkowski's avatar
Thomas Witkowski committed
426
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
427
428
    inline void resize(int n) 
    { 
429
      FUNCNAME("DOFVector<T>::resize()");
430
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
431
      vec.resize(n);
432
    } 
433

Thomas Witkowski's avatar
Thomas Witkowski committed
434
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
435
436
    inline void resize(int n, T init) 
    { 
437
      FUNCNAME("DOFVector<T>::resize()");
438
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
439
      vec.resize(n, init);
440
    } 
441

Thomas Witkowski's avatar
Thomas Witkowski committed
442
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
443
444
    inline const T& operator[](DegreeOfFreedom i) const 
    {
445
      FUNCNAME("DOFVector<T>::operator[]");
446
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
447
	("Illegal vector index %d.\n", i);
448
      return vec[i];
449
    } 
450

Thomas Witkowski's avatar
Thomas Witkowski committed
451
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
    inline T& operator[](DegreeOfFreedom i) 
    {
454
      FUNCNAME("DOFVector<T>::operator[]");
455
456

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

459
      return vec[i];
460
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
461
 
Thomas Witkowski's avatar
Thomas Witkowski committed
462
    /// Calculates Integral of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
    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;

480

Thomas Witkowski's avatar
Thomas Witkowski committed
481
482
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
483
484
485
486
487
488
    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
489

Thomas Witkowski's avatar
Thomas Witkowski committed
490
491
492
    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
493
494
495
496
497
498
499
    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
500
    /// Calculates L1 norm of this DOFVector
501
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
502
 
Thomas Witkowski's avatar
Thomas Witkowski committed
503
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
    inline double L2Norm(Quadrature* q = NULL) const 
    {
506
      return sqrt(L2NormSquare());
507
    }
508

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

Thomas Witkowski's avatar
Thomas Witkowski committed
512
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
513
514
    inline double H1Norm(Quadrature* q = NULL) const 
    {
515
      return sqrt(H1NormSquare());
516
    }
517

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

Thomas Witkowski's avatar
Thomas Witkowski committed
521
    /// Calculates euclidian norm of this DOFVector
522
523
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
524
    /// Returns square of the euclidian norm.
525
526
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
527
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
528
529
    inline double l2norm() const 
    { 
530
      return nrm2();
531
    }
532

Thomas Witkowski's avatar
Thomas Witkowski committed
533
    /// Calculates the absolute sum of this DOFVector
534
535
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
536
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
537
538
    inline double l1norm() const 
    { 
539
      return asum();
540
    } 
541

Thomas Witkowski's avatar
Thomas Witkowski committed
542
    /// Calculates doublewell of this DOFVector
543
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
544
 
Thomas Witkowski's avatar
Thomas Witkowski committed
545
    /// Calculates the sum of this DOFVector
546
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
547
 
Thomas Witkowski's avatar
Thomas Witkowski committed
548
    /// Sets \ref vec[i] = val, i=0 , ... , size
549
550
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
551
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
552
553
    inline DOFVector<T>& operator=(T d) 
    {
554
555
      set(d); 
      return *this;
556
    } 
557

Thomas Witkowski's avatar
Thomas Witkowski committed
558
    /// Assignment operator between two vectors
559
    DOFVector<T>& operator=(const DOFVector<T>&); 
560

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

564
    /// Returns minimum of DOFVector.
565
566
    T min() const; 

567
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
568
569
    T max() const;

570
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
571
    T absMax() const;
572

573
574
575
576
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
577
578
    void print() const; 

579
    ///
580
581
    int calcMemoryUsage() const;

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

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

588

Thomas Witkowski's avatar
Thomas Witkowski committed
589
590
591
592
593
594
    /// 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 
595
596
597
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
Praetorius, Simon's avatar
Praetorius, Simon committed
598
      return *value;
599
    }
600

Thomas Witkowski's avatar
Thomas Witkowski committed
601
602
603
604
605
606
607
    /// 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;
608

609
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
610
    void serialize(ostream &out) 
Thomas Witkowski's avatar
Thomas Witkowski committed
611
    {
612
613
614
      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));
615
    }
616

617
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
618
    void deserialize(istream &in) 
Thomas Witkowski's avatar
Thomas Witkowski committed
619
    {
620
621
622
623
      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));
624
    }
625

626
627
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
628
629
630

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

631
632
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
633
634

  protected: 
635

636
    /// Data container
Thomas Witkowski's avatar
Thomas Witkowski committed
637
    vector<T> vec; 
638
 
639
    /// Specifies what operation should be performed after coarsening
640
641
642
    CoarsenOperation coarsenOperation;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
643

644
  template<>
645
646
647
648
649
650
651
652
  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
653
654
  const double& DOFVector<double>::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
655
656

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

660
  template<>
661
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
662
663

  template<>
664
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
665

Thomas Witkowski's avatar
Thomas Witkowski committed
666
667
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
668
    return v.min();
669
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
670

Thomas Witkowski's avatar
Thomas Witkowski committed
671
672
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
673
    return v.max();
674
  }
675
676
677
678
679
680
681
682
683
684


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Thomas Witkowski's avatar
Thomas Witkowski committed
685
686
687
    /// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
    /// as DOFContainer at DOFAdmin
    DOFVectorDOF(const FiniteElemSpace* feSpace_, string name_)
688
689
690
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
691
    }
692
  
693
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
694
695
    ~DOFVectorDOF() 
    {
696
      feSpace->getAdmin()->removeDOFContainer(this);
697
    }
698

Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
701
702
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
703
      return DOFVector<DegreeOfFreedom>::operator[](i);
704
    }
705

Thomas Witkowski's avatar
Thomas Witkowski committed
706
707
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
708
709
710
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

711
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
712
713
    int getSize() const 
    {
714
      return DOFVector<DegreeOfFreedom>::getSize();
715
    }
716

717
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
    void resize(int size) 
    {
720
721
722
723
724
725
726
727
728
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
729

730
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
731
732
  double norm(DOFVector<T> *vec) 
  {
733
    return vec->nrm2();
734
  }
735
736

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
737
738
  double L2Norm(DOFVector<T> *vec) 
  {
739
    return vec->L2Norm();
740
  }
741
742

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
749
  void print(DOFVector<T> *vec) 
750
  {
751
    vec->print();
752
  }
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784

  // 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>
785
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
786
787
788
789
790
791
792
793
794
795
796
797
798

  // 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
799
800
  inline void scal(T a, DOFVector<T>& y) 
  {
801
802
    y *= a;
  }
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818

  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);

819
820
821
822
823
824
825
826
827
  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);


828