DOFVector.h 21.3 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20
21
22
23
24
25

/** \file DOFVector.h */

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

30
#include "AMDiS_fwd.h"
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 "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
Thomas Witkowski's avatar
Thomas Witkowski committed
42
#include "BasisFunction.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
43
#include "FiniteElemSpace.h"
44

45
46
47
48
49
50
namespace AMDiS {

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

52
53
    DOFVectorBase() 
      : feSpace(NULL),
54
	elementVector(3),
Thomas Witkowski's avatar
Thomas Witkowski committed
55
56
        boundaryManager(NULL),
        nBasFcts(0)
57
58
59
60
61
    {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      rankDofs = NULL;
#endif
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    
63
    DOFVectorBase(const FiniteElemSpace *f, std::string n);
64

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

67
68
69
    /// Creates temporary used data structures.
    void createTempData();

Thomas Witkowski's avatar
Thomas Witkowski committed
70
71
72
73
    /** \brief 
     * For the given element, this function returns an array of all DOFs of this
     * DOFVector that are defined on this element.
     */
74
75
    virtual void getLocalVector(const Element *el, 
				mtl::dense_vector<T>& localVec) const;
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
78
79
80
    /** \brief
     * Evaluates the DOF vector at a set of quadrature points defined on the 
     * given element.
     */
81
82
83
84
    void getVecAtQPs(const ElInfo *elInfo, 
		     const Quadrature *quad,
		     const FastQuadrature *quadFast,
		     mtl::dense_vector<T>& vecAtQPs) const;
85

86
87
88
89
90
    void getVecAtQPs(const ElInfo *smallElInfo, 
		     const ElInfo *largeElInfo,
		     const Quadrature *quad,
		     const FastQuadrature *quadFast,
		     mtl::dense_vector<T>& vecAtQPs) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
91

92
93
    const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo, 
				      const Quadrature *quad,
94
				      const FastQuadrature *quadFast,
95
				      WorldVector<T> *grdAtQPs) const;
96

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

103
104
    const WorldMatrix<T> *getD2AtQPs(const ElInfo *elInfo, 
				     const Quadrature *quad,
105
				     const FastQuadrature *quadFast,
106
				     WorldMatrix<T> *d2AtQPs) const;
107

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

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

    void addElementVector(T sign,
129
			  const ElementVector& elVec, 
130
			  const BoundaryType *bound,
131
			  ElInfo *elInfo,
132
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
135
136
137
138
139

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

Thomas Witkowski's avatar
Thomas Witkowski committed
150
151
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
152
      return operatorFactor.begin();
153
    }
154

Thomas Witkowski's avatar
Thomas Witkowski committed
155
156
    inline std::vector<double*>::iterator getOperatorFactorEnd() 
    {
157
      return operatorFactor.end();
158
    }
159

Thomas Witkowski's avatar
Thomas Witkowski committed
160
161
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
162
      return operatorEstFactor.begin();
163
    }
164

Thomas Witkowski's avatar
Thomas Witkowski committed
165
166
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() 
    {
167
      return operatorEstFactor.end();
168
    }
169
170


Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
173
      return operators.begin();
174
    }
175

Thomas Witkowski's avatar
Thomas Witkowski committed
176
177
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
178
      return operators.end();
179
    }
180
181
182
183
184
185
186
187
188
189
190

    Flag getAssembleFlag();

    /** \brief
     * Evaluates \f[ u_h(x(\lambda)) = \sum_{i=0}^{m-1} vec[ind[i]] * 
     * \varphi^i(\lambda) \f] where \f$ \varphi^i \f$ is the i-th basis function,
     * \f$ x(\lambda) \f$ are the world coordinates of lambda and
     * \f$ m \f$ is the number of basis functions
     */
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

Thomas Witkowski's avatar
Thomas Witkowski committed
191
192
    inline std::vector<Operator*>& getOperators() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
193
      return operators; 
194
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
195

