DOFVector.h 21.1 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
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
337
    {}
338

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

585
    /// Used by \ref interpol
586
587
    AbstractFunction<T, WorldVector<double> > *interFct;

588
    /// Used for mesh traversal
589
590
591
    static DOFVector<T> *traverseVector;

  protected:
592
    /// Used while calculating vector norms
593
594
    static FastQuadrature *quad_fast;

595
    /// Stores the last calculated vector norm
596
597
598
    static double norm;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
599

600
601
602
603
604
605
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

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

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


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

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
673

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

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

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

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

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

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

763
764
765
766
767
768
769
770
771
  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);


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

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

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

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

800
801
802
803
804
805
806
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_