DOFVector.h 20.8 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

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

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

83
84
85
86
87
    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
88

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

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

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

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

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

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

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

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

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

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

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


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

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

    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
188
189
    inline std::vector<Operator*>& getOperators() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
190
      return operators; 
191
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
192

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

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

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

209
210
    inline void setName(std::string n)
    {
211
      name = n;
212
213
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
216
      return boundaryManager; 
217
    }
218

Thomas Witkowski's avatar
Thomas Witkowski committed
219
220
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
221
      boundaryManager = bm;
222
    }
223

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

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

234
      return (*rankDofs)[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
237
    }
#endif

238
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
239
    ///
240
241
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
242
    ///
243
    std::string name;
244

Thomas Witkowski's avatar
Thomas Witkowski committed
245
    ///
246
    ElementVector elementVector;
247

Thomas Witkowski's avatar
Thomas Witkowski committed
248
    ///
249
    std::vector<Operator*> operators;
250

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
257
    ///
258
259
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
260
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
261
    int nBasFcts;
262

Thomas Witkowski's avatar
Thomas Witkowski committed
263
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
264
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
265
266

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
267
    std::map<DegreeOfFreedom, bool> *rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
268
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
269
  };
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
  /// Specifies which operation should be done after coarsening
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
299
300
301
302
  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)
303
      {}
304
305
306

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
307
      {}
308
309
310
311
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
312
313
314
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
315

Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
      DOFVector<T> *create() 
      { 
318
	return new DOFVector<T>(feSpace, ""); 
319
      }
320

Thomas Witkowski's avatar
Thomas Witkowski committed
321
322
      void free(DOFVector<T> *vec) 
      { 
323
	delete vec; 
324
      }
325
326
327
328
329
330

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
331
    /// Empty constructor. No initialization!
332
333
334
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
335
    {}
336

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

Thomas Witkowski's avatar
Thomas Witkowski committed
340
    /// Initialization.
341
    void init(const FiniteElemSpace* f, std::string n);
342

Thomas Witkowski's avatar
Thomas Witkowski committed
343
    /// Copy Constructor
Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
    DOFVector(const DOFVector& rhs) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
346
      *this = rhs;   
347
348
349
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
350
    }
351

Thomas Witkowski's avatar
Thomas Witkowski committed
352
    /// Destructor
353
354
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
    typename std::vector<T>::iterator begin() 
    { 
358
      return vec.begin(); 
359
    }
360

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
    typename std::vector<T>::iterator end() 
    { 
364
      return vec.end(); 
365
    }
366
367
368
369
370
371
  
    /** \brief
     * Used by DOFAdmin to compress this DOFVector. Implementation of
     * DOFIndexedBase::compress()
     */
    virtual void compressDOFIndexed(int first, int last,
372
				    std::vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
373
374

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
375
376
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
377
      coarsenOperation = op; 
378
    }
379

Thomas Witkowski's avatar
Thomas Witkowski committed
380
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
381
382
    inline CoarsenOperation getCoarsenOperation() 
    { 
383
      return coarsenOperation; 
384
    }
385

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

Thomas Witkowski's avatar
Thomas Witkowski committed
389
    /// Interpolation after refinement.
390
    inline void refineInterpol(RCNeighbourList&, int) {}
391

Thomas Witkowski's avatar
Thomas Witkowski committed
392
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
393
394
    std::vector<T>& getVector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
395
      return vec;
396
    }
397

Thomas Witkowski's avatar
Thomas Witkowski committed
398
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
399
400
    inline int getSize() const 
    { 
401
      return vec.size();
402
    } 
403

Thomas Witkowski's avatar
Thomas Witkowski committed
404
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
    inline int getUsedSize() const 
    { 
407
      return this->feSpace->getAdmin()->getUsedSize(); 
408
    }
409

Thomas Witkowski's avatar
Thomas Witkowski committed
410
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
411
412
    inline void resize(int n) 
    { 
413
      FUNCNAME("DOFVector<T>::resize()");
414
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
415
      vec.resize(n);
416
    } 
417

Thomas Witkowski's avatar
Thomas Witkowski committed
418
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
    inline void resize(int n, T init) 
    { 
421
      FUNCNAME("DOFVector<T>::resize()");
422
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
423
      vec.resize(n, init);
424
    } 
425

Thomas Witkowski's avatar
Thomas Witkowski committed
426
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
427
428
    inline const T& operator[](DegreeOfFreedom i) const 
    {
429
      FUNCNAME("DOFVector<T>::operator[]");
430
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
431
	("Illegal vector index %d.\n", i);
432
      return vec[i];
433
    } 
434

