DOFVector.h 21.2 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
// ==                                                                        ==
// ============================================================================

/** \file DOFVector.h */

#ifndef AMDIS_DOFVECTOR_H 
#define AMDIS_DOFVECTOR_H 
 
25
26
27
28
#include <vector> 
#include <memory> 
#include <list> 

29
#include "AMDiS_fwd.h"
30
31
32
33
34
35
36
37
38
39
40
#include "FixVec.h"
#include "Global.h" 
#include "Flag.h" 
#include "RCNeighbourList.h" 
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "DOFContainer.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
Thomas Witkowski's avatar
Thomas Witkowski committed
41
#include "BasisFunction.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
42
#include "FiniteElemSpace.h"
43

44
45
46
47
48
49
namespace AMDiS {

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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
105
106
107
108
109
110
111
112
113
114
115
116
117
    /** \brief
     * Assembles the element vector for the given ellement and adds the
     * element matrix to the current DOF vector.
     */ 
    void assemble(T factor, ElInfo *elInfo,			    
		  const BoundaryType *bound, 
		  Operator *op = NULL);

    void assemble2(T factor, 
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
		   const BoundaryType *bound, 
		   Operator *op = NULL);
118
119

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

    /* \brief
     * That function must be called after the matrix assembling has been finished. 
     * This makes it possible to start some cleanup or matrix data compressing 
     * procedures.
     */
    void finishAssembling();
131
132
133
134
135
136
137
138
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
139
    }
140

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

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

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

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


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

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

    Flag getAssembleFlag();

    /** \brief
     * Evaluates \f[ u_h(x(\lambda)) = \sum_{i=0}^{m-1} vec[ind[i]] * 
     * \varphi^i(\lambda) \f] where \f$ \varphi^i \f$ is the i-th basis function,
     * \f$ x(\lambda) \f$ are the world coordinates of lambda and
     * \f$ m \f$ is the number of basis functions
     */
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
276
  /// Specifies which operation should be done after coarsening
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
  typedef enum{
    NO_OPERATION = 0,   
    COARSE_RESTRICT = 1,
    COARSE_INTERPOL = 2 
  } CoarsenOperation;
 

  /** \ingroup DOFAdministration 
   * \brief
   * The DOFs described above are just integers that can be used as indices into 
   * vectors and matrices. During refinement and coarsening of the mesh, the 
   * number of used DOFs, the meaning of one integer index, and even the total 
   * range of DOFs change. To be able to handle these changes automatically for 
   * all vectors, which are indexed by the DOFs, special data structures are 
   * used which contain such vector data. Lists of these structures are kept in 
   * DOFAdmin, so that all vectors in the lists can be resized together with the
   * range of DOFs. During refinement and coarsening of elements, values can be
   * interpolated automatically to new DOFs, and restricted from old DOFs.
   */
  template<typename T> 
  class DOFVector : public DOFVectorBase<T>, public Serializable
  {  
  public:
    /** \ingroup DOFAdministration
     * \brief
     * Enables the access of DOFVector<T>::Iterator. Alias for DOFIterator<T>
     */
    class Iterator : public DOFIterator<T> {
    public:
      Iterator(DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(c, type)
308
      {}
309
310
311

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

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

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

542
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
543
544
    void serialize(std::ostream &out) 
    {
545
546
547
      unsigned int size = vec.size();
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
548
    }
549

550
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
    void deserialize(std::istream &in) 
    {
553
554
555
556
      unsigned int size;
      in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int));
      vec.resize(size);
      in.read(reinterpret_cast<char*>(&(vec[0])), size * sizeof(T));
557
    }
558
559
560
561
562

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

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

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

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

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

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

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

581
  protected: 
582
    /// Data container
583
    std::vector<T> vec; 
584
 
585
    /// Specifies what operation should be performed after coarsening
586
587
    CoarsenOperation coarsenOperation;

588
    /// Used by \ref interpol
589
590
    AbstractFunction<T, WorldVector<double> > *interFct;

591
    /// Used for mesh traversal
592
593
594
    static DOFVector<T> *traverseVector;

  protected:
595
    /// Used while calculating vector norms
596
597
    static FastQuadrature *quad_fast;

598
    /// Stores the last calculated vector norm
599
600
601
    static double norm;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
602

603
604
605
606
607
608
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
609
610
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
611
    return v.min();
612
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
613

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


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
655
656
657
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

658
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
659
660
    int getSize() const 
    {
661
      return DOFVector<DegreeOfFreedom>::getSize();
662
    }
663

664
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
665
666
    void resize(int size) 
    {
667
668
669
670
671
672
673
674
675
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
676

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
684
685
  double L2Norm(DOFVector<T> *vec) 
  {
686
    return vec->L2Norm();
687
  }
688
689

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
696
697
  void print(DOFVector<T> *vec) 
{
698
    vec->print();
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
733
734
735
736
737
738
739
740
741
742
743
744
745

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

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

766
767
768
769
770
771
772
773
774
  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);


775
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
776
777
  inline void set(DOFVector<T>& vec, T d) 
  {
778
    vec.set(d);
779
  }
780
781

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
788
789
  inline int size(DOFVector<T> *vec) 
  {
790
    return vec->getUsedSize();
791
  }
792

Thomas Witkowski's avatar
Thomas Witkowski committed
793
794
795
796
797
798
799
800
801
802
  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());
  }

803
804
805
806
807
808
809
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_