DOFVector.h 21.9 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
// ==                                                                        ==
// ============================================================================

/** \file DOFVector.h */

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

29
#include "AMDiS_fwd.h"
30
31
32
33
34
35
36
37
38
39
40
#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
41
#include "BasisFunction.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
42
#include "FiniteElemSpace.h"
43

44
45
46
47
48
49
namespace AMDiS {

  template<typename T> 
  class DOFVectorBase : public DOFIndexed<T>
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
50

51
52
    DOFVectorBase() 
      : feSpace(NULL),
53
	elementVector(3),
Thomas Witkowski's avatar
Thomas Witkowski committed
54
55
        boundaryManager(NULL),
        nBasFcts(0)
56
57
58
59
60
    {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      rankDofs = NULL;
#endif
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
61
    
62
    DOFVectorBase(const FiniteElemSpace *f, std::string n);
63

Thomas Witkowski's avatar
Thomas Witkowski committed
64
    virtual ~DOFVectorBase();
65

66
67
68
    /// Creates temporary used data structures.
    void createTempData();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
76
77
78
79
    /** \brief
     * 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

91
92
    const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo, 
				      const Quadrature *quad,
93
				      const FastQuadrature *quadFast,
94
				      WorldVector<T> *grdAtQPs) const;
95

Thomas Witkowski's avatar
Thomas Witkowski committed
96
97
98
99
100
101
    const WorldVector<T> *getGrdAtQPs(const ElInfo *smallElInfo, 
				      const ElInfo *largeElInfo,
				      const Quadrature *quad,
				      const FastQuadrature *quadFast,
				      WorldVector<T> *grdAtQPs) const;

102
103
    const WorldMatrix<T> *getD2AtQPs(const ElInfo *elInfo, 
				     const Quadrature *quad,
104
				     const FastQuadrature *quadFast,
105
				     WorldMatrix<T> *d2AtQPs) const;
106

Thomas Witkowski's avatar
Thomas Witkowski committed
107
    /// Returns the FE space the DOF vector is defined on.
108
    inline const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
109
    {
110
      return feSpace;
111
    }
112

Thomas Witkowski's avatar
Thomas Witkowski committed
113
114
115
116
117
118
119
120
121
122
123
124
125
    /** \brief
     * Assembles the element vector for the given ellement and adds the
     * element matrix to the current DOF vector.
     */ 
    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);
126
127

    void addElementVector(T sign,
128
			  const ElementVector& elVec, 
129
			  const BoundaryType *bound,
130
			  ElInfo *elInfo,
131
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
132
133
134
135
136
137
138

    /* \brief
     * 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.
     */
    void finishAssembling();
139
140
141
142
143
144
145
146
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
147
    }
148

Thomas Witkowski's avatar
Thomas Witkowski committed
149
150
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
151
      return operatorFactor.begin();
152
    }
153

Thomas Witkowski's avatar
Thomas Witkowski committed
154
155
    inline std::vector<double*>::iterator getOperatorFactorEnd() 
    {
156
      return operatorFactor.end();
157
    }
158

Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
161
      return operatorEstFactor.begin();
162
    }
163

Thomas Witkowski's avatar
Thomas Witkowski committed
164
165
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() 
    {
166
      return operatorEstFactor.end();
167
    }
168
169


Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
172
      return operators.begin();
173
    }
174

Thomas Witkowski's avatar
Thomas Witkowski committed
175
176
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
177
      return operators.end();
178
    }
179
180
181
182
183
184
185
186
187
188
189

    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
190
191
    inline std::vector<Operator*>& getOperators() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
192
      return operators; 
193
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
194

Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
    inline std::vector<double*>& getOperatorFactor() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
197
      return operatorFactor; 
198
    }
199

Thomas Witkowski's avatar
Thomas Witkowski committed
200
201
    inline std::vector<double*>& getOperatorEstFactor() 
    { 
202
      return operatorEstFactor; 
203
    }
204

Thomas Witkowski's avatar
Thomas Witkowski committed
205
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
206
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
207
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
208
      return name; 
209
    } 
210

211
212
    inline void setName(std::string n)
    {
213
      name = n;
214
215
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
216
217
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
218
      return boundaryManager; 
219
    }
220

Thomas Witkowski's avatar
Thomas Witkowski committed
221
222
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
223
      boundaryManager = bm;
224
    }
225

Thomas Witkowski's avatar
Thomas Witkowski committed
226
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
227
    inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
