DOFVector.h 21.6 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
    virtual const T *getLocalVector(const Element *el, T* localVec) const;
74

Thomas Witkowski's avatar
Thomas Witkowski committed
75
76
77
78
    /** \brief
     * Evaluates the DOF vector at a set of quadrature points defined on the 
     * given element.
     */
79
80
    const T *getVecAtQPs(const ElInfo *elInfo, 
			 const Quadrature *quad,
81
			 const FastQuadrature *quadFast,
82
			 T *vecAtQPs) const;
83

Thomas Witkowski's avatar
Thomas Witkowski committed
84
85
86
87
88
89
    const T *getVecAtQPs(const ElInfo *smallElInfo, 
			 const ElInfo *largeElInfo,
			 const Quadrature *quad,
			 const FastQuadrature *quadFast,
			 T *vecAtQPs) const;

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
112
113
114
115
116
117
118
119
120
121
122
123
124
    /** \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);
125
126

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

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

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

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

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

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


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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
210
211
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
212
      return boundaryManager; 
213
    }
214

Thomas Witkowski's avatar
Thomas Witkowski committed
215
216
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
217
      boundaryManager = bm;
218
    }
219

Thomas Witkowski's avatar
Thomas Witkowski committed
220
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
221
    inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
222
    {
223
      rankDofs = &dofmap;
Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
226
227
    }

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

230
      return (*rankDofs)[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
233
    }
#endif

234
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
235
    ///
236
237
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
238
    ///
239
    std::string name;
240

Thomas Witkowski's avatar
Thomas Witkowski committed
241
    ///
242
    ElementVector elementVector;
243

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    ///
245
    std::vector<Operator*> operators;
246

Thomas Witkowski's avatar
Thomas Witkowski committed
247
    ///
248
    std::vector<double*> operatorFactor;
249

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

Thomas Witkowski's avatar
Thomas Witkowski committed
253
    ///
254
255
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
256
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
257
    int nBasFcts;
258

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

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    /// Are used to store temporary local values of an element. Thread safe.
Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
    std::vector<T*> localVectors;

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
274
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
275
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
276
277

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
278
    std::map<DegreeOfFreedom, bool> *rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
279
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
280
  };
281

Thomas Witkowski's avatar
Thomas Witkowski committed
282
  /// Specifies which operation should be done after coarsening
283
284
285
286
287
288
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
  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)
314
      {}
315
316
317

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
318
      {}
319
320
321
322
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
323
324
325
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
326

Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
      DOFVector<T> *create() 
      { 
329
	return new DOFVector<T>(feSpace, ""); 
330
      }
331

Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
      void free(DOFVector<T> *vec) 
      { 
334
	delete vec; 
335
      }
336
337
338
339
340
341

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
342
    /// Empty constructor. No initialization!
343
344
345
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
346
    {}
347

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

Thomas Witkowski's avatar
Thomas Witkowski committed
351
    /// Initialization.
352
    void init(const FiniteElemSpace* f, std::string n);
353

Thomas Witkowski's avatar
Thomas Witkowski committed
354
    /// Copy Constructor
Thomas Witkowski's avatar
Thomas Witkowski committed
355
356
    DOFVector(const DOFVector& rhs) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
357
      *this = rhs;   
358
359
360
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
361
362

      this->createTempData();
363
    }
364

Thomas Witkowski's avatar
Thomas Witkowski committed
365
    /// Destructor
366
367
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
368
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
    typename std::vector<T>::iterator begin() 
    { 
371
      return vec.begin(); 
372
    }
373

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
393
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
394
395
    inline CoarsenOperation getCoarsenOperation() 
    { 
396
      return coarsenOperation; 
397
    }
398

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

Thomas Witkowski's avatar
Thomas Witkowski committed
402
    /// Interpolation after refinement.
403
    inline void refineInterpol(RCNeighbourList&, int) {}
404

Thomas Witkowski's avatar
Thomas Witkowski committed
405
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
406
407
    std::vector<T>& getVector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
408
      return vec;
409
    }
410

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

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

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

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

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

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

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

456
      return vec[i];
457
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
458
 
Thomas Witkowski's avatar
Thomas Witkowski committed
459
    /// Calculates Integral of this DOFVector
460
461
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
462
    /// Calculates L1 norm of this DOFVector
463
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
464
 
Thomas Witkowski's avatar
Thomas Witkowski committed
465
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
466
467
    inline double L2Norm(Quadrature* q = NULL) const 
    {
468
      return sqrt(L2NormSquare());
469
    }
470

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

Thomas Witkowski's avatar
Thomas Witkowski committed
474
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
475
476
    inline double H1Norm(Quadrature* q = NULL) const 
    {
477
      return sqrt(H1NormSquare());
478
    }
479

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

Thomas Witkowski's avatar
Thomas Witkowski committed
483
    /// Calculates euclidian norm of this DOFVector
484
485
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
486
    /// Returns square of the euclidian norm.
487
488
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
489
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
490
491
    inline double l2norm() const 
    { 
492
      return nrm2();
493
    }
494

Thomas Witkowski's avatar
Thomas Witkowski committed
495
    /// Calculates the absolute sum of this DOFVector
496
497
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
498
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
499
500
    inline double l1norm() const 
    { 
501
      return asum();
502
    } 
503

Thomas Witkowski's avatar
Thomas Witkowski committed
504
    /// Calculates doublewell of this DOFVector
505
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
506
 
Thomas Witkowski's avatar
Thomas Witkowski committed
507
    /// Calculates the sum of this DOFVector
508
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
509
 
Thomas Witkowski's avatar
Thomas Witkowski committed
510
    /// Sets \ref vec[i] = val, i=0 , ... , size
511
512
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
513
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
514
515
    inline DOFVector<T>& operator=(T d) 
    {
516
517
      set(d); 
      return *this;
518
    } 
519

Thomas Witkowski's avatar
Thomas Witkowski committed
520
    /// Assignment operator between two vectors
521
    DOFVector<T>& operator=(const DOFVector<T>&); 
522

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

526
    /// Returns minimum of DOFVector.
527
528
    T min() const; 

529
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
530
531
    T max() const;

532
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
533
    T absMax() const;
534

535
536
537
538
    /// Returns the average value of the DOFVector.
    T average() const;

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

541
    /// Prints \ref vec to stdout.
542
543
    void print() const; 

544
    ///
545
546
    int calcMemoryUsage() const;

547
548
549
550
551
552
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

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

555
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
556
557
    void serialize(std::ostream &out) 
    {
558
559
560
      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));
561
    }
562

563
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
564
565
    void deserialize(std::istream &in) 
    {
566
567
568
569
      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));
570
    }
571
572
573
574
575

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

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

576
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
577
578

  protected: 
579
    /// Data container
580
    std::vector<T> vec; 
581
 
582
    /// Specifies what operation should be performed after coarsening
583
584
    CoarsenOperation coarsenOperation;

585
    /// Used by \ref interpol
586
587
    AbstractFunction<T, WorldVector<double> > *interFct;

588
    /// Used for mesh traversal
589
590
591
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
592

593
594
595
596
597
598
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
599
600
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
601
    return v.min();
602
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
603

Thomas Witkowski's avatar
Thomas Witkowski committed
604
605
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
606
    return v.max();
607
  }
608
609
610
611
612
613
614
615
616
617
618
619
620
621


  /** \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
     */
