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

    Flag getAssembleFlag();

197
198
199
200
    /// 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
201
202
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

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

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

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

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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    string name;
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    ///
272
    ElementVector elementVector;
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
275
    vector<Operator*> operators;
276

Thomas Witkowski's avatar
Thomas Witkowski committed
277
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    vector<double*> operatorFactor;
279

Thomas Witkowski's avatar
Thomas Witkowski committed
280
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    vector<double*> operatorEstFactor;
282

Thomas Witkowski's avatar
Thomas Witkowski committed
283
    ///
284
285
    BoundaryManager *boundaryManager;

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

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

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
293
    FeSpaceDofMap *dofMap;
294

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

299

Thomas Witkowski's avatar
Thomas Witkowski committed
300
  /// Specifies which operation should be done after coarsening
301
302
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
  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)
332
      {}
333
334
335

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

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

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
381
    /// Destructor
382
383
    virtual ~DOFVector();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

470
      return vec[i];
471
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
472
 
Thomas Witkowski's avatar
Thomas Witkowski committed
473
    /// Calculates Integral of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
    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;

491

Thomas Witkowski's avatar
Thomas Witkowski committed
492
493
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
494
495
496
497
498
499
    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
500

Thomas Witkowski's avatar
Thomas Witkowski committed
501
502
503
    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
504
505
506
507
508
509
510
    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
511
    /// Calculates L1 norm of this DOFVector
512
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
513
 
Thomas Witkowski's avatar
Thomas Witkowski committed
514
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
515
516
    inline double L2Norm(Quadrature* q = NULL) const 
    {
517
      return sqrt(L2NormSquare());
518
    }
519

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

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

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

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

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

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

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

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

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

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

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

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

575
    /// Returns minimum of DOFVector.
576
577
    T min() const; 

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

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

584
585
586
587
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
588
589
    void print() const; 

590
    ///
591
592
    int calcMemoryUsage() const;

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

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

599

Thomas Witkowski's avatar
Thomas Witkowski committed
600
601
602
603
604
605
    /// 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 
606
607
608
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
Praetorius, Simon's avatar
Praetorius, Simon committed
609
      return *value;
610
    }
611

Thomas Witkowski's avatar
Thomas Witkowski committed
612
613
614
615
616
617
618
    /// 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;
619

620
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
621
    void serialize(ostream &out) 
Thomas Witkowski's avatar
Thomas Witkowski committed
622
    {
623
624
625
      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));
626
    }
627

628
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
629
    void deserialize(istream &in) 
Thomas Witkowski's avatar
Thomas Witkowski committed
630
    {
631
632
633
634
      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));
635
    }
636

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

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

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

  protected: 
646

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

Thomas Witkowski's avatar
Thomas Witkowski committed
654

655
  template<>
656
657
658
659
660
661
662
663
  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
664
665
  const double& DOFVector<double>::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
666
667

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

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

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

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

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


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

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
740

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

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
760
  void print(DOFVector<T> *vec) 
761
  {
762
    vec->print();
763
  }
764
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

  // 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>
796
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
797
798
799
800
801
802
803
804
805
806
807
808
809

  // 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
810
811
  inline void scal(T a, DOFVector<T>& y) 
  {
812
813
    y *= a;
  }
814
815
816
817
818
819
820