DOFVector.h 26.6 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
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
395
396
397
398
399
400
401
402
403
404
405
406
407

    /// Returns const_iterator to the begin of \ref vec
    typename vector<T>::const_iterator begin() const
    {
      return vec.begin();
    }

    /// Returns const_iterator to the end of \ref vec
    typename vector<T>::const_iterator end() const
    {
      return vec.end();
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
    /// Used by DOFAdmin to compress this DOFVector. Implementation of
    /// \ref DOFIndexedBase::compress()
410
    virtual void compressDOFIndexed(int first, int last,
Thomas Witkowski's avatar
Thomas Witkowski committed
411
				    vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
412
413

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
414
415
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
416
      coarsenOperation = op; 
417
    }
418

Thomas Witkowski's avatar
Thomas Witkowski committed
419
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
420
421
    inline CoarsenOperation getCoarsenOperation() 
    { 
422
      return coarsenOperation; 
423
    }
424

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

428
    /// Interpolation after refinement. Implemented for DOFVector<double>
429
    inline void refineInterpol(RCNeighbourList&, int) {}
430

Thomas Witkowski's avatar
Thomas Witkowski committed
431
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
432
    vector<T>& getVector() 
Thomas Witkowski's avatar
Thomas Witkowski committed
433
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
434
      return vec;
435
    }
436

Thomas Witkowski's avatar
Thomas Witkowski committed
437
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
438
439
    inline int getSize() const 
    { 
440
      return vec.size();
441
    } 
442

Thomas Witkowski's avatar
Thomas Witkowski committed
443
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
    inline int getUsedSize() const 
    { 
446
      return this->feSpace->getAdmin()->getUsedSize(); 
447
    }
448

Thomas Witkowski's avatar
Thomas Witkowski committed
449
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
450
451
    inline void resize(int n) 
    { 
452
      FUNCNAME("DOFVector<T>::resize()");
453
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
454
      vec.resize(n);
455
    } 
456

Thomas Witkowski's avatar
Thomas Witkowski committed
457
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
    inline void resize(int n, T init) 
    { 
460
      FUNCNAME("DOFVector<T>::resize()");
461
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
462
      vec.resize(n, init);
463
    } 
464

Thomas Witkowski's avatar
Thomas Witkowski committed
465
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
466
467
    inline const T& operator[](DegreeOfFreedom i) const 
    {
468
      FUNCNAME("DOFVector<T>::operator[]");
469
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
470
	("Illegal vector index %d.\n", i);
471
      return vec[i];
472
    } 
473

Thomas Witkowski's avatar
Thomas Witkowski committed
474
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
475
476
    inline T& operator[](DegreeOfFreedom i) 
    {
477
      FUNCNAME("DOFVector<T>::operator[]");
478
479

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

482
      return vec[i];
483
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
484
 
Thomas Witkowski's avatar
Thomas Witkowski committed
485
    /// Calculates Integral of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
    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;

503

Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
506
507
508
509
510
511
    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
512

Thomas Witkowski's avatar
Thomas Witkowski committed
513
514
515
    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
516
517
518
519
520
521
522
    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
523
    /// Calculates L1 norm of this DOFVector
524
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
525
 
Thomas Witkowski's avatar
Thomas Witkowski committed
526
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
527
528
    inline double L2Norm(Quadrature* q = NULL) const 
    {
529
      return sqrt(L2NormSquare());
530
    }
531

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

Thomas Witkowski's avatar
Thomas Witkowski committed
535
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
536
537
    inline double H1Norm(Quadrature* q = NULL) const 
    {
538
      return sqrt(H1NormSquare());
539
    }
540

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

Thomas Witkowski's avatar
Thomas Witkowski committed
544
    /// Calculates euclidian norm of this DOFVector
545
546
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
547
    /// Returns square of the euclidian norm.
548
549
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
550
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
    inline double l2norm() const 
    { 
553
      return nrm2();
554
    }
555

Thomas Witkowski's avatar
Thomas Witkowski committed
556
    /// Calculates the absolute sum of this DOFVector
557
558
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
559
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
560
561
    inline double l1norm() const 
    { 
562
      return asum();
563
    } 
564

Thomas Witkowski's avatar
Thomas Witkowski committed
565
    /// Calculates doublewell of this DOFVector
566
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
567
 
Thomas Witkowski's avatar
Thomas Witkowski committed
568
    /// Calculates the sum of this DOFVector
569
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
570
 
Thomas Witkowski's avatar
Thomas Witkowski committed
571
    /// Sets \ref vec[i] = val, i=0 , ... , size
572
573
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
574
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
575
576
    inline DOFVector<T>& operator=(T d) 
    {
577
578
      set(d); 
      return *this;
579
    } 
580

Thomas Witkowski's avatar
Thomas Witkowski committed
581
    /// Assignment operator between two vectors
582
    DOFVector<T>& operator=(const DOFVector<T>&); 
583

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

587
    /// Returns minimum of DOFVector.
588
589
    T min() const; 

590
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
591
592
    T max() const;

593
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
594
    T absMax() const;
595

596
597
598
599
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
600
601
    void print() const; 

602
    ///
603
604
    int calcMemoryUsage() const;

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

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

611

Thomas Witkowski's avatar
Thomas Witkowski committed
612
613
614
615
616
617
    /// 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 
618
619
620
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
Praetorius, Simon's avatar
Praetorius, Simon committed
621
      return *value;
622
    }