Thomas Witkowski's avatar
Thomas Witkowski committed
435
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
    inline T& operator[](DegreeOfFreedom i) 
    {
438
      FUNCNAME("DOFVector<T>::operator[]");
439
440

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

443
      return vec[i];
444
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
445
 
Thomas Witkowski's avatar
Thomas Witkowski committed
446
    /// Calculates Integral of this DOFVector
447
448
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
449
    /// Calculates L1 norm of this DOFVector
450
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
451
 
Thomas Witkowski's avatar
Thomas Witkowski committed
452
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
453
454
    inline double L2Norm(Quadrature* q = NULL) const 
    {
455
      return sqrt(L2NormSquare());
456
    }
457

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

Thomas Witkowski's avatar
Thomas Witkowski committed
461
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
462
463
    inline double H1Norm(Quadrature* q = NULL) const 
    {
464
      return sqrt(H1NormSquare());
465
    }
466

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

Thomas Witkowski's avatar
Thomas Witkowski committed
470
    /// Calculates euclidian norm of this DOFVector
471
472
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
473
    /// Returns square of the euclidian norm.
474
475
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
476
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
477
478
    inline double l2norm() const 
    { 
479
      return nrm2();
480
    }
481

Thomas Witkowski's avatar
Thomas Witkowski committed
482
    /// Calculates the absolute sum of this DOFVector
483
484
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
485
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
486
487
    inline double l1norm() const 
    { 
488
      return asum();
489
    } 
490

Thomas Witkowski's avatar
Thomas Witkowski committed
491
    /// Calculates doublewell of this DOFVector
492
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
493
 
Thomas Witkowski's avatar
Thomas Witkowski committed
494
    /// Calculates the sum of this DOFVector
495
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
496
 
Thomas Witkowski's avatar
Thomas Witkowski committed
497
    /// Sets \ref vec[i] = val, i=0 , ... , size
498
499
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
500
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
501
502
    inline DOFVector<T>& operator=(T d) 
    {
503
504
      set(d); 
      return *this;
505
    } 
506

Thomas Witkowski's avatar
Thomas Witkowski committed
507
    /// Assignment operator between two vectors
508
    DOFVector<T>& operator=(const DOFVector<T>&); 
509

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

513
    /// Returns minimum of DOFVector.
514
515
    T min() const; 

516
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
517
518
    T max() const;

519
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
520
    T absMax() const;
521

522
523
524
525
    /// Returns the average value of the DOFVector.
    T average() const;

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

528
    /// Prints \ref vec to stdout.
529
530
    void print() const; 

531
    ///
532
533
    int calcMemoryUsage() const;

534
535
536
537
538
539
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

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

542
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
543
544
    void serialize(std::ostream &out) 
    {
545
546
547
      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));
548
    }
549

550
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
    void deserialize(std::istream &in) 
    {
553
554
555
556
      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));
557
    }
558
559
560
561
562

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

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

563
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
564
565

  protected: 
566
    /// Data container
567
    std::vector<T> vec; 
568
 
569
    /// Specifies what operation should be performed after coarsening
570
571
    CoarsenOperation coarsenOperation;

572
    /// Used by \ref interpol
573
574
    AbstractFunction<T, WorldVector<double> > *interFct;

575
    /// Used for mesh traversal
576
577
578
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
579

580
581
582
583
584
585
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
586
587
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
588
    return v.min();
589
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
590

Thomas Witkowski's avatar
Thomas Witkowski committed
591
592
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
593
    return v.max();
594
  }
595
596
597
598
599
600
601
602
603
604
605
606
607
608


  /** \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
     */
609
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
610
611
612
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
613
    }
614
  
615
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
616
617
    ~DOFVectorDOF() 
    {
618
      feSpace->getAdmin()->removeDOFContainer(this);
619
    }
620
621
622

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
623
     * DOFVector<DegreeOfFreedom>::operator[]()
624
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
625
626
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
627
      return DOFVector<DegreeOfFreedom>::operator[](i);
628
    }
629

Thomas Witkowski's avatar
Thomas Witkowski committed
630
631
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
632
633
634
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

635
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
636
637
    int getSize() const 
    {
638
      return DOFVector<DegreeOfFreedom>::getSize();
639
    }
640

641
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
642
643
    void resize(int size) 
    {
644
645
646
647
648
649
650
651
652
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
653

654
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
655
656
  double norm(DOFVector<T> *vec) 
  {
657
    return vec->nrm2();
658
  }
659
660

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
661
662
  double L2Norm(DOFVector<T> *vec) 
  {
663
    return vec->L2Norm();
664
  }
665
666

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
667
668
  double H1Norm(DOFVector<T> *vec) 
  {
669
    return vec->H1Norm();
670
  }
671
672

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
673
674
  void print(DOFVector<T> *vec) 
{
675
    vec->print();
676
  }
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722

  // 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
723
724
  inline void scal(T a, DOFVector<T>& y) 
  {
725
726
    y *= a;
  }
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742

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

743
744
745
746
747
748
749
750
751
  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);


752
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
753
754
  inline void set(DOFVector<T>& vec, T d) 
  {
755
    vec.set(d);
756
  }
757
758

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
  inline void setValue(DOFVector<T>& vec, T d) 
  {
761
    vec.set(d);
762
  }
763
764

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
765
766
  inline int size(DOFVector<T> *vec) 
  {
767
    return vec->getUsedSize();
768
  }
769

Thomas Witkowski's avatar
Thomas Witkowski committed
770
771
772
773
774
775
776
777
778
779
  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());
  }

780
781
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801

  /** \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
802
803
804
  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct);

805
806
807
808
  double integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<double, WorldVector<double> > *fct);


809
810
811
812
813
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_