DOFVector.h 20.7 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
41
#include "FixVec.h"
#include "Global.h" 
#include "Flag.h" 
#include "RCNeighbourList.h" 
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "DOFContainer.h"
#include "MemoryManager.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
Thomas Witkowski's avatar
Thomas Witkowski committed
42
#include "BasisFunction.h"
43
44
45
46

#ifdef HAVE_PARALLEL_AMDIS
#include "petscao.h"
#endif
47
48
49
50
51
52
53
 
namespace AMDiS {

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

55
56
57
    DOFVectorBase() 
      : feSpace(NULL),
	elementVector(NULL),
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
        boundaryManager(NULL),
        nBasFcts(0)
60
    {}
Thomas Witkowski's avatar
Thomas Witkowski committed
61
    
62
    DOFVectorBase(const FiniteElemSpace *f, std::string n);
63

Thomas Witkowski's avatar
Thomas Witkowski committed
64
    virtual ~DOFVectorBase();
65

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

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

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

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

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

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

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

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

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

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

143
    inline std::vector<double*>::iterator getOperatorFactorBegin() {
144
      return operatorFactor.begin();
145
    }
146

147
    inline std::vector<double*>::iterator getOperatorFactorEnd() {
148
      return operatorFactor.end();
149
    }
150

151
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
152
      return operatorEstFactor.begin();
153
    }
154

155
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
156
      return operatorEstFactor.end();
157
    }
158
159


160
    inline std::vector<Operator*>::iterator getOperatorsBegin() {
161
      return operators.begin();
162
    }
163

164
    inline std::vector<Operator*>::iterator getOperatorsEnd() {
165
      return operators.end();
166
    }
167
168
169
170
171
172
173
174
175
176
177

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

178
    inline std::vector<Operator*>& getOperators() { 
Thomas Witkowski's avatar
Thomas Witkowski committed
179
      return operators; 
180
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
181

182
    inline std::vector<double*>& getOperatorFactor() { 
Thomas Witkowski's avatar
Thomas Witkowski committed
183
      return operatorFactor; 
184
    }
185

186
    inline std::vector<double*>& getOperatorEstFactor() { 
187
      return operatorEstFactor; 
188
    }
189

Thomas Witkowski's avatar
Thomas Witkowski committed
190
    /// Returns \ref name
191
    inline const std::string& getName() const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
192
      return name; 
193
    } 
194

Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
    inline BoundaryManager* getBoundaryManager() const { 
      return boundaryManager; 
197
    }
198
199
200

    inline void setBoundaryManager(BoundaryManager *bm) {
      boundaryManager = bm;
201
    }
202
203

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
204
    ///
205
206
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
207
    ///
208
    std::string name;
209

Thomas Witkowski's avatar
Thomas Witkowski committed
210
    ///
211
212
    ElementVector *elementVector;

Thomas Witkowski's avatar
Thomas Witkowski committed
213
    ///
214
    std::vector<Operator*> operators;
215

Thomas Witkowski's avatar
Thomas Witkowski committed
216
    ///
217
    std::vector<double*> operatorFactor;
218

Thomas Witkowski's avatar
Thomas Witkowski committed
219
    ///
220
    std::vector<double*> operatorEstFactor;
221

Thomas Witkowski's avatar
Thomas Witkowski committed
222
    ///
223
224
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
225
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
226
    int nBasFcts;
227

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
243
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
244
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
245
  };
246

