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

/** \file DOFVector.h */

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

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

43
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
44
45
#include "petscao.h"
#endif
46
47
48
49
50
51
52
 
namespace AMDiS {

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
65
66
67
68
69
    /** \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;
70

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

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

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

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

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

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

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

    void addElementVector(T sign,
122
			  const ElementVector& elVec, 
123
			  const BoundaryType *bound,
124
			  ElInfo *elInfo,
125
			  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
    ElementVector elementVector;
212

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
  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)
279
      {}
280
281
282

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
283
      {}
284
285
286
287
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
288
289
290
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
291
292
293

      DOFVector<T> *create() { 
	return NEW DOFVector<T>(feSpace, ""); 
294
      }
295
296
297

      void free(DOFVector<T> *vec) { 
	DELETE vec; 
298
      }
299
300
301
302
303
304

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
305
    /// Empty constructor. No initialization!
306
307
308
309
    DOFVector() 
      : DOFVectorBase<T>(), 
	feSpace(NULL),
	coarsenOperation(NO_OPERATION)
310
    {}
311

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

Thomas Witkowski's avatar
Thomas Witkowski committed
315
    /// Initialization.
316
    void init(const FiniteElemSpace* f, std::string n);
317

Thomas Witkowski's avatar
Thomas Witkowski committed
318
    /// Copy Constructor
319
    DOFVector(const DOFVector& rhs) {
Thomas Witkowski's avatar
Thomas Witkowski committed
320
321
322
      *this = rhs;   
      name = rhs.name + "copy";
      if (feSpace && feSpace->getAdmin()) {
323
324
	(dynamic_cast<DOFAdmin*>(feSpace->getAdmin()))->addDOFIndexed(this);
      }
325
    }
326

Thomas Witkowski's avatar
Thomas Witkowski committed
327
    /// Destructor
328
329
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
330
    /// Returns iterator to the begin of \ref vec
331
    typename std::vector<T>::iterator begin() { 
332
      return vec.begin(); 
333
    }
334

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

    /// Sets \ref coarsenOperation to op
348
349
    inline void setCoarsenOperation(CoarsenOperation op) { 
      coarsenOperation = op; 
350
    }
351

Thomas Witkowski's avatar
Thomas Witkowski committed
352
    /// Returns \ref coarsenOperation
353
354
    inline CoarsenOperation getCoarsenOperation() { 
      return coarsenOperation; 
355
    }
356

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

Thomas Witkowski's avatar
Thomas Witkowski committed
360
    /// Interpolation after refinement.
361
    inline void refineInterpol(RCNeighbourList&, int) {}
362

Thomas Witkowski's avatar
Thomas Witkowski committed
363
    /// Returns \ref vec
364
    std::vector<T>& getVector() { 
Thomas Witkowski's avatar
Thomas Witkowski committed
365
      return vec;
366
    }
367

Thomas Witkowski's avatar
Thomas Witkowski committed
368
    /// Returns size of \ref vec
369
370
    inline int getSize() const { 
      return vec.size();
371
    } 
372

Thomas Witkowski's avatar
Thomas Witkowski committed
373
    /// Returns used size of the vector.
374
375
    inline int getUsedSize() const { 
      return feSpace->getAdmin()->getUsedSize(); 
376
    }
377

Thomas Witkowski's avatar
Thomas Witkowski committed
378
    /// Resizes \ref vec to n
379
    inline void resize(int n) { 
380
381
      FUNCNAME("DOFVector<T>::resize()");
      TEST_EXIT_DBG((n >= 0))("Can't resize DOFVector to negative size\n"); 
382
      vec.resize(n);
383
    } 
384

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
411
    /// Calculates L1 norm of this DOFVector
412
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
413
 
Thomas Witkowski's avatar
Thomas Witkowski committed
414
    /// Calculates L2 norm of this DOFVector
415
416
    inline double L2Norm(Quadrature* q = NULL) const {
      return sqrt(L2NormSquare());
417
    }
418

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

Thomas Witkowski's avatar
Thomas Witkowski committed
422
    /// Calculates H1 norm of this DOFVector
423
424
    inline double H1Norm(Quadrature* q = NULL) const {
      return sqrt(H1NormSquare());
425
    }
426

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

Thomas Witkowski's avatar
Thomas Witkowski committed
430
    /// Calculates euclidian norm of this DOFVector
431
432
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
433
    /// Returns square of the euclidian norm.
434
435
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
436
    /// Calculates l2 norm of this DOFVector
437
438
    inline double l2norm() const { 
      return nrm2();
439
    }
440

Thomas Witkowski's avatar
Thomas Witkowski committed
441
    /// Calculates the absolute sum of this DOFVector
442
443
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
444
    /// Calculates the l1 norm of this DOFVector
445
446
    inline double l1norm() const { 
      return asum();
447
    } 
448

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

Thomas Witkowski's avatar
Thomas Witkowski committed
458
    /// Assignment operator for setting complete vector to a certain value d
459
460
461
    inline DOFVector<T>& operator=(T d) {
      set(d); 
      return *this;
462
    } 
463

Thomas Witkowski's avatar
Thomas Witkowski committed
464
    /// Assignment operator between two vectors
465
    DOFVector<T>& operator=(const DOFVector<T>&); 
466

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

Thomas Witkowski's avatar
Thomas Witkowski committed
470
    /// Returns minimum of DOFVector
471
472
    T min() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
473
    /// Returns maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
474
475
    T max() const;

476
    /// Returns absolute maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
477
    T absMax() const;
478

479
    /// Used by interpol while mesh traversal
480
    static int interpolFct(ElInfo* elinfo);
481

482
    /// Prints \ref vec to stdout
483
484
    void print() const; 

485
    ///
486
487
    int calcMemoryUsage() const;

488
489
490
491
492
493
494
495
    /** \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);

496
    void serialize(std::ostream &out) {
497
498
499
      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));
500
    }
501

502
    void deserialize(std::istream &in) {
503
504
505
506
      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));
507
    }
508
509
510
511
512

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

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

513
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
514

515
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
516
517
518
519
520
521
    /// Sets the petsc application ordering object to map dof indices.
    void useApplicationOrdering(AO *ao) {
      applicationOrdering = ao;
    }
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
522
  protected: 
Thomas Witkowski's avatar
Thomas Witkowski committed
523
    /// Used by Int while mesh traversal
524
525
    static int Int_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
526
    /// Used by L1Norm while mesh traversal
527
    static int L1Norm_fct(ElInfo* elinfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
528
 
Thomas Witkowski's avatar
Thomas Witkowski committed
529
    /// Used by L2Norm while mesh traversal
530
531
    static int L2Norm_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
532
    /// Used by H1Norm while mesh traversal
533
534
    static int H1Norm_fct(ElInfo* elinfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
535
    /// Used by DoubleWell while mesh traversal
536
537
    static int DoubleWell_fct(ElInfo* elinfo);

538
  protected: 
539
    /// Name of this DOFVector
540
    std::string name; 
541

542
    /// FiniteElemSpace of the vector
543
544
    const FiniteElemSpace *feSpace; 

545
    /// Data container
546
    std::vector<T> vec; 
547
 
548
    /// Specifies what operation should be performed after coarsening
549
550
    CoarsenOperation coarsenOperation;

551
    /// Used by \ref interpol
552
553
    AbstractFunction<T, WorldVector<double> > *interFct;

554
    /// Used for mesh traversal
555
556
    static DOFVector<T> *traverseVector;

557
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
558
559
560
561
    /// Petsc application ordering to map dof indices.
    AO *applicationOrdering;
#endif

562
  protected:
563
    /// Used while calculating vector norms
564
565
    static FastQuadrature *quad_fast;

566
    /// Stores the last calculated vector norm
567
568
569
    static double norm;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
570

571
572
573
574
575
576
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
577
578
  inline double min(const DOFVector<double>& v) {
    return v.min();
579
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
580
581
582

  inline double max(const DOFVector<double>& v) {
    return v.max();
583
  }
584
585
586
587
588
589
590
591
592
593
594
595
596
597


  /** \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
     */
598
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
599
600
601
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
602
    }
