DOFVector.h 20.7 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

Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
68
69
    /** \brief 
     * For the given element, this function returns an array of all DOFs of this
     * DOFVector that are defined on this element.
     */
70
71
72
    const T *getLocalVector(const Element *el, T* localVec) const;

    const T *getLocalVector(const ElInfo *elInfo, T* localVec) const;
73

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
281
  /// Specifies which operation should be done after coarsening
282
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
  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)
313
      {}
314
315
316

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

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

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
362
    /// Destructor
363
364
    virtual ~DOFVector();

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
    inline CoarsenOperation getCoarsenOperation() 
    { 
393
      return coarsenOperation; 
394
    }
395

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

Thomas Witkowski's avatar
Thomas Witkowski committed
399
    /// Interpolation after refinement.
400
    inline void refineInterpol(RCNeighbourList&, int) {}
401

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
445
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
    inline T& operator[](DegreeOfFreedom i) 
    {
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
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
453
 
Thomas Witkowski's avatar
Thomas Witkowski committed
454
    /// Calculates Integral of this DOFVector
455
456
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
457
    /// Calculates L1 norm of this DOFVector
458
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
459
 
Thomas Witkowski's avatar
Thomas Witkowski committed
460
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
    inline double L2Norm(Quadrature* q = NULL) const 
    {
463
      return sqrt(L2NormSquare());
464
    }
465

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

Thomas Witkowski's avatar
Thomas Witkowski committed
469
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
    inline double H1Norm(Quadrature* q = NULL) const 
    {
472
      return sqrt(H1NormSquare());
473
    }
474

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

Thomas Witkowski's avatar
Thomas Witkowski committed
478
    /// Calculates euclidian norm of this DOFVector
479
480
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
481
    /// Returns square of the euclidian norm.
482
483
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
484
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
    inline double l2norm() const 
    { 
487
      return nrm2();
488
    }
489

Thomas Witkowski's avatar
Thomas Witkowski committed
490
    /// Calculates the absolute sum of this DOFVector
491
492
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
493
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
494
495
    inline double l1norm() const 
    { 
496
      return asum();
497
    } 
498

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

Thomas Witkowski's avatar
Thomas Witkowski committed
508
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
509
510
    inline DOFVector<T>& operator=(T d) 
    {
511
512
      set(d); 
      return *this;
513
    } 
514

Thomas Witkowski's avatar
Thomas Witkowski committed
515
    /// Assignment operator between two vectors
516
    DOFVector<T>& operator=(const DOFVector<T>&); 
517

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

Thomas Witkowski's avatar
Thomas Witkowski committed
521
    /// Returns minimum of DOFVector
522
523
    T min() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
524
    /// Returns maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
525
526
    T max() const;

527
    /// Returns absolute maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
528
    T absMax() const;
529

530
    /// Used by interpol while mesh traversal
531
    void interpolFct(ElInfo* elinfo);
532

533
    /// Prints \ref vec to stdout
534
535
    void print() const; 

536
    ///
537
538
    int calcMemoryUsage() const;

539
540
541
542
543
544
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

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

547
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
548
549
    void serialize(std::ostream &out) 
    {
550
551
552
      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));
553
    }
554

555
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
556
557
    void deserialize(std::istream &in) 
    {
558
559
560
561
      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));
562
    }
563
564
565
566
567

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

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

568
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
569
570

  protected: 
571
    /// Data container
572
    std::vector<T> vec; 
573
 
574
    /// Specifies what operation should be performed after coarsening
575
576
    CoarsenOperation coarsenOperation;

577
    /// Used by \ref interpol
578
579
    AbstractFunction<T, WorldVector<double> > *interFct;

580
    /// Used for mesh traversal
581
582
583
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
584

585
586
587
588
589
590
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
591
592
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
593
    return v.min();
594
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
595

Thomas Witkowski's avatar
Thomas Witkowski committed
596
597
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
598
    return v.max();
599
  }
600
601
602
603
604
605
606
607
608
609
610
611
612
613


  /** \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
     */
614
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
615
616
617
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
618
    }
619
  
620
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
621
622
    ~DOFVectorDOF() 
    {
623
      feSpace->getAdmin()->removeDOFContainer(this);
624
    }
625
626
627
628
629

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVEctor<DegreeOfFreedom>::operator[]()
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
630
631
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
632
      return DOFVector<DegreeOfFreedom>::operator[](i);
633
    }
634

Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
637
638
639
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

640
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
641
642
    int getSize() const 
    {
643
      return DOFVector<DegreeOfFreedom>::getSize();
644
    }
645

646
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
647
648
    void resize(int size) 
    {
649
650
651
652
653
654
655
656
657
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
658

659
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
660
661
  double norm(DOFVector<T> *vec) 
  {
662
    return vec->nrm2();
663
  }
664
665

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
678
679
  void print(DOFVector<T> *vec) 
{
680
    vec->print();
681
  }
682
683
684
685
686
687
688
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

  // 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
728
729
  inline void scal(T a, DOFVector<T>& y) 
  {
730
731
    y *= a;
  }
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747

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

748
749
750
751
752
753
754
755
756
  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);


757
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
758
759
  inline void set(DOFVector<T>& vec, T d) 
  {
760
    vec.set(d);
761
  }
762
763

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
770
771
  inline int size(DOFVector<T> *vec) 
  {
772
    return vec->getUsedSize();
773
  }
774

Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
777
778
779
780
781
782
783
784
  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());
  }

785
786
787
788
789
790
791
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_