DOFVector.h 20.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
    {}
Thomas Witkowski's avatar
Thomas Witkowski committed
57
    
58
    DOFVectorBase(const FiniteElemSpace *f, std::string n);
59

Thomas Witkowski's avatar
Thomas Witkowski committed
60
    virtual ~DOFVectorBase();
61

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

Thomas Witkowski's avatar
Thomas Witkowski committed
68
69
70
71
    /** \brief
     * Evaluates the DOF vector at a set of quadrature points defined on the 
     * given element.
     */
72
73
    const T *getVecAtQPs(const ElInfo *elInfo, 
			 const Quadrature *quad,
74
			 const FastQuadrature *quadFast,
75
			 T *vecAtQPs) const;
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
78
79
80
81
82
    const T *getVecAtQPs(const ElInfo *smallElInfo, 
			 const ElInfo *largeElInfo,
			 const Quadrature *quad,
			 const FastQuadrature *quadFast,
			 T *vecAtQPs) const;

83
84
    const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo, 
				      const Quadrature *quad,
85
				      const FastQuadrature *quadFast,
86
				      WorldVector<T> *grdAtQPs) const;
87

Thomas Witkowski's avatar
Thomas Witkowski committed
88
89
90
91
92
93
    const WorldVector<T> *getGrdAtQPs(const ElInfo *smallElInfo, 
				      const ElInfo *largeElInfo,
				      const Quadrature *quad,
				      const FastQuadrature *quadFast,
				      WorldVector<T> *grdAtQPs) const;

94
95
    const WorldMatrix<T> *getD2AtQPs(const ElInfo *elInfo, 
				     const Quadrature *quad,
96
				     const FastQuadrature *quadFast,
97
				     WorldMatrix<T> *d2AtQPs) const;
98

Thomas Witkowski's avatar
Thomas Witkowski committed
99
    /// Returns the FE space the DOF vector is defined on.
Thomas Witkowski's avatar
Thomas Witkowski committed
100
101
    inline const FiniteElemSpace* getFESpace() const 
    {
102
      return feSpace;
103
    }
104