228
    {
229
      rankDofs = &dofmap;
Thomas Witkowski's avatar
Thomas Witkowski committed
230
231
232
233
    }

    inline bool isRankDof(DegreeOfFreedom dof)
    {
234
      TEST_EXIT_DBG(rankDofs)("No rank dofs set!\n");
235

236
      return (*rankDofs)[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
237
238
239
    }
#endif

240
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
241
    ///
242
243
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    ///
245
    std::string name;
246

Thomas Witkowski's avatar
Thomas Witkowski committed
247
    ///
248
    ElementVector elementVector;
249

Thomas Witkowski's avatar
Thomas Witkowski committed
250
    ///
251
    std::vector<Operator*> operators;
252

Thomas Witkowski's avatar
Thomas Witkowski committed
253
    ///
254
    std::vector<double*> operatorFactor;
255

Thomas Witkowski's avatar
Thomas Witkowski committed
256
    ///
257
    std::vector<double*> operatorEstFactor;
258

Thomas Witkowski's avatar
Thomas Witkowski committed
259
    ///
260
261
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
263
    int nBasFcts;
264

Thomas Witkowski's avatar
Thomas Witkowski committed
265
    /// Are used to store temporary local dofs of an element. Thread safe.
266
267
    std::vector<DegreeOfFreedom*> localIndices;

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    /// Are used to store temporary local values of an element. Thread safe.
269
    std::vector<mtl::dense_vector<T> > localVectors;
Thomas Witkowski's avatar
Thomas Witkowski committed
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    /// Temporarly used in \ref getGrdAtQPs. Thread safe.
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
    std::vector<DimVec<double>*> grdTmp;

Thomas Witkowski's avatar
Thomas Witkowski committed
274
    /// Temporarly used in \ref getGrdAtQPs. Thread safe.
275
276
    std::vector<DimVec<double>*> grdPhis;

Thomas Witkowski's avatar
Thomas Witkowski committed
277
    /// Temporarly used in \ref getD2AtQPs. Thread safe.
278
    std::vector<DimMat<double>*> D2Phis;
Thomas Witkowski's avatar
Thomas Witkowski committed
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
284
    std::map<DegreeOfFreedom, bool> *rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
285
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
286
  };
287

Thomas Witkowski's avatar
Thomas Witkowski committed
288
  /// Specifies which operation should be done after coarsening
289
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
  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)
320
      {}
321
322
323

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

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

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

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

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

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

      this->createTempData();
369
    }
370

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
399
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
    inline CoarsenOperation getCoarsenOperation() 
    { 
402
      return coarsenOperation; 
403
    }
404

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

Thomas Witkowski's avatar
Thomas Witkowski committed
408
    /// Interpolation after refinement.
409
    inline void refineInterpol(RCNeighbourList&, int) {}
410

Thomas Witkowski's avatar
Thomas Witkowski committed
411
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
    std::vector<T>& getVector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
414
      return vec;
415
    }
416

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

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

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

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

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

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

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

462
      return vec[i];
463
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
464
 
Thomas Witkowski's avatar
Thomas Witkowski committed
465
    /// Calculates Integral of this DOFVector
466
467
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
468
    /// Calculates L1 norm of this DOFVector
469
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
470
 
Thomas Witkowski's avatar
Thomas Witkowski committed
471
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
472
473
    inline double L2Norm(Quadrature* q = NULL) const 
    {
474
      return sqrt(L2NormSquare());
475
    }
476

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

Thomas Witkowski's avatar
Thomas Witkowski committed
480
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
481
482
    inline double H1Norm(Quadrature* q = NULL) const 
    {
483
      return sqrt(H1NormSquare());
484
    }
485

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

Thomas Witkowski's avatar
Thomas Witkowski committed
489
    /// Calculates euclidian norm of this DOFVector
490
491
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
492
    /// Returns square of the euclidian norm.
493
494
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
495
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
496
497
    inline double l2norm() const 
    { 
498
      return nrm2();
499
    }
500

Thomas Witkowski's avatar
Thomas Witkowski committed
501
    /// Calculates the absolute sum of this DOFVector
502
503
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
504
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
505
506
    inline double l1norm() const 
    { 
507
      return asum();
508
    } 
509

Thomas Witkowski's avatar
Thomas Witkowski committed
510
    /// Calculates doublewell of this DOFVector
511
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
512
 
Thomas Witkowski's avatar
Thomas Witkowski committed
513
    /// Calculates the sum of this DOFVector
514
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
515
 
Thomas Witkowski's avatar
Thomas Witkowski committed
516
    /// Sets \ref vec[i] = val, i=0 , ... , size
517
518
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
519
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
520
521
    inline DOFVector<T>& operator=(T d) 
    {
522
523
      set(d); 
      return *this;
524
    } 
525

