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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
235
    ///
236
    std::string name;
237

Thomas Witkowski's avatar
Thomas Witkowski committed
238
    ///
239
    ElementVector elementVector;
240

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
250
    ///
251
252
    BoundaryManager *boundaryManager;

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
279
  /// Specifies which operation should be done after coarsening
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
308
309
310
  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)
311
      {}
312
313
314

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

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

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
360
    /// Destructor
361
362
    virtual ~DOFVector();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

531
    /// Prints \ref vec to stdout
532
533
    void print() const; 

534
    ///
535
536
    int calcMemoryUsage() const;

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

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

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

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

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

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

566
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
567
568

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
582

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

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

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

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


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

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
656

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

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
  void print(DOFVector<T> *vec) 
{
678
    vec->print();
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
723
724
725

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

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

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


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

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

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

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

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

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_