Thomas Witkowski's avatar
Thomas Witkowski committed
196
197
    inline std::vector<double*>& getOperatorFactor() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
198
      return operatorFactor; 
199
    }
200

Thomas Witkowski's avatar
Thomas Witkowski committed
201
202
    inline std::vector<double*>& getOperatorEstFactor() 
    { 
203
      return operatorEstFactor; 
204
    }
205

Thomas Witkowski's avatar
Thomas Witkowski committed
206
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
207
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
208
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
209
      return name; 
210
    } 
211

212
213
    inline void setName(std::string n)
    {
214
      name = n;
215
216
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
217
218
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
219
      return boundaryManager; 
220
    }
221

Thomas Witkowski's avatar
Thomas Witkowski committed
222
223
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
224
      boundaryManager = bm;
225
    }
226

Thomas Witkowski's avatar
Thomas Witkowski committed
227
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
228
    inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
229
    {
230
      rankDofs = &dofmap;
Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
233
234
    }

    inline bool isRankDof(DegreeOfFreedom dof)
    {
235
      TEST_EXIT_DBG(rankDofs)("No rank dofs set!\n");
236

237
      return (*rankDofs)[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
238
239
240
    }
#endif

241
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
242
    ///
243
244
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
245
    ///
246
    std::string name;
247

Thomas Witkowski's avatar
Thomas Witkowski committed
248
    ///
249
    ElementVector elementVector;
250

Thomas Witkowski's avatar
Thomas Witkowski committed
251
    ///
252
    std::vector<Operator*> operators;
253

Thomas Witkowski's avatar
Thomas Witkowski committed
254
    ///
255
    std::vector<double*> operatorFactor;
256

Thomas Witkowski's avatar
Thomas Witkowski committed
257
    ///
258
    std::vector<double*> operatorEstFactor;
259

Thomas Witkowski's avatar
Thomas Witkowski committed
260
    ///
261
262
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
263
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
264
    int nBasFcts;
265

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

Thomas Witkowski's avatar
Thomas Witkowski committed
269
    /// Are used to store temporary local values of an element. Thread safe.
270
    std::vector<mtl::dense_vector<T> > localVectors;
Thomas Witkowski's avatar
Thomas Witkowski committed
271

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
281
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
282
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
283
284

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
285
    std::map<DegreeOfFreedom, bool> *rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
286
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
287
  };
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
  /// Specifies which operation should be done after coarsening
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
  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)
321
      {}
322
323
324

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
325
      {}
326
327
328
329
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
330
331
332
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
333

Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
      DOFVector<T> *create() 
      { 
336
	return new DOFVector<T>(feSpace, ""); 
337
      }
338

Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
      void free(DOFVector<T> *vec) 
      { 
341
	delete vec; 
342
      }
343
344
345
346
347
348

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
349
    /// Empty constructor. No initialization!
350
351
352
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
353
    {}
354

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

Thomas Witkowski's avatar
Thomas Witkowski committed
358
    /// Initialization.
359
    void init(const FiniteElemSpace* f, std::string n);
360

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    /// Copy Constructor
Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
    DOFVector(const DOFVector& rhs) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
364
      *this = rhs;   
365
366
367
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
368
369

      this->createTempData();
370
    }
371

Thomas Witkowski's avatar
Thomas Witkowski committed
372
    /// Destructor
373
374
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
375
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
376
377
    typename std::vector<T>::iterator begin() 
    { 
378
      return vec.begin(); 
379
    }
380

Thomas Witkowski's avatar
Thomas Witkowski committed
381
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
    typename std::vector<T>::iterator end() 
    { 
384
      return vec.end(); 
385
    }
386
387
388
389
390
391
  
    /** \brief
     * Used by DOFAdmin to compress this DOFVector. Implementation of
     * DOFIndexedBase::compress()
     */
    virtual void compressDOFIndexed(int first, int last,
392
				    std::vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
393
394

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
397
      coarsenOperation = op; 
398
    }
399