622
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
623
624
625
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
626
    }
627
  
628
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
    ~DOFVectorDOF() 
    {
631
      feSpace->getAdmin()->removeDOFContainer(this);
632
    }
633
634
635
636
637

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVEctor<DegreeOfFreedom>::operator[]()
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
638
639
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
640
      return DOFVector<DegreeOfFreedom>::operator[](i);
641
    }
642

Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
645
646
647
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

648
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
649
650
    int getSize() const 
    {
651
      return DOFVector<DegreeOfFreedom>::getSize();
652
    }
653

654
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
655
656
    void resize(int size) 
    {
657
658
659
660
661
662
663
664
665
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
666

667
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
668
669
  double norm(DOFVector<T> *vec) 
  {
670
    return vec->nrm2();
671
  }
672
673

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
686
687
  void print(DOFVector<T> *vec) 
{
688
    vec->print();
689
  }
690
691
692
693
694
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

  // 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
736
737
  inline void scal(T a, DOFVector<T>& y) 
  {
738
739
    y *= a;
  }
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755

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

756
757
758
759
760
761
762
763
764
  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);


765
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
766
767
  inline void set(DOFVector<T>& vec, T d) 
  {
768
    vec.set(d);
769
  }
770
771

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
772
773
  inline void setValue(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 int size(DOFVector<T> *vec) 
  {
780
    return vec->getUsedSize();
781
  }
782

Thomas Witkowski's avatar
Thomas Witkowski committed
783
784
785
786
787
788
789
790
791
792
  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());
  }

793
794
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814

  /** \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);

815
816
817
818
819
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_