DOFVector.h 21.2 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 getD2AtQPs. Thread safe.
273
    std::vector<DimMat<double>*> D2Phis;
Thomas Witkowski's avatar
Thomas Witkowski committed
274

Thomas Witkowski's avatar
Thomas Witkowski committed
275
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
276
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
277
278

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
279
    std::map<DegreeOfFreedom, bool> *rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
280
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
281
  };
282

Thomas Witkowski's avatar
Thomas Witkowski committed
283
  /// Specifies which operation should be done after coarsening
284
285
286
287
288
289
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
  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)
315
      {}
316
317
318

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
319
      {}
320
321
322
323
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
324
325
326
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
327

Thomas Witkowski's avatar
Thomas Witkowski committed
328
329
      DOFVector<T> *create() 
      { 
330
	return new DOFVector<T>(feSpace, ""); 
331
      }
332

Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
      void free(DOFVector<T> *vec) 
      { 
335
	delete vec; 
336
      }
337
338
339
340
341
342

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
343
    /// Empty constructor. No initialization!
344
345
346
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
347
    {}
348

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

Thomas Witkowski's avatar
Thomas Witkowski committed
352
    /// Initialization.
353
    void init(const FiniteElemSpace* f, std::string n);
354

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

      this->createTempData();
364
    }
365

Thomas Witkowski's avatar
Thomas Witkowski committed
366
    /// Destructor
367
368
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
369
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
    typename std::vector<T>::iterator begin() 
    { 
372
      return vec.begin(); 
373
    }
374

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

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
391
      coarsenOperation = op; 
392
    }
393

Thomas Witkowski's avatar
Thomas Witkowski committed
394
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
    inline CoarsenOperation getCoarsenOperation() 
    { 
397
      return coarsenOperation; 
398
    }
399

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

Thomas Witkowski's avatar
Thomas Witkowski committed
403
    /// Interpolation after refinement.
404
    inline void refineInterpol(RCNeighbourList&, int) {}
405

Thomas Witkowski's avatar
Thomas Witkowski committed
406
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
407
408
    std::vector<T>& getVector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
409
      return vec;
410
    }
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
    inline int getSize() const 
    { 
415
      return vec.size();
416
    } 
417

Thomas Witkowski's avatar
Thomas Witkowski committed
418
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
    inline int getUsedSize() const 
    { 
421
      return this->feSpace->getAdmin()->getUsedSize(); 
422
    }
423

Thomas Witkowski's avatar
Thomas Witkowski committed
424
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
    inline void resize(int n) 
    { 
427
      FUNCNAME("DOFVector<T>::resize()");
428
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
429
      vec.resize(n);
430
    } 
431

Thomas Witkowski's avatar
Thomas Witkowski committed
432
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
433
434
    inline void resize(int n, T init) 
    { 
435
      FUNCNAME("DOFVector<T>::resize()");
436
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
437
      vec.resize(n, init);
438
    } 
439

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