603
  
604
    /// Deregisters itself at DOFAdmin.
605
606
    ~DOFVectorDOF() {
      feSpace->getAdmin()->removeDOFContainer(this);
607
    }
608
609
610
611
612
613
614

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

617
618
619
620
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const {
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

621
    /// Implements DOFIndexedBase::getSize()
622
623
    int getSize() const {
      return DOFVector<DegreeOfFreedom>::getSize();
624
    }
625

626
    /// Implements DOFIndexedBase::resize()
627
628
629
630
631
632
633
634
635
636
    void resize(int size) {
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
637

638
639
640
  template<typename T>
  double norm(DOFVector<T> *vec) {
    return vec->nrm2();
641
  }
642
643
644
645

  template<typename T>
  double L2Norm(DOFVector<T> *vec) {
    return vec->L2Norm();
646
  }
647
648
649
650

  template<typename T>
  double H1Norm(DOFVector<T> *vec) {
    return vec->H1Norm();
651
  }
652
653
654
655

  template<typename T>
  void print(DOFVector<T> *vec) {
    vec->print();
656
  }
657
658
659
660
661
662
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

  // 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>
703
704
705
  inline void scal(T a, DOFVector<T>& y) {
    y *= a;
  }
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721

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

722
723
724
725
726
727
728
729
730
  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);


731
  template<typename T>
732
  inline void set(DOFVector<T>& vec, T d) {
733
    vec.set(d);
734
  }
735
736

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

  template<typename T>
  inline int size(DOFVector<T> *vec) {
    return vec->getUsedSize();
744
  }
745
746
747
748
749
750
751
752

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

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_