DOFVector.h 22.9 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
25
26
27
28
// ============================================================================
// ==                                                                        ==
// == 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 
 
// ===========================================================================
// ===== includes ============================================================
// ===========================================================================

29
30
31
32
#include <vector> 
#include <memory> 
#include <list> 

33
#include "AMDiS_fwd.h"
34
35
36
37
38
39
40
41
42
43
44
45
#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
46
#include "BasisFunction.h"
47
48
49
50

#ifdef HAVE_PARALLEL_AMDIS
#include "petscao.h"
#endif
51
52
53
54
55
56
57
58
59
60
61
 
namespace AMDiS {

  // ===========================================================================
  // ===== forward declarations ================================================
  // ===========================================================================

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

63
64
65
    DOFVectorBase() 
      : feSpace(NULL),
	elementVector(NULL),
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
        boundaryManager(NULL),
        nBasFcts(0)
68
    {}
Thomas Witkowski's avatar
Thomas Witkowski committed
69
    
70
    DOFVectorBase(const FiniteElemSpace *f, std::string n);
71

Thomas Witkowski's avatar
Thomas Witkowski committed
72
    virtual ~DOFVectorBase();
73

Thomas Witkowski's avatar
Thomas Witkowski committed
74
75
76
77
78
    /** \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;
79

Thomas Witkowski's avatar
Thomas Witkowski committed
80
81
82
83
    /** \brief
     * Evaluates the DOF vector at a set of quadrature points defined on the 
     * given element.
     */
84
85
    const T *getVecAtQPs(const ElInfo *elInfo, 
			 const Quadrature *quad,
86
			 const FastQuadrature *quadFast,
87
			 T *vecAtQPs) const;
88

Thomas Witkowski's avatar
Thomas Witkowski committed
89
90
91
92
93
94
    const T *getVecAtQPs(const ElInfo *smallElInfo, 
			 const ElInfo *largeElInfo,
			 const Quadrature *quad,
			 const FastQuadrature *quadFast,
			 T *vecAtQPs) const;

95
96
    const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo, 
				      const Quadrature *quad,
97
				      const FastQuadrature *quadFast,
98
				      WorldVector<T> *grdAtQPs) const;
99

Thomas Witkowski's avatar
Thomas Witkowski committed
100
101
102
103
104
105
    const WorldVector<T> *getGrdAtQPs(const ElInfo *smallElInfo, 
				      const ElInfo *largeElInfo,
				      const Quadrature *quad,
				      const FastQuadrature *quadFast,
				      WorldVector<T> *grdAtQPs) const;

106
107
    const WorldMatrix<T> *getD2AtQPs(const ElInfo *elInfo, 
				     const Quadrature *quad,
108
				     const FastQuadrature *quadFast,
109
				     WorldMatrix<T> *d2AtQPs) const;
110

Thomas Witkowski's avatar
Thomas Witkowski committed
111
    /// Returns the FE space the DOF vector is defined on.
112
    inline const FiniteElemSpace* getFESpace() const {
113
      return feSpace;
114
    }
115