Thomas Witkowski's avatar
Thomas Witkowski committed
449
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
450
451
    inline T& operator[](DegreeOfFreedom i) 
    {
452
      FUNCNAME("DOFVector<T>::operator[]");
453
454

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

457
      return vec[i];
458
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
459
 
Thomas Witkowski's avatar
Thomas Witkowski committed
460
    /// Calculates Integral of this DOFVector
461
462
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
463
    /// Calculates L1 norm of this DOFVector
464
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
465
 
Thomas Witkowski's avatar
Thomas Witkowski committed
466
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
467
468
    inline double L2Norm(Quadrature* q = NULL) const 
    {
469
      return sqrt(L2NormSquare());
470
    }
471

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

Thomas Witkowski's avatar
Thomas Witkowski committed
475
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
    inline double H1Norm(Quadrature* q = NULL) const 
    {
478
      return sqrt(H1NormSquare());
479
    }
480

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

Thomas Witkowski's avatar
Thomas Witkowski committed
484
    /// Calculates euclidian norm of this DOFVector
485
486
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
487
    /// Returns square of the euclidian norm.
488
489
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
490
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
491
492
    inline double l2norm() const 
    { 
493
      return nrm2();
494
    }
495

Thomas Witkowski's avatar
Thomas Witkowski committed
496
    /// Calculates the absolute sum of this DOFVector
497
498
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
499
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
500
501
    inline double l1norm() const 
    { 
502
      return asum();
503
    } 
504

Thomas Witkowski's avatar
Thomas Witkowski committed
505
    /// Calculates doublewell of this DOFVector
506
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
507
 
Thomas Witkowski's avatar
Thomas Witkowski committed
508
    /// Calculates the sum of this DOFVector
509
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
510
 
Thomas Witkowski's avatar
Thomas Witkowski committed
511
    /// Sets \ref vec[i] = val, i=0 , ... , size
512
513
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
514
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
515
516
    inline DOFVector<T>& operator=(T d) 
    {
517
518
      set(d); 
      return *this;
519
    } 
520

Thomas Witkowski's avatar
Thomas Witkowski committed
521
    /// Assignment operator between two vectors
522
    DOFVector<T>& operator=(const DOFVector<T>&); 
523

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

527
    /// Returns minimum of DOFVector.
528
529
    T min() const; 

530
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
531
532
    T max() const;

533
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
534
    T absMax() const;
535

536
537
538
539
    /// Returns the average value of the DOFVector.
    T average() const;

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

542
    /// Prints \ref vec to stdout.
543
544
    void print() const; 

545
    ///
546
547
    int calcMemoryUsage() const;

548
549
550
551
552
553
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

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

556
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
557
558
    void serialize(std::ostream &out) 
    {
559
560
561
      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));
562
    }
563

564
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
565
566
    void deserialize(std::istream &in) 
    {
567
568
569
570
      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));
571
    }
572
573
574
575
576

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

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

577
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
578
579

  protected: 
580
    /// Data container
581
    std::vector<T> vec; 
582
 
583
    /// Specifies what operation should be performed after coarsening
584
585
    CoarsenOperation coarsenOperation;

586
    /// Used by \ref interpol
587
588
    AbstractFunction<T, WorldVector<double> > *interFct;

589
    /// Used for mesh traversal
590
591
592
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
593

594
595
596
597
598
599
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

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

Thomas Witkowski's avatar
Thomas Witkowski committed
600
601
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
602
    return v.min();
603
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
604

Thomas Witkowski's avatar
Thomas Witkowski committed
605
606
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
607
    return v.max();
608
  }
609
610
611
612
613
614
615
616
617
618
619
620
621
622


  /** \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
     */
623
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
624
625
626
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
627
    }
628
  
629
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
630
631
    ~DOFVectorDOF() 
    {
632
      feSpace->getAdmin()->removeDOFContainer(this);
633
    }
634
635
636

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
637
     * DOFVector<DegreeOfFreedom>::operator[]()
638
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
639
640
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
641
      return DOFVector<DegreeOfFreedom>::operator[](i);
642
    }
643

Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
646
647
648
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

649
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
650
651
    int getSize() const 
    {
652
      return DOFVector<DegreeOfFreedom>::getSize();
653
    }
654

655
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
656
657
    void resize(int size) 
    {
658
659
660
661
662
663
664
665
666
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
667

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

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
687
688
  void print(DOFVector<T> *vec) 
{
689
    vec->print();
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
723
724
725
726
727
728
729
730
731
732
733
734
735
736

  // 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
737
738
  inline void scal(T a, DOFVector<T>& y) 
  {
739
740
    y *= a;
  }
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756

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

757
758
759
760
761
762
763
764
765
  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);


766
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
767
768
  inline void set(DOFVector<T>& vec, T d) 
  {
769
    vec.set(d);
770
  }
771
772

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
773
774
  inline void setValue(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 int size(DOFVector<T> *vec) 
  {
781
    return vec->getUsedSize();
782
  }
783

Thomas Witkowski's avatar
Thomas Witkowski committed
784
785
786
787
788
789
790
791
792
793
  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());
  }

794
795
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815

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

819
820
821
822
  double integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<double, WorldVector<double> > *fct);


823
824
825
826
827
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_