DOFVector.h 20.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \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"
42

43
44
45
46
47
48
namespace AMDiS {

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

50
51
    DOFVectorBase() 
      : feSpace(NULL),
52
	elementVector(3),
Thomas Witkowski's avatar
Thomas Witkowski committed
53
54
        boundaryManager(NULL),
        nBasFcts(0)
55
    {}
Thomas Witkowski's avatar
Thomas Witkowski committed
56
    
57
    DOFVectorBase(const FiniteElemSpace *f, std::string n);
58

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

Thomas Witkowski's avatar
Thomas Witkowski committed
61
62
63
64
65
    /** \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;
66

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
213
    ///
214
215
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
216
    ///
217
    std::string name;
218

Thomas Witkowski's avatar
Thomas Witkowski committed
219
    ///
220
    ElementVector elementVector;
221

Thomas Witkowski's avatar
Thomas Witkowski committed
222
    ///
223
    std::vector<Operator*> operators;
224

Thomas Witkowski's avatar
Thomas Witkowski committed
225
    ///
226
    std::vector<double*> operatorFactor;
227

Thomas Witkowski's avatar
Thomas Witkowski committed
228
    ///
229
    std::vector<double*> operatorEstFactor;
230

Thomas Witkowski's avatar
Thomas Witkowski committed
231
    ///
232
233
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
234
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
235
    int nBasFcts;
236

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
252
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
253
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
254
  };
255

Thomas Witkowski's avatar
Thomas Witkowski committed
256
  /// Specifies which operation should be done after coarsening
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
  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)
288
      {}
289
290
291

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
292
      {}
293
294
295
296
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
297
298
299
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
300

Thomas Witkowski's avatar
Thomas Witkowski committed
301
302
      DOFVector<T> *create() 
      { 
303
	return new DOFVector<T>(feSpace, ""); 
304
      }
305

Thomas Witkowski's avatar
Thomas Witkowski committed
306
307
      void free(DOFVector<T> *vec) 
      { 
308
	delete vec; 
309
      }
310
311
312
313
314
315

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
316
    /// Empty constructor. No initialization!
317
318
319
320
    DOFVector() 
      : DOFVectorBase<T>(), 
	feSpace(NULL),
	coarsenOperation(NO_OPERATION)
321
    {}
322

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

Thomas Witkowski's avatar
Thomas Witkowski committed
326
    /// Initialization.
327
    void init(const FiniteElemSpace* f, std::string n);
328

Thomas Witkowski's avatar
Thomas Witkowski committed
329
    /// Copy Constructor
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
    DOFVector(const DOFVector& rhs) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
      *this = rhs;   
      name = rhs.name + "copy";
Thomas Witkowski's avatar
Thomas Witkowski committed
334
      if (feSpace && feSpace->getAdmin())
335
	(dynamic_cast<DOFAdmin*>(feSpace->getAdmin()))->addDOFIndexed(this);
336
    }
337

Thomas Witkowski's avatar
Thomas Witkowski committed
338
    /// Destructor
339
340
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
341
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
    typename std::vector<T>::iterator begin() 
    { 
344
      return vec.begin(); 
345
    }
346

Thomas Witkowski's avatar
Thomas Witkowski committed
347
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
348
349
    typename std::vector<T>::iterator end() 
    { 
350
      return vec.end(); 
351
    }
352
353
354
355
356
357
  
    /** \brief
     * Used by DOFAdmin to compress this DOFVector. Implementation of
     * DOFIndexedBase::compress()
     */
    virtual void compressDOFIndexed(int first, int last,
358
				    std::vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
363
      coarsenOperation = op; 
364
    }
365

Thomas Witkowski's avatar
Thomas Witkowski committed
366
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
367
368
    inline CoarsenOperation getCoarsenOperation() 
    { 
369
      return coarsenOperation; 
370
    }
371

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

Thomas Witkowski's avatar
Thomas Witkowski committed
375
    /// Interpolation after refinement.
376
    inline void refineInterpol(RCNeighbourList&, int) {}
377

Thomas Witkowski's avatar
Thomas Witkowski committed
378
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
    std::vector<T>& getVector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
381
      return vec;
382
    }
383

Thomas Witkowski's avatar
Thomas Witkowski committed
384
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
    inline int getSize() const 
    { 
387
      return vec.size();
388
    } 
389

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
    inline int getUsedSize() const 
    { 
393
      return feSpace->getAdmin()->getUsedSize(); 
394
    }
395

Thomas Witkowski's avatar
Thomas Witkowski committed
396
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
397
398
    inline void resize(int n) 
    { 
399
400
      FUNCNAME("DOFVector<T>::resize()");
      TEST_EXIT_DBG((n >= 0))("Can't resize DOFVector to negative size\n"); 
401
      vec.resize(n);
402
    } 
403