Thomas Witkowski's avatar
Thomas Witkowski committed
116
117
118
119
120
121
122
123
124
125
126
127
128
    /** \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);
129
130
131
132
133

    void addElementVector(T sign,
			  const ElementVector &elVec, 
			  const BoundaryType *bound,
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
134
135
136
137
138
139
140

    /* \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();
141
142
143
144
145
146
147
148
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
149
    }
150

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

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

159
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
160
      return operatorEstFactor.begin();
161
    }
162

163
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
164
      return operatorEstFactor.end();
165
    }
166
167


168
    inline std::vector<Operator*>::iterator getOperatorsBegin() {
169
      return operators.begin();
170
    }
171

172
    inline std::vector<Operator*>::iterator getOperatorsEnd() {
173
      return operators.end();
174
    }
175
176
177
178
179
180
181
182
183
184
185

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

186
    inline std::vector<Operator*>& getOperators() { 
Thomas Witkowski's avatar
Thomas Witkowski committed
187
      return operators; 
188
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
189

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

194
    inline std::vector<double*>& getOperatorEstFactor() { 
195
      return operatorEstFactor; 
196
    }
197

Thomas Witkowski's avatar
Thomas Witkowski committed
198
    /// Returns \ref name
199
    inline const std::string& getName() const { 
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 { 
      return boundaryManager; 
205
    }
206
207
208

    inline void setBoundaryManager(BoundaryManager *bm) {
      boundaryManager = bm;
209
    }
210
211

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
212
    ///
213
214
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
215
    ///
216
    std::string name;
217

Thomas Witkowski's avatar
Thomas Witkowski committed
218
    ///
219
220
    ElementVector *elementVector;

Thomas Witkowski's avatar
Thomas Witkowski committed
221
    ///
222
    std::vector<Operator*> operators;
223

Thomas Witkowski's avatar
Thomas Witkowski committed
224
    ///
225
    std::vector<double*> operatorFactor;
226

Thomas Witkowski's avatar
Thomas Witkowski committed
227
    ///
228
    std::vector<double*> operatorEstFactor;
229

Thomas Witkowski's avatar
Thomas Witkowski committed
230
    ///
231
232
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
233
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
234
    int nBasFcts;
235

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
251
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
252
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
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
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298


  // ===========================================================================
  // ===== defs ================================================================
  // ===========================================================================

  /** \brief
   * Specifies which operation should be done after coarsening
   */
  typedef enum{
    NO_OPERATION = 0,   
    COARSE_RESTRICT = 1,
    COARSE_INTERPOL = 2 
  } CoarsenOperation;
 
  // ===========================================================================
  // ===== class DOFVector =====================================================
  // ===========================================================================

  /** \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)
299
      {}
300
301
302

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
303
      {}
304
305
306
307
308
309
    };

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

310
311
312
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
313
314
315

      DOFVector<T> *create() { 
	return NEW DOFVector<T>(feSpace, ""); 
316
      }
317
318
319

      void free(DOFVector<T> *vec) { 
	DELETE vec; 
320
      }
321
322
323
324
325
326
327
328
329
330
331

    private:
      FiniteElemSpace *feSpace;
    };

  public:
    /** \brief
     * Empty constructor. No initialization!
     */
    DOFVector() 
      : DOFVectorBase<T>(), 
332
	//elementVector(NULL), 
333
334
	feSpace(NULL),
	coarsenOperation(NO_OPERATION)
335
    {}
336
337
338
339

    /** \brief
     * Constructs a DOFVector with name n belonging to FiniteElemSpace f
     */
340
    DOFVector(const FiniteElemSpace* f, std::string n); 
341
342
343
344

    /** \brief
     * Initialization.
     */
345
    void init(const FiniteElemSpace* f, std::string n);
346
347
348
349
350

    /** \brief
     * Copy Constructor
     */
    DOFVector(const DOFVector& rhs) {
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353
      *this = rhs;   
      name = rhs.name + "copy";
      if (feSpace && feSpace->getAdmin()) {
354
355
	(dynamic_cast<DOFAdmin*>(feSpace->getAdmin()))->addDOFIndexed(this);
      }
356
    }
357
358
359
360
361
362
363
364
365

    /** \brief
     * Destructor
     */
    virtual ~DOFVector();

    /** \brief
     * Returns iterator to the begin of \ref vec
     */
366
    typename std::vector<T>::iterator begin() { 
367
      return vec.begin(); 
368
    }
369
370
371
372

    /** \brief
     * Returns iterator to the end of \ref vec
     */
373
    typename std::vector<T>::iterator end() { 
374
      return vec.end(); 
375
    }
376
377
378
379
380
381
  
    /** \brief
     * Used by DOFAdmin to compress this DOFVector. Implementation of
     * DOFIndexedBase::compress()
     */
    virtual void compressDOFIndexed(int first, int last,
382
				    std::vector<DegreeOfFreedom> &newDof);
383
384
385
386
387
    /** \brief
     * Sets \ref coarsenOperation to op
     */
    inline void setCoarsenOperation(CoarsenOperation op) { 
      coarsenOperation = op; 
388
    }
389
390
391
392

    /** \brief
     * Returns \ref coarsenOperation
     */
393
394
    inline CoarsenOperation getCoarsenOperation() { 
      return coarsenOperation; 
395
    }
396
397
398
399

    /** \brief
     * Restriction after coarsening. Implemented for DOFVector<double>
     */
400
    inline void coarseRestrict(RCNeighbourList&, int) {}
401
402
403
404

    /** \brief
     * Interpolation after refinement.
     */
405
    inline void refineInterpol(RCNeighbourList&, int) {}
406
407
408
409

    /** \brief
     * Returns \ref vec
     */
410
    std::vector<T>& getVector() { 
Thomas Witkowski's avatar
Thomas Witkowski committed
411
      return vec;
412
    }
413
414
415
416

    /** \brief
     * Returns size of \ref vec
     */
417
418
    inline int getSize() const { 
      return vec.size();
419
    } 