623

Thomas Witkowski's avatar
Thomas Witkowski committed
624
625
626
627
628
629
630
    /// 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;
631

632
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
633
    void serialize(ostream &out) 
Thomas Witkowski's avatar
Thomas Witkowski committed
634
    {
635
636
637
      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));
638
    }
639

640
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
641
    void deserialize(istream &in) 
Thomas Witkowski's avatar
Thomas Witkowski committed
642
    {
643
644
645
646
      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));
647
    }
648

649
650
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
651
652
653

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

654
655
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
656
657

  protected: 
658

659
    /// Data container
Thomas Witkowski's avatar
Thomas Witkowski committed
660
    vector<T> vec; 
661
 
662
    /// Specifies what operation should be performed after coarsening
663
664
665
    CoarsenOperation coarsenOperation;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
666

667
  template<>
668
669
670
671
672
673
674
675
  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
676
677
  const double& DOFVector<double>::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
678
679

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

683
  template<>
684
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
685
686

  template<>
687
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
688

Thomas Witkowski's avatar
Thomas Witkowski committed
689
690
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
691
    return v.min();
692
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
693

Thomas Witkowski's avatar
Thomas Witkowski committed
694
695
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
696
    return v.max();
697
  }
698
699
700
701
702
703
704
705
706
707


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Thomas Witkowski's avatar
Thomas Witkowski committed
708
709
710
    /// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
    /// as DOFContainer at DOFAdmin
    DOFVectorDOF(const FiniteElemSpace* feSpace_, string name_)
711
712
713
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
714
    }
715
  
716
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
717
718
    ~DOFVectorDOF() 
    {
719
      feSpace->getAdmin()->removeDOFContainer(this);
720
    }
721

Thomas Witkowski's avatar
Thomas Witkowski committed
722
723
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
724
725
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
726
      return DOFVector<DegreeOfFreedom>::operator[](i);
727
    }
728

Thomas Witkowski's avatar
Thomas Witkowski committed
729
730
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
731
732
733
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

734
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
    int getSize() const 
    {
737
      return DOFVector<DegreeOfFreedom>::getSize();
738
    }
739

740
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
741
742
    void resize(int size) 
    {
743
744
745
746
747
748
749
750
751
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
752

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
760
761
  double L2Norm(DOFVector<T> *vec) 
  {
762
    return vec->L2Norm();
763
  }
764
765

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
766
767
  double H1Norm(DOFVector<T> *vec) 
  {
768
    return vec->H1Norm();
769
  }
770
771

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
772
  void print(DOFVector<T> *vec) 
773
  {
774
    vec->print();
775
  }
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807

  // 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>
808
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
809
810
811
812
813
814
815
816
817
818
819
820
821

  // 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
822
823
  inline void scal(T a, DOFVector<T>& y) 
  {
824
825
    y *= a;
  }
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841

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

842
843
844
845
846
847
848
849
850
  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);


851
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
852
853
  inline void set(DOFVector<T>& vec, T d) 
  {
854
    vec.set(d);
855
  }
856
857

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
858
859
  inline void setValue(DOFVector<T>& vec, T d) 
  {
860
    vec.set(d);
861
  }
862
863

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
864
865
  inline int size(DOFVector<T> *vec) 
  {
866
    return vec->getUsedSize();
867
  }
868

869
870
871
872
873
874
  template<typename T>
  inline int size(const DOFVector<T>& vec) 
  {
    return vec.getUsedSize();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
875
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
876
  inline void checkFeSpace(const FiniteElemSpace* feSpace, const vector<T>& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
877
878
879
880
881
882
883
884
  {
    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());
  }

885
886
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
887
888
  
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
889
890
  vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     vector<DOFVector<double>*> *res);
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915

  
  /// 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);
  }
916
  
917
  /** \brief
918
919
   * 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
920
921
   * This function works also for the case that the DOFVectors are defined on 
   * two different meshes.
922
   */
923