Thomas Witkowski's avatar
Thomas Witkowski committed
526
    /// Assignment operator between two vectors
527
    DOFVector<T>& operator=(const DOFVector<T>&); 
528

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

532
    /// Returns minimum of DOFVector.
533
534
    T min() const; 

535
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
536
537
    T max() const;

538
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
539
    T absMax() const;
540

541
542
543
544
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Used by interpol while mesh traversal.
545
    void interpolFct(ElInfo* elinfo);
546

547
    /// Prints \ref vec to stdout.
548
549
    void print() const; 

550
    ///
551
552
    int calcMemoryUsage() const;

553
554
555
556
557
558
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

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

561
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
562
563
    void serialize(std::ostream &out) 
    {
564
565
566
      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));
567
    }
568

569
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
570
571
    void deserialize(std::istream &in) 
    {
572
573
574
575
      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));
576
    }
577
578
579
580
581

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

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

582
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
583
584

  protected: 
585
    /// Data container
586
    std::vector<T> vec; 
587
 
588
    /// Specifies what operation should be performed after coarsening
589
590
    CoarsenOperation coarsenOperation;

591
    /// Used by \ref interpol
592
593
    AbstractFunction<T, WorldVector<double> > *interFct;

594
    /// Used for mesh traversal
595
596
597
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
598

599
600
601
602
603
604
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
605
606
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
607
    return v.min();
608
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
609

Thomas Witkowski's avatar
Thomas Witkowski committed
610
611
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
612
    return v.max();
613
  }
614
615
616
617
618
619
620
621
622
623
624
625
626
627


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
    /** \brief
     * Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
     * as DOFContainer at DOFAdmin
     */
628
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
629
630
631
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
632
    }
633
  
634
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
    ~DOFVectorDOF() 
    {
637
      feSpace->getAdmin()->removeDOFContainer(this);
638
    }
639
640
641
642
643

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVEctor<DegreeOfFreedom>::operator[]()
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
646
      return DOFVector<DegreeOfFreedom>::operator[](i);
647
    }
648

Thomas Witkowski's avatar
Thomas Witkowski committed
649
650
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
651
652
653
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

654
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
655
656
    int getSize() const 
    {
657
      return DOFVector<DegreeOfFreedom>::getSize();
658
    }
659

660
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
661
662
    void resize(int size) 
    {
663
664
665
666
667
668
669
670
671
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
672

673
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
674
675
  double norm(DOFVector<T> *vec) 
  {
676
    return vec->nrm2();
677
  }
678
679

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
680
681
  double L2Norm(DOFVector<T> *vec) 
  {
682
    return vec->L2Norm();
683
  }
684
685

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
686
687
  double H1Norm(DOFVector<T> *vec) 
  {
688
    return vec->H1Norm();
689
  }
690
691

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
  void print(DOFVector<T> *vec) 
{
694
    vec->print();
695
  }
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741

  // 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>
  void axpy(double a,const DOFVector<T>& x, DOFVector<T>& y);

  // 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
742
743
  inline void scal(T a, DOFVector<T>& y) 
  {
744
745
    y *= a;
  }
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761

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

762
763
764
765
766
767
768
769
770
  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);


771
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
772
773
  inline void set(DOFVector<T>& vec, T d) 
  {
774
    vec.set(d);
775
  }
776
777

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
778
779
  inline void setValue(DOFVector<T>& vec, T d) 
  {
780
    vec.set(d);
781
  }
782
783

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
784
785
  inline int size(DOFVector<T> *vec) 
  {
786
    return vec->getUsedSize();
787
  }
788

Thomas Witkowski's avatar
Thomas Witkowski committed
789
790
791
792
793
794
795
796
797
798
  template<typename T>
  inline void checkFeSpace(const FiniteElemSpace* feSpace, const std::vector<T>& vec)
  {
    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());
  }

799
800
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820

  /** \brief
   * Computes the integral of a function that includes two different DOFVectors. This
   * function works also for the case that the DOFVectors are defined on two different
   * meshes.
   */
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct);

  /// Computes the integral of a function with two DOFVectors defined on the same mesh.
  double intSingleMesh(const DOFVector<double> &vec1,
		       const DOFVector<double> &vec2,
		       BinaryAbstractFunction<double, double, double> *fct);

  /// Computes the integral of a function with two DOFVectors defined on different meshes.
  double intDualMesh(const DOFVector<double> &vec1,
		     const DOFVector<double> &vec2,
		     BinaryAbstractFunction<double, double, double> *fct);

Thomas Witkowski's avatar
Thomas Witkowski committed
821
822
823
  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct);

824
825
826
827
828
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_