420
421
422
423

    /** \brief
     * Returns used size of the vector.
     */
424
425
    inline int getUsedSize() const { 
      return feSpace->getAdmin()->getUsedSize(); 
426
    }
427
428
429
430
431

    /** \brief
     * Resizes \ref vec to n
     */
    inline void resize(int n) { 
432
433
      FUNCNAME("DOFVector<T>::resize()");
      TEST_EXIT_DBG((n >= 0))("Can't resize DOFVector to negative size\n"); 
434
      vec.resize(n);
435
    } 
436
437
438
439
440

    /** \brief
     * Resizes \ref vec to n and inits new values with init
     */
    inline void resize(int n, T init) { 
441
442
      FUNCNAME("DOFVector<T>::resize()");
      TEST_EXIT_DBG((n >= 0))("Can't resize DOFVector to negative size\n"); 
443
      vec.resize(n, init);
444
    } 
445
446
447
448
449
450

    /** \brief
     * Returns \ref vec[i]
     */
    inline const T& operator[](DegreeOfFreedom i) const {
      FUNCNAME("DOFVector<T>::operator[]");
451
452
      TEST_EXIT_DBG((i>= 0) && (i < static_cast<int>(vec.size())))
	("Illegal vector index %d.\n",i);
453
      return vec[i];
454
    } 
