DOFVector.h 21.2 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"
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
214
215
216
217
218
219
220
221
222
223
224
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    void setRankDofs(std::map<DegreeOfFreedom, bool>& dofmap)
    {
      rankDofs = dofmap;
    }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
229
    ///
230
    std::string name;
231

Thomas Witkowski's avatar
Thomas Witkowski committed
232
    ///
233
    ElementVector elementVector;
234

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    ///
245
246
    BoundaryManager *boundaryManager;

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
273
  /// Specifies which operation should be done after coarsening
274
275
276
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
  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)
305
      {}
306
307
308

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
309
      {}
310
311
312
313
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
314
315
316
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
317

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    /// Destructor
356
357
    virtual ~DOFVector();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

523
    /// Used by interpol while mesh traversal
524
    static int interpolFct(ElInfo* elinfo);
525

526
    /// Prints \ref vec to stdout
527
528
    void print() const; 

529
    ///
530
531
    int calcMemoryUsage() const;

532
533
534
535
536
537
538
539
    /** \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);

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

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

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

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

561
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
562

Thomas Witkowski's avatar
Thomas Witkowski committed
563
  protected: 
Thomas Witkowski's avatar
Thomas Witkowski committed
564
    /// Used by Int while mesh traversal
565
566
    static int Int_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
567
    /// Used by L1Norm while mesh traversal
568
    static int L1Norm_fct(ElInfo* elinfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
569
 
Thomas Witkowski's avatar
Thomas Witkowski committed
570
    /// Used by L2Norm while mesh traversal
571
572
    static int L2Norm_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
573
    /// Used by H1Norm while mesh traversal
574
575
    static int H1Norm_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
576
    /// Used by DoubleWell while mesh traversal
577
578
    static int DoubleWell_fct(ElInfo* elinfo);

579
  protected: 
580
    /// Name of this DOFVector
581
    std::string name; 
582

583
    /// FiniteElemSpace of the vector
584
585
    const FiniteElemSpace *feSpace; 

586
    /// Data container
587
    std::vector<T> vec; 
588
 
589
    /// Specifies what operation should be performed after coarsening
590
591
    CoarsenOperation coarsenOperation;

592
    /// Used by \ref interpol
593
594
    AbstractFunction<T, WorldVector<double> > *interFct;

595
    /// Used for mesh traversal
596
597
598
    static DOFVector<T> *traverseVector;

  protected:
599
    /// Used while calculating vector norms
600
601
    static FastQuadrature *quad_fast;

602
    /// Stores the last calculated vector norm
603
604
605
    static double norm;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
606

607
608
609
610
611
612
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
613
614
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
615
    return v.min();
616
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
617

Thomas Witkowski's avatar
Thomas Witkowski committed
618
619
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
620
    return v.max();
621
  }
622
623
624
625
626
627
628
629
630
631
632
633
634
635


  /** \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
     */
636
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
637
638
639
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
640
    }
641
  
642
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
    ~DOFVectorDOF() 
    {
645
      feSpace->getAdmin()->removeDOFContainer(this);
646
    }
647
648
649
650
651

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVEctor<DegreeOfFreedom>::operator[]()
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
652
653
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
654
      return DOFVector<DegreeOfFreedom>::operator[](i);
655
    }
656

Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
659
660
661
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

662
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
663
664
    int getSize() const 
    {
665
      return DOFVector<DegreeOfFreedom>::getSize();
666
    }
667

668
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
669
670
    void resize(int size) 
    {
671
672
673
674
675
676
677
678
679
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
680

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
688
689
  double L2Norm(DOFVector<T> *vec) 
  {
690
    return vec->L2Norm();
691
  }
692
693

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
694
695
  double H1Norm(DOFVector<T> *vec) 
  {
696
    return vec->H1Norm();
697
  }
698
699

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
700
701
  void print(DOFVector<T> *vec) 
{
702
    vec->print();
703
  }
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749

  // 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
750
751
  inline void scal(T a, DOFVector<T>& y) 
  {
752
753
    y *= a;
  }
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769

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

770
771
772
773
774
775
776
777
778
  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);


779
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
780
781
  inline void set(DOFVector<T>& vec, T d) 
  {
782
    vec.set(d);
783
  }
784
785

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
786
787
  inline void setValue(DOFVector<T>& vec, T d) 
  {
788
    vec.set(d);
789
  }
790
791

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
792
793
  inline int size(DOFVector<T> *vec) 
  {
794
    return vec->getUsedSize();
795
  }
796

Thomas Witkowski's avatar
Thomas Witkowski committed
797
798
799
800
801
802
803
804
805
806
  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());
  }

807
808
809
810
811
812
813
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_