Thomas Witkowski's avatar
Thomas Witkowski committed
105
106
107
108
109
110
111
112
113
114
115
116
117
    /** \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);
118
119

    void addElementVector(T sign,
120
			  const ElementVector& elVec, 
121
			  const BoundaryType *bound,
122
			  ElInfo *elInfo,
123
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
124
125
126
127
128
129
130

    /* \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();
131
132
133
134
135
136
137
138
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
139
    }
140

Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
143
      return operatorFactor.begin();
144
    }
145

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

Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
153
      return operatorEstFactor.begin();
154
    }
155

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


Thomas Witkowski's avatar
Thomas Witkowski committed
162
163
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
164
      return operators.begin();
165
    }
166

Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
169
      return operators.end();
170
    }
171
172
173
174
175
176
177
178
179
180
181

    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
182
183
    inline std::vector<Operator*>& getOperators() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
184
      return operators; 
185
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
186

Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
    inline std::vector<double*>& getOperatorFactor() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
189
      return operatorFactor; 
190
    }
191

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

Thomas Witkowski's avatar
Thomas Witkowski committed
197
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
198
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
199
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
200
      return name; 
201
    } 
202

Thomas Witkowski's avatar
Thomas Witkowski committed
203
204
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
205
      return boundaryManager; 
206
    }
207

Thomas Witkowski's avatar
Thomas Witkowski committed
208
209
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
210
      boundaryManager = bm;
211
    }
212

Thomas Witkowski's avatar
Thomas Witkowski committed
213
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
214
    inline void setRankDofs(std::map<DegreeOfFreedom, bool> dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
215
    {
216
217
218
      rankDofs.clear();

      //      rankDofs = dofmap;
Thomas Witkowski's avatar
Thomas Witkowski committed
219
220
221
222
223
224
225
226
    }

    inline bool isRankDof(DegreeOfFreedom dof)
    {
      return rankDofs[dof];
    }
#endif

227
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
228
    ///
229
230
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
231
    ///
232
    std::string name;
233

Thomas Witkowski's avatar
Thomas Witkowski committed
234
    ///
235
    ElementVector elementVector;
236

Thomas Witkowski's avatar
Thomas Witkowski committed
237
    ///
238
    std::vector<Operator*> operators;
239

Thomas Witkowski's avatar
Thomas Witkowski committed
240
    ///
241
    std::vector<double*> operatorFactor;
242

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

Thomas Witkowski's avatar
Thomas Witkowski committed
246
    ///
247
248
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
249
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
250
    int nBasFcts;
251

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
267
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
268
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
271
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
    std::map<DegreeOfFreedom, bool> rankDofs;
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
274
  };
275

Thomas Witkowski's avatar
Thomas Witkowski committed
276
  /// Specifies which operation should be done after coarsening
277
278
279
280
281
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
  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)
308
      {}
309
310
311

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
312
      {}
313
314
315
316
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
317
318
319
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
320

Thomas Witkowski's avatar
Thomas Witkowski committed
321
322
      DOFVector<T> *create() 
      { 
323
	return new DOFVector<T>(feSpace, ""); 
324
      }
325

Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
      void free(DOFVector<T> *vec) 
      { 
328
	delete vec; 
329
      }
330
331
332
333
334
335

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    /// Empty constructor. No initialization!
337
338
339
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
340
    {}
341

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

Thomas Witkowski's avatar
Thomas Witkowski committed
345
    /// Initialization.
346
    void init(const FiniteElemSpace* f, std::string n);
347

Thomas Witkowski's avatar
Thomas Witkowski committed
348
    /// Copy Constructor
Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
    DOFVector(const DOFVector& rhs) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
351
      *this = rhs;   
352
353
354
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
355
    }
356

Thomas Witkowski's avatar
Thomas Witkowski committed
357
    /// Destructor
358
359
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
360
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
    typename std::vector<T>::iterator begin() 
    { 
363
      return vec.begin(); 
364
    }
365

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

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
380
381
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
382
      coarsenOperation = op; 
383
    }
384

Thomas Witkowski's avatar
Thomas Witkowski committed
385
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
    inline CoarsenOperation getCoarsenOperation() 
    { 
388
      return coarsenOperation; 
389
    }
390

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

Thomas Witkowski's avatar
Thomas Witkowski committed
394
    /// Interpolation after refinement.
395
    inline void refineInterpol(RCNeighbourList&, int) {}
396

Thomas Witkowski's avatar
Thomas Witkowski committed
397
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
    std::vector<T>& getVector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
400
      return vec;
401
    }
402

Thomas Witkowski's avatar
Thomas Witkowski committed
403
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
404
405
    inline int getSize() const 
    { 
406
      return vec.size();
407
    } 
408

Thomas Witkowski's avatar
Thomas Witkowski committed
409
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
410
411
    inline int getUsedSize() const 
    { 
412
      return this->feSpace->getAdmin()->getUsedSize(); 
413
    }
414

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
440
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
441
442
    inline T& operator[](DegreeOfFreedom i) 
    {
443
      FUNCNAME("DOFVector<T>::operator[]");
444
445
      TEST_EXIT_DBG((i >= 0) && (i < static_cast<int>(vec.size())))
	("Illegal vector index %d.\n",i);
446
      return vec[i];
447
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
448
 
Thomas Witkowski's avatar
Thomas Witkowski committed
449
    /// Calculates Integral of this DOFVector
450
451
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
452
    /// Calculates L1 norm of this DOFVector
453
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
454
 
Thomas Witkowski's avatar
Thomas Witkowski committed
455
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
    inline double L2Norm(Quadrature* q = NULL) const 
    {
458
      return sqrt(L2NormSquare());
459
    }
460

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

Thomas Witkowski's avatar
Thomas Witkowski committed
464
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
465
466
    inline double H1Norm(Quadrature* q = NULL) const 
    {
467
      return sqrt(H1NormSquare());
468
    }
469

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

Thomas Witkowski's avatar
Thomas Witkowski committed
473
    /// Calculates euclidian norm of this DOFVector
474
475
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
476
    /// Returns square of the euclidian norm.
477
478
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
479
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
480
481
    inline double l2norm() const 
    { 
482
      return nrm2();
483
    }
484

Thomas Witkowski's avatar
Thomas Witkowski committed
485
    /// Calculates the absolute sum of this DOFVector
486
487
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
488
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
489
490
    inline double l1norm() const 
    { 
491
      return asum();
492
    } 
493

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

Thomas Witkowski's avatar
Thomas Witkowski committed
503
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
    inline DOFVector<T>& operator=(T d) 
    {
506
507
      set(d); 
      return *this;
508
    } 
509

Thomas Witkowski's avatar
Thomas Witkowski committed
510
    /// Assignment operator between two vectors
511
    DOFVector<T>& operator=(const DOFVector<T>&); 
512

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

Thomas Witkowski's avatar
Thomas Witkowski committed
516
    /// Returns minimum of DOFVector
517
518
    T min() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
519
    /// Returns maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
520
521
    T max() const;

522
    /// Returns absolute maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
523
    T absMax() const;
524

525
    /// Used by interpol while mesh traversal
526
    void interpolFct(ElInfo* elinfo);
527

528
    /// Prints \ref vec to stdout
529
530
    void print() const; 

531
    ///
532
533
    int calcMemoryUsage() const;

534
535
536
537
538
539
540
541
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

    void interpol(DOFVector<T> *v, double factor);

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

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

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

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

563
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
564
565

  protected: 
566
    /// Data container
567
    std::vector<T> vec; 
568
 
569
    /// Specifies what operation should be performed after coarsening
570
571
    CoarsenOperation coarsenOperation;

572
    /// Used by \ref interpol
573
574
    AbstractFunction<T, WorldVector<double> > *interFct;

575
    /// Used for mesh traversal
576
577
578
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
579

580
581
582
583
584
585
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
586
587
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
588
    return v.min();
589
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
590

Thomas Witkowski's avatar
Thomas Witkowski committed
591
592
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
593
    return v.max();
594
  }
595
596
597
598
599
600
601
602
603
604
605
606
607
608


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

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVEctor<DegreeOfFreedom>::operator[]()
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
625
626
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
627
      return DOFVector<DegreeOfFreedom>::operator[](i);
628
    }
629

Thomas Witkowski's avatar
Thomas Witkowski committed
630
631
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
632
633
634
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

635
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
636
637
    int getSize() const 
    {
638
      return DOFVector<DegreeOfFreedom>::getSize();
639
    }
640

641
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
642
643
    void resize(int size) 
    {
644
645
646
647
648
649
650
651
652
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
653

654
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
655
656
  double norm(DOFVector<T> *vec) 
  {
657
    return vec->nrm2();
658
  }
659
660

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
673
674
  void print(DOFVector<T> *vec) 
{
675
    vec->print();
676
  }
677
678
679
680
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

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

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

743
744
745
746
747
748
749
750
751
  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);


752
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
753
754
  inline void set(DOFVector<T>& vec, T d) 
  {
755
    vec.set(d);
756
  }
757
758

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
765
766
  inline int size(DOFVector<T> *vec) 
  {
767
    return vec->getUsedSize();
768
  }
769

Thomas Witkowski's avatar
Thomas Witkowski committed
770
771
772
773
774
775
776
777
778
779
  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());
  }

780
781
782
783
784
785
786
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_