Thomas Witkowski's avatar
Thomas Witkowski committed
400
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
401
402
    inline CoarsenOperation getCoarsenOperation() 
    { 
403
      return coarsenOperation; 
404
    }
405

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

Thomas Witkowski's avatar
Thomas Witkowski committed
409
    /// Interpolation after refinement.
410
    inline void refineInterpol(RCNeighbourList&, int) {}
411

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

Thomas Witkowski's avatar
Thomas Witkowski committed
418
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
    inline int getSize() const 
    { 
421
      return vec.size();
422
    } 
423

Thomas Witkowski's avatar
Thomas Witkowski committed
424
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
    inline int getUsedSize() const 
    { 
427
      return this->feSpace->getAdmin()->getUsedSize(); 
428
    }
429

Thomas Witkowski's avatar
Thomas Witkowski committed
430
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
    inline void resize(int n) 
    { 
433
      FUNCNAME("DOFVector<T>::resize()");
434
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
435
      vec.resize(n);
436
    } 
437

Thomas Witkowski's avatar
Thomas Witkowski committed
438
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
439
440
    inline void resize(int n, T init) 
    { 
441
      FUNCNAME("DOFVector<T>::resize()");
442
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
443
      vec.resize(n, init);
444
    } 
445

Thomas Witkowski's avatar
Thomas Witkowski committed
446
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
    inline const T& operator[](DegreeOfFreedom i) const 
    {
449
      FUNCNAME("DOFVector<T>::operator[]");
450
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
451
	("Illegal vector index %d.\n", i);
452
      return vec[i];
453
    } 
454

Thomas Witkowski's avatar
Thomas Witkowski committed
455
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
    inline T& operator[](DegreeOfFreedom i) 
    {
458
      FUNCNAME("DOFVector<T>::operator[]");
459
460

      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size())) 
461
 	("Illegal vector index %d.\n", i); 
462

463
      return vec[i];
464
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
465
 
Thomas Witkowski's avatar
Thomas Witkowski committed
466
    /// Calculates Integral of this DOFVector
467
468
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
469
    /// Calculates L1 norm of this DOFVector
470
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
471
 
Thomas Witkowski's avatar
Thomas Witkowski committed
472
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
473
474
    inline double L2Norm(Quadrature* q = NULL) const 
    {
475
      return sqrt(L2NormSquare());
476
    }
477

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

Thomas Witkowski's avatar
Thomas Witkowski committed
481
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
482
483
    inline double H1Norm(Quadrature* q = NULL) const 
    {
484
      return sqrt(H1NormSquare());
485
    }
486

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

Thomas Witkowski's avatar
Thomas Witkowski committed
490
    /// Calculates euclidian norm of this DOFVector
491
492
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
493
    /// Returns square of the euclidian norm.
494
495
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
496
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
497
498
    inline double l2norm() const 
    { 
499
      return nrm2();
500
    }
501

Thomas Witkowski's avatar
Thomas Witkowski committed
502
    /// Calculates the absolute sum of this DOFVector
503
504
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
505
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
506
507
    inline double l1norm() const 
    { 
508
      return asum();
509
    } 
510

Thomas Witkowski's avatar
Thomas Witkowski committed
511
    /// Calculates doublewell of this DOFVector
512
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
513
 
Thomas Witkowski's avatar
Thomas Witkowski committed
514
    /// Calculates the sum of this DOFVector
515
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
516
 
Thomas Witkowski's avatar
Thomas Witkowski committed
517
    /// Sets \ref vec[i] = val, i=0 , ... , size
518
519
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
520
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
    inline DOFVector<T>& operator=(T d) 
    {
523
524
      set(d); 
      return *this;
525
    } 
526

Thomas Witkowski's avatar
Thomas Witkowski committed
527
    /// Assignment operator between two vectors
528
    DOFVector<T>& operator=(const DOFVector<T>&); 
529

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

533
    /// Returns minimum of DOFVector.
534
535
    T min() const; 