Thomas Witkowski's avatar
Thomas Witkowski committed
404
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
    inline void resize(int n, T init) 
    { 
407
408
      FUNCNAME("DOFVector<T>::resize()");
      TEST_EXIT_DBG((n >= 0))("Can't resize DOFVector to negative size\n"); 
409
      vec.resize(n, init);
410
    } 
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
    inline const T& operator[](DegreeOfFreedom i) const 
    {
415
      FUNCNAME("DOFVector<T>::operator[]");
416
417
      TEST_EXIT_DBG((i>= 0) && (i < static_cast<int>(vec.size())))
	("Illegal vector index %d.\n",i);
418
      return vec[i];
419
    } 
420

Thomas Witkowski's avatar
Thomas Witkowski committed
421
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
    inline T& operator[](DegreeOfFreedom i) 
    {
424
      FUNCNAME("DOFVector<T>::operator[]");
425
426
      TEST_EXIT_DBG((i >= 0) && (i < static_cast<int>(vec.size())))
	("Illegal vector index %d.\n",i);
427
      return vec[i];
428
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
429
 
Thomas Witkowski's avatar
Thomas Witkowski committed
430
    /// Calculates Integral of this DOFVector
431
432
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
433
    /// Calculates L1 norm of this DOFVector
434
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
435
 
Thomas Witkowski's avatar
Thomas Witkowski committed
436
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
437
438
    inline double L2Norm(Quadrature* q = NULL) const 
    {
439
      return sqrt(L2NormSquare());
440
    }
441

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

Thomas Witkowski's avatar
Thomas Witkowski committed
445
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
    inline double H1Norm(Quadrature* q = NULL) const 
    {
448
      return sqrt(H1NormSquare());
449
    }
450

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

Thomas Witkowski's avatar
Thomas Witkowski committed
454
    /// Calculates euclidian norm of this DOFVector
455
456
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
457
    /// Returns square of the euclidian norm.
458
459
    double squareNrm2() const;

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() const 
    { 
463
      return nrm2();
464
    }
465

Thomas Witkowski's avatar
Thomas Witkowski committed
466
    /// Calculates the absolute sum of this DOFVector
467
468
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
469
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
    inline double l1norm() const 
    { 
472
      return asum();
473
    } 
474

Thomas Witkowski's avatar
Thomas Witkowski committed
475
    /// Calculates doublewell of this DOFVector
476
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
477
 
Thomas Witkowski's avatar
Thomas Witkowski committed
478
    /// Calculates the sum of this DOFVector
479
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
480
 
Thomas Witkowski's avatar
Thomas Witkowski committed
481
    /// Sets \ref vec[i] = val, i=0 , ... , size
482
483
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
484
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
    inline DOFVector<T>& operator=(T d) 
    {
487
488
      set(d); 
      return *this;
489
    } 
490

Thomas Witkowski's avatar
Thomas Witkowski committed
491
    /// Assignment operator between two vectors
492
    DOFVector<T>& operator=(const DOFVector<T>&); 
493

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

Thomas Witkowski's avatar
Thomas Witkowski committed
497
    /// Returns minimum of DOFVector
498
499
    T min() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
500
    /// Returns maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
501
502
    T max() const;

503
    /// Returns absolute maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
504
    T absMax() const;
505

506
    /// Used by interpol while mesh traversal
507
    static int interpolFct(ElInfo* elinfo);
508

509
    /// Prints \ref vec to stdout
510
511
    void print() const; 

512
    ///
513
514
    int calcMemoryUsage() const;

515
516
517
518
519
520
521
522
    /** \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);

523
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
524
525
    void serialize(std::ostream &out) 
    {
526
527
528
      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));
529
    }
530

531
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
532
533
    void deserialize(std::istream &in) 
    {
534
535
536
537
      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));
538
    }
539
540
541
542
543

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

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

544
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
545

Thomas Witkowski's avatar
Thomas Witkowski committed
546
  protected: 
Thomas Witkowski's avatar
Thomas Witkowski committed
547
    /// Used by Int while mesh traversal
548
549
    static int Int_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
550
    /// Used by L1Norm while mesh traversal
551
    static int L1Norm_fct(ElInfo* elinfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
552
 
Thomas Witkowski's avatar
Thomas Witkowski committed
553
    /// Used by L2Norm while mesh traversal
554
555
    static int L2Norm_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
556
    /// Used by H1Norm while mesh traversal
557
558
    static int H1Norm_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
559
    /// Used by DoubleWell while mesh traversal
560
561
    static int DoubleWell_fct(ElInfo* elinfo);

562
  protected: 
563
    /// Name of this DOFVector
564
    std::string name; 
565

566
    /// FiniteElemSpace of the vector
567
568
    const FiniteElemSpace *feSpace; 

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;

  protected:
582
    /// Used while calculating vector norms
583
584
    static FastQuadrature *quad_fast;

585
    /// Stores the last calculated vector norm
586
587
588
    static double norm;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
589

590
591
592
593
594
595
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

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

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


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

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
663

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

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
683
684
  void print(DOFVector<T> *vec) 
{
685
    vec->print();
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
728
729
730
731
732

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

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

753
754
755
756
757
758
759
760
761
  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);


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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
  inline int size(DOFVector<T> *vec) 
  {
777
    return vec->getUsedSize();
778
  }
779
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_