Thomas Witkowski's avatar
Thomas Witkowski committed
247
  /// Specifies which operation should be done after coarsening
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
  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:
    MEMORY_MANAGED(DOFVector<T>);

    /** \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)
281
      {}
282
283
284

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
285
      {}
286
287
288
289
290
291
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
      MEMORY_MANAGED(Creator);

292
293
294
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
295
296
297

      DOFVector<T> *create() { 
	return NEW DOFVector<T>(feSpace, ""); 
298
      }
299
300
301

      void free(DOFVector<T> *vec) { 
	DELETE vec; 
302
      }
303
304
305
306
307
308

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
309
    /// Empty constructor. No initialization!
310
311
    DOFVector() 
      : DOFVectorBase<T>(), 
312
	//elementVector(NULL), 
313
314
	feSpace(NULL),
	coarsenOperation(NO_OPERATION)
315
    {}
316

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

Thomas Witkowski's avatar
Thomas Witkowski committed
320
    /// Initialization.
321
    void init(const FiniteElemSpace* f, std::string n);
322

Thomas Witkowski's avatar
Thomas Witkowski committed
323
    /// Copy Constructor
324
    DOFVector(const DOFVector& rhs) {
Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
327
      *this = rhs;   
      name = rhs.name + "copy";
      if (feSpace && feSpace->getAdmin()) {
328
329
	(dynamic_cast<DOFAdmin*>(feSpace->getAdmin()))->addDOFIndexed(this);
      }
330
    }
331

Thomas Witkowski's avatar
Thomas Witkowski committed
332
    /// Destructor
333
334
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
335
    /// Returns iterator to the begin of \ref vec
336
    typename std::vector<T>::iterator begin() { 
337
      return vec.begin(); 
338
    }
339

Thomas Witkowski's avatar
Thomas Witkowski committed
340
    /// Returns iterator to the end of \ref vec
341
    typename std::vector<T>::iterator end() { 
342
      return vec.end(); 
343
    }
344
345
346
347
348
349
  
    /** \brief
     * Used by DOFAdmin to compress this DOFVector. Implementation of
     * DOFIndexedBase::compress()
     */
    virtual void compressDOFIndexed(int first, int last,
350
				    std::vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352

    /// Sets \ref coarsenOperation to op
353
354
    inline void setCoarsenOperation(CoarsenOperation op) { 
      coarsenOperation = op; 
355
    }
356

Thomas Witkowski's avatar
Thomas Witkowski committed
357
    /// Returns \ref coarsenOperation
358
359
    inline CoarsenOperation getCoarsenOperation() { 
      return coarsenOperation; 
360
    }
361

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

Thomas Witkowski's avatar
Thomas Witkowski committed
365
    /// Interpolation after refinement.
366
    inline void refineInterpol(RCNeighbourList&, int) {}
367

Thomas Witkowski's avatar
Thomas Witkowski committed
368
    /// Returns \ref vec
369
    std::vector<T>& getVector() { 
Thomas Witkowski's avatar
Thomas Witkowski committed
370
      return vec;
371
    }
372

Thomas Witkowski's avatar
Thomas Witkowski committed
373
    /// Returns size of \ref vec
374
375
    inline int getSize() const { 
      return vec.size();
376
    } 
377

Thomas Witkowski's avatar
Thomas Witkowski committed
378
    /// Returns used size of the vector.
379
380
    inline int getUsedSize() const { 
      return feSpace->getAdmin()->getUsedSize(); 
381
    }
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
    /// Resizes \ref vec to n
384
    inline void resize(int n) { 
385
386
      FUNCNAME("DOFVector<T>::resize()");
      TEST_EXIT_DBG((n >= 0))("Can't resize DOFVector to negative size\n"); 
387
      vec.resize(n);
388
    } 
389

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    /// Resizes \ref vec to n and inits new values with init
391
    inline void resize(int n, T init) { 
392
393
      FUNCNAME("DOFVector<T>::resize()");
      TEST_EXIT_DBG((n >= 0))("Can't resize DOFVector to negative size\n"); 
394
      vec.resize(n, init);
395
    } 
396

Thomas Witkowski's avatar
Thomas Witkowski committed
397
    /// Returns \ref vec[i]
398
399
    inline const T& operator[](DegreeOfFreedom i) const {
      FUNCNAME("DOFVector<T>::operator[]");
400
401
      TEST_EXIT_DBG((i>= 0) && (i < static_cast<int>(vec.size())))
	("Illegal vector index %d.\n",i);
402
      return vec[i];
403
    } 
404

Thomas Witkowski's avatar
Thomas Witkowski committed
405
    /// Returns \ref vec[i]
406
407
    inline T& operator[](DegreeOfFreedom i) {
      FUNCNAME("DOFVector<T>::operator[]");
408
409
      TEST_EXIT_DBG((i >= 0) && (i < static_cast<int>(vec.size())))
	("Illegal vector index %d.\n",i);
410
      return vec[i];
411
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
412
 
Thomas Witkowski's avatar
Thomas Witkowski committed
413
    /// Calculates Integral of this DOFVector
414
415
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
416
    /// Calculates L1 norm of this DOFVector
417
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
418
 
Thomas Witkowski's avatar
Thomas Witkowski committed
419
    /// Calculates L2 norm of this DOFVector
420
421
    inline double L2Norm(Quadrature* q = NULL) const {
      return sqrt(L2NormSquare());
422
    }
423

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

Thomas Witkowski's avatar
Thomas Witkowski committed
427
    /// Calculates H1 norm of this DOFVector
428
429
    inline double H1Norm(Quadrature* q = NULL) const {
      return sqrt(H1NormSquare());
430
    }
431

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

Thomas Witkowski's avatar
Thomas Witkowski committed
435
    /// Calculates euclidian norm of this DOFVector
436
437
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
438
    /// Returns square of the euclidian norm.
439
440
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
441
    /// Calculates l2 norm of this DOFVector
442
443
    inline double l2norm() const { 
      return nrm2();
444
    }
445

Thomas Witkowski's avatar
Thomas Witkowski committed
446
    /// Calculates the absolute sum of this DOFVector
447
448
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
449
    /// Calculates the l1 norm of this DOFVector
450
451
    inline double l1norm() const { 
      return asum();
452
    } 
453

Thomas Witkowski's avatar
Thomas Witkowski committed
454
    /// Calculates doublewell of this DOFVector
455
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
456
 
Thomas Witkowski's avatar
Thomas Witkowski committed
457
    /// Calculates the sum of this DOFVector
458
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
459
 
Thomas Witkowski's avatar
Thomas Witkowski committed
460
    /// Sets \ref vec[i] = val, i=0 , ... , size
461
462
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
463
    /// Assignment operator for setting complete vector to a certain value d
464
465
466
    inline DOFVector<T>& operator=(T d) {
      set(d); 
      return *this;
467
    } 
468

Thomas Witkowski's avatar
Thomas Witkowski committed
469
    /// Assignment operator between two vectors
470
    DOFVector<T>& operator=(const DOFVector<T>&); 
471

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

Thomas Witkowski's avatar
Thomas Witkowski committed
475
    /// Returns minimum of DOFVector
476
477
    T min() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
478
    /// Returns maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
479
480
    T max() const;

481
    /// Returns absolute maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
482
    T absMax() const;
483

484
    /// Used by interpol while mesh traversal
485
    static int interpolFct(ElInfo* elinfo);
486

487
    /// Prints \ref vec to stdout
488
489
    void print() const; 

490
    ///
491
492
    int calcMemoryUsage() const;

493
494
495
496
497
498
499
500
    /** \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);

501
    void serialize(std::ostream &out) {
502
503
504
      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));
505
    }
506

507
    void deserialize(std::istream &in) {
508
509
510
511
      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));
512
    }
513
514
515
516
517

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

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

518
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
519

520
521
522
523
524
525
526
#ifdef HAVE_PARALLEL_AMDIS
    /// Sets the petsc application ordering object to map dof indices.
    void useApplicationOrdering(AO *ao) {
      applicationOrdering = ao;
    }
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
527
  protected: 
Thomas Witkowski's avatar
Thomas Witkowski committed
528
    /// Used by Int while mesh traversal
529
530
    static int Int_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
531
    /// Used by L1Norm while mesh traversal
532
    static int L1Norm_fct(ElInfo* elinfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
533
 
Thomas Witkowski's avatar
Thomas Witkowski committed
534
    /// Used by L2Norm while mesh traversal
535
536
    static int L2Norm_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
537
    /// Used by H1Norm while mesh traversal
538
539
    static int H1Norm_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
540
    /// Used by DoubleWell while mesh traversal
541
542
    static int DoubleWell_fct(ElInfo* elinfo);

543
  protected: 
544
    /// Name of this DOFVector
545
    std::string name; 
546

547
    /// FiniteElemSpace of the vector
548
549
    const FiniteElemSpace *feSpace; 

550
    /// Data container
551
    std::vector<T> vec; 
552
 
553
    /// Specifies what operation should be performed after coarsening
554
555
    CoarsenOperation coarsenOperation;

556
    /// Used by \ref interpol
557
558
    AbstractFunction<T, WorldVector<double> > *interFct;

559
    /// Used for mesh traversal
560
561
    static DOFVector<T> *traverseVector;

562
563
564
565
566
#ifdef HAVE_PARALLEL_AMDIS
    /// Petsc application ordering to map dof indices.
    AO *applicationOrdering;
#endif

567
  protected:
568
    /// Used while calculating vector norms
569
570
    static FastQuadrature *quad_fast;

571
    /// Stores the last calculated vector norm
572
573
574
    static double norm;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
575

576
577
578
579
580
581
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
582
583
  inline double min(const DOFVector<double>& v) {
    return v.min();
584
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
585
586
587

  inline double max(const DOFVector<double>& v) {
    return v.max();
588
  }
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
    MEMORY_MANAGED(DOFVectorDOF);

    /** \brief
     * Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
     * as DOFContainer at DOFAdmin
     */
605
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
606
607
608
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
609
    }
610
  
611
    /// Deregisters itself at DOFAdmin.
612
613
    ~DOFVectorDOF() {
      feSpace->getAdmin()->removeDOFContainer(this);
614
    }
615
616
617
618
619
620
621

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVEctor<DegreeOfFreedom>::operator[]()
     */
    DegreeOfFreedom& operator[](DegreeOfFreedom i) {
      return DOFVector<DegreeOfFreedom>::operator[](i);
622
    }
623

624
625
626
627
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const {
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

628
    /// Implements DOFIndexedBase::getSize()
629
630
    int getSize() const {
      return DOFVector<DegreeOfFreedom>::getSize();
631
    }
632

633
    /// Implements DOFIndexedBase::resize()
634
635
636
637
638
639
640
641
642
643
    void resize(int size) {
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
644

645
646
647
  template<typename T>
  double norm(DOFVector<T> *vec) {
    return vec->nrm2();
648
  }
649
650
651
652

  template<typename T>
  double L2Norm(DOFVector<T> *vec) {
    return vec->L2Norm();
653
  }
654
655
656
657

  template<typename T>
  double H1Norm(DOFVector<T> *vec) {
    return vec->H1Norm();
658
  }
659
660
661
662

  template<typename T>
  void print(DOFVector<T> *vec) {
    vec->print();
663
  }
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709

  // 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>
710
711
712
  inline void scal(T a, DOFVector<T>& y) {
    y *= a;
  }
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728

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

729
730
731
732
733
734
735
736
737
  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);


738
  template<typename T>
739
  inline void set(DOFVector<T>& vec, T d) {
740
    vec.set(d);
741
  }
742
743

  template<typename T>
744
  inline void setValue(DOFVector<T>& vec, T d) {
745
    vec.set(d);
746
  }
747
748
749
750

  template<typename T>
  inline int size(DOFVector<T> *vec) {
    return vec->getUsedSize();
751
  }
752
753
754
755
756
757
758
759

  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_