536
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
537
538
    T max() const;

539
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
540
    T absMax() const;
541

542
543
544
545
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Used by interpol while mesh traversal.
546
    void interpolFct(ElInfo* elinfo);
547

548
    /// Prints \ref vec to stdout.
549
550
    void print() const; 

551
    ///
552
553
    int calcMemoryUsage() const;

554
555
556
557
558
559
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

560
    void interpol(DOFVector<T> *v, double factor = 1.0);
561

562
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
563
564
    void serialize(std::ostream &out) 
    {
565
566
567
      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));
568
    }
569

570
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
571
572
    void deserialize(std::istream &in) 
    {
573
574
575
576
      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));
577
    }
578
579
580
581
582

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

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

583
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
584
585

  protected: 
586
    /// Data container
587
    std::vector<T> vec; 
588
 
589
    /// Specifies what operation should be performed after coarsening
590
591
    CoarsenOperation coarsenOperation;

592
    /// Used by \ref interpol
593
594
    AbstractFunction<T, WorldVector<double> > *interFct;

595
    /// Used for mesh traversal
596
597
598
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
599

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

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

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

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


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
    /** \brief
     * Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
     * as DOFContainer at DOFAdmin
     */
629
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
630
631
632
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
633
    }
634
  
635
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
636
637
    ~DOFVectorDOF() 
    {
638
      feSpace->getAdmin()->removeDOFContainer(this);
639
    }
640
641
642

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
673

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

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
693
694
  void print(DOFVector<T> *vec) 
{
695
    vec->print();
696
  }
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742

  // point wise multiplication
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y);

  // multiplication with scalar
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal);

  // scalar product
  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y);

  // addition
  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y);

  // subtraction
  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y);

  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d);

  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v);

  template<typename T>
  const DOFVector<T>& operator+(const DOFVector<T>&v1 , const DOFVector<T>& v2);

  // y = a*x + y
  template<typename T>
  void axpy(double a,const DOFVector<T>& x, DOFVector<T>& y);

  // matrix vector product
  template<typename T>
  void mv(MatrixTranspose transpose, 
	  const DOFMatrix &a, 
	  const DOFVector<T> &x,
	  DOFVector<T> &result,
	  bool add = false); 

  template<typename T>
  void xpay(double a,const DOFVector<T>& x,DOFVector<T>& y);

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
  inline void scal(T a, DOFVector<T>& y) 
  {
745
746
    y *= a;
  }
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762

  template<typename T>
  inline const DOFVector<T>& mult(double scal,
				  const DOFVector<T>& v,
				  DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v1,
				 const DOFVector<T>& v2,
				 DOFVector<T>& result);

763
764
765
766
767
768
769
770
771
  template<typename T>
  inline const DOFVector<T>& mod(const DOFVector<T>& v,
				 DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& Tanh(const DOFVector<T>& v,
				  DOFVector<T>& result);


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
790
791
792
793
794
795
796
797
798
799
  template<typename T>
  inline void checkFeSpace(const FiniteElemSpace* feSpace, const std::vector<T>& vec)
  {
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       feSpace->getAdmin()->getUsedSize());
  }

800
801
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821

  /** \brief
   * Computes the integral of a function that includes two different DOFVectors. This
   * function works also for the case that the DOFVectors are defined on two different
   * meshes.
   */
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct);

  /// Computes the integral of a function with two DOFVectors defined on the same mesh.
  double intSingleMesh(const DOFVector<double> &vec1,
		       const DOFVector<double> &vec2,
		       BinaryAbstractFunction<double, double, double> *fct);

  /// Computes the integral of a function with two DOFVectors defined on different meshes.
  double intDualMesh(const DOFVector<double> &vec1,
		     const DOFVector<double> &vec2,
		     BinaryAbstractFunction<double, double, double> *fct);

Thomas Witkowski's avatar
Thomas Witkowski committed
822
823
824
  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct);

825
826
827
828
829
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_