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
#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

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

      void free(DOFVector<T> *vec) { 
297
	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
    /// Writes the data of the DOFVector to an output stream.
497
    void serialize(std::ostream &out) {
498
499
500
      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));
501
    }
502

503
    /// Reads data of an DOFVector from an input stream.
504
    void deserialize(std::istream &in) {
505
506
507
508
      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));
509
    }
510
511
512
513
514

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

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

515
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
516

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

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

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

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

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

540
  protected: 
541
    /// Name of this DOFVector
542
    std::string name; 
543

544
    /// FiniteElemSpace of the vector
545
546
    const FiniteElemSpace *feSpace; 

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

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

556
    /// Used for mesh traversal
557
558
    static DOFVector<T> *traverseVector;

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

564
  protected:
565
    /// Used while calculating vector norms
566
567
    static FastQuadrature *quad_fast;

568
    /// Stores the last calculated vector norm
569
570
571
    static double norm;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
572

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

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

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

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


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

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
639

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

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

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

  template<typename T>
  void print(DOFVector<T> *vec) {
    vec->print();
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
703
704

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

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

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


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

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

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

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

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_