455
456
457
458
459
460

    /** \brief
     * Returns \ref vec[i]
     */
    inline T& operator[](DegreeOfFreedom i) {
      FUNCNAME("DOFVector<T>::operator[]");
461
462
      TEST_EXIT_DBG((i >= 0) && (i < static_cast<int>(vec.size())))
	("Illegal vector index %d.\n",i);
463
      return vec[i];
464
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
465
 
466
467
468
469
470
471
472
473
474
    /** \brief
     * Calculates Integral of this DOFVector
     */
    double Int(Quadrature* q = NULL) const;

    /** \brief
     * Calculates L1 norm of this DOFVector
     */
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
475
 
476
477
478
479
480
    /** \brief
     * Calculates L2 norm of this DOFVector
     */
    inline double L2Norm(Quadrature* q = NULL) const {
      return sqrt(L2NormSquare());
481
    }
482
483
484
485
486
487
488
489
490
491
492

    /** \brief
     * Calculates square of L2 norm of this DOFVector
     */
    double L2NormSquare(Quadrature* q = NULL) const;

    /** \brief
     * Calculates H1 norm of this DOFVector
     */
    inline double H1Norm(Quadrature* q = NULL) const {
      return sqrt(H1NormSquare());
493
    }
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514

    /** \brief
     * Calculates square of H1 norm of this DOFVector
     */
    double H1NormSquare(Quadrature* q = NULL) const;  

    /** \brief
     * Calculates euclidian norm of this DOFVector
     */
    double nrm2() const; 

    /** \brief
     * Returns square of the euclidian norm.
     */
    double squareNrm2() const;

    /** \brief
     * Calculates l2 norm of this DOFVector
     */
    inline double l2norm() const { 
      return nrm2();
515
    }
516
517
518
519
520
521
522
523
524
525
526

    /** \brief
     * Calculates the absolute sum of this DOFVector
     */
    T asum() const; 

    /** \brief
     * Calculates the l1 norm of this DOFVector
     */
    inline double l1norm() const { 
      return asum();
527
    } 
528
529
530
531
532

    /** \brief
     * Calculates doublewell of this DOFVector
     */
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
533
 
534
535
536
537
    /** \brief
     * Calculates the sum of this DOFVector
     */
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
538
 
539
540
541
542
543
544
545
546
547
548
549
    /** \brief
     * Sets \ref vec[i] = val, i=0 , ... , size
     */
    void set(T val); 

    /** \brief
     * Assignment operator for setting complete vector to a certain value d
     */
    inline DOFVector<T>& operator=(T d) {
      set(d); 
      return *this;
550
    } 
551
552
553
554

    /** \brief
     * Assignment operator between two vectors
     */
555
    DOFVector<T>& operator=(const DOFVector<T>&); 
556
557
558
559
560
561
562
563
564
565
566
567
568
569

    /** \brief
     * vec[i] = v.vec[i]
     */
    void copy(const DOFVector<T>& v); 

    /** \brief
     * Returns minimum of DOFVector
     */
    T min() const; 

    /** \brief
     * Returns maximum of DOFVector
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
570
571
    T max() const;

572
    /// Returns absolute maximum of DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
573
    T absMax() const;
574

575
    /// Used by interpol while mesh traversal
576
    static int interpolFct(ElInfo* elinfo);
577

578
    /// Prints \ref vec to stdout
579
580
    void print() const; 

581
    ///
582
583
    int calcMemoryUsage() const;

584
585
586
587
588
589
590
591
    /** \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);

592
    void serialize(std::ostream &out) {
593
594
595
      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));
596
    }
597

598
    void deserialize(std::istream &in) {
599
600
601
602
      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));
603
    }
604
605
606
607
608

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

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

609
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
610

611
612
613
614
615
616
617
#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
618
  protected: 
619
620
621
622
623
624
625
626
627
    /** \brief
     * Used by Int while mesh traversal
     */
    static int Int_fct(ElInfo* elinfo);

    /** \brief
     * Used by L1Norm while mesh traversal
     */
    static int L1Norm_fct(ElInfo* elinfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
628
 
629
630
631
632
633
634
635
636
637
638
    /** \brief
     * Used by L2Norm while mesh traversal
     */
    static int L2Norm_fct(ElInfo* elinfo);

    /** \brief
     * Used by H1Norm while mesh traversal
     */
    static int H1Norm_fct(ElInfo* elinfo);

639
640
641
642
643
    /** \brief
     * Used by DoubleWell while mesh traversal
     */
    static int DoubleWell_fct(ElInfo* elinfo);

644
  protected: 
645
    /// Name of this DOFVector
646
    std::string name; 
647

648
    /// FiniteElemSpace of the vector
649
650
    const FiniteElemSpace *feSpace; 

651
    /// Data container
652
    std::vector<T> vec; 
653
 
654
    /// Specifies what operation should be performed after coarsening
655
656
    CoarsenOperation coarsenOperation;

657
    /// Used by \ref interpol
658
659
    AbstractFunction<T, WorldVector<double> > *interFct;

660
    /// Used for mesh traversal
661
662
    static DOFVector<T> *traverseVector;

663
664
665
666
667
#ifdef HAVE_PARALLEL_AMDIS
    /// Petsc application ordering to map dof indices.
    AO *applicationOrdering;
#endif

668
  protected:
669
    /// Used while calculating vector norms
670
671
    static FastQuadrature *quad_fast;

672
    /// Stores the last calculated vector norm
673
674
675
    static double norm;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
  // ===============================================================================

678
679
680
681
682
683
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
684
685
  inline double min(const DOFVector<double>& v) {
    return v.min();
686
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
687
688
689

  inline double max(const DOFVector<double>& v) {
    return v.max();
690
  }
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709

  // ===========================================================================
  // ===== class DOFVectorDOF ==================================================
  // ===========================================================================

  /** \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
     */
710
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
711
712
713
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
714
    }
715
  
716
    /// Deregisters itself at DOFAdmin.
717
718
    ~DOFVectorDOF() {
      feSpace->getAdmin()->removeDOFContainer(this);
719
    }
720
721
722
723
724
725
726

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

729
730
731
732
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const {
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

733
    /// Implements DOFIndexedBase::getSize()
734
735
    int getSize() const {
      return DOFVector<DegreeOfFreedom>::getSize();
736
    }
737

738
    /// Implements DOFIndexedBase::resize()
739
740
741
742
743
744
745
746
747
748
    void resize(int size) {
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
749
750
751
752
753


  // ===============================================================================


754
755
756
  template<typename T>
  double norm(DOFVector<T> *vec) {
    return vec->nrm2();
757
  }
758
759
760
761

  template<typename T>
  double L2Norm(DOFVector<T> *vec) {
    return vec->L2Norm();
762
  }
763
764
765
766

  template<typename T>
  double H1Norm(DOFVector<T> *vec) {
    return vec->H1Norm();
767
  }
768
769
770
771

  template<typename T>
  void print(DOFVector<T> *vec) {
    vec->print();
772
  }
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825

  // 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>
826
827
828
  inline void scal(T a, DOFVector<T>& y) {
    y *= a;
  }
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844

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

845
846
847
848
849
850
851
852
853
  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);


854
  template<typename T>
855
  inline void set(DOFVector<T>& vec, T d) {
856
    vec.set(d);
857
  }
858
859

  template<typename T>
860
  inline void setValue(DOFVector<T>& vec, T d) {
861
    vec.set(d);
862
  }
863
864
865
866

  template<typename T>
  inline int size(DOFVector<T> *vec) {
    return vec->getUsedSize();
867
  }
868
869
870
871
872
873
874
875

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

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_