DOFVector.h 26.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
#include "SurfaceQuadrature.h"
45

46
47
48
49
50
51
namespace AMDiS {

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

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

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

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
90
91
92
93
    /** \brief
     * Evaluates the gradient of a DOF vector at a set of quadrature points defined on the
     * given element.
     */
94
95
96
97
98
99
100
101
102
103
104
    void getGrdAtQPs( const ElInfo *elInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;

    void getGrdAtQPs( const ElInfo *smallElInfo,
		      const ElInfo *largeElInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;

Praetorius, Simon's avatar
Praetorius, Simon committed
105
106
107
108
    /** \brief
     * Evaluates the comp'th component of the derivative of a DOF vector at a set of quadrature points defined on the
     * given element.
     */
109
110
111
112
113
114
115
116
117
118
119
120
121
    void getDerivativeAtQPs( const ElInfo *elInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      int comp,
		      mtl::dense_vector<T> &derivativeAtQPs) const;

    void getDerivativeAtQPs( const ElInfo *smallElInfo,
		      const ElInfo *largeElInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      int comp,
		      mtl::dense_vector<T> &derivativeAtQPs) const;

Praetorius, Simon's avatar
Praetorius, Simon committed
122
123
124
125
    /** \brief
     * Evaluates the jacobian of a DOF vector at a set of quadrature points defined on the
     * given element.
     */
126
127
128
129
    void getD2AtQPs(const ElInfo *elInfo,
		    const Quadrature *quad,
		    const FastQuadrature *quadFast,
		    mtl::dense_vector<typename D2Type<T>::type> &d2AtQPs) const;
130

Thomas Witkowski's avatar
Thomas Witkowski committed
131
    /// Returns the FE space the DOF vector is defined on.
132
    inline const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
133
    {
134
      return feSpace;
135
    }
136

Praetorius, Simon's avatar
Praetorius, Simon committed
137
138
139
140
    /** \brief
     * Assembles the element vector for the given ellement and adds the
     * element matrix to the current DOF vector.
     */ 
Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
143
144
145
146
147
148
149
    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);
150
151

    void addElementVector(T sign,
152
			  const ElementVector& elVec, 
153
			  const BoundaryType *bound,
154
			  ElInfo *elInfo,
155
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
156

Praetorius, Simon's avatar
Praetorius, Simon committed
157
158
159
160
161
    /* \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.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
162
    void finishAssembling();
163
164
165
166
167
168
169
170
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
171
    }
172

Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
175
      return operatorFactor.begin();
176
    }
177

Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
    inline std::vector<double*>::iterator getOperatorFactorEnd() 
    {
180
      return operatorFactor.end();
181
    }
182

Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
185
      return operatorEstFactor.begin();
186
    }
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188
189
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() 
    {
190
      return operatorEstFactor.end();
191
    }
192
193


Thomas Witkowski's avatar
Thomas Witkowski committed
194
195
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
196
      return operators.begin();
197
    }
198

Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
201
      return operators.end();
202
    }
203
204
205
206
207
208
209
210
211
212
213

    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
214
215
    inline std::vector<Operator*>& getOperators() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
216
      return operators; 
217
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
218

Thomas Witkowski's avatar
Thomas Witkowski committed
219
220
    inline std::vector<double*>& getOperatorFactor() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
221
      return operatorFactor; 
222
    }
223

Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
    inline std::vector<double*>& getOperatorEstFactor() 
    { 
226
      return operatorEstFactor; 
227
    }
228

Thomas Witkowski's avatar
Thomas Witkowski committed
229
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
230
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
231
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
232
      return name; 
233
    } 
234

235
236
    inline void setName(std::string n)
    {
237
      name = n;
238
239
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
240
241
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
242
      return boundaryManager; 
243
    }
244

Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
247
      boundaryManager = bm;
248
    }
249

Thomas Witkowski's avatar
Thomas Witkowski committed
250
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
251
    inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
252
    {
253
      rankDofs = &dofmap;
Thomas Witkowski's avatar
Thomas Witkowski committed
254
255
256
257
    }

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

260
      return (*rankDofs)[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
263
    }
#endif

264
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
265
    ///
266
267
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    ///
269
    std::string name;
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    ///
272
    ElementVector elementVector;
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
    ///
275
    std::vector<Operator*> operators;
276

Thomas Witkowski's avatar
Thomas Witkowski committed
277
    ///
278
    std::vector<double*> operatorFactor;
279

Thomas Witkowski's avatar
Thomas Witkowski committed
280
    ///
281
    std::vector<double*> operatorEstFactor;
282

Thomas Witkowski's avatar
Thomas Witkowski committed
283
    ///
284
285
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
286
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
287
    int nBasFcts;
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
290
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
291
292

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
293
    std::map<DegreeOfFreedom, bool> *rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
294
#endif
295

Thomas Witkowski's avatar
Thomas Witkowski committed
296
  };
297

Thomas Witkowski's avatar
Thomas Witkowski committed
298
  /// Specifies which operation should be done after coarsening
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
  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)
330
      {}
331
332
333

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
334
      {}
335
336
337
338
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
339
340
341
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
342

Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
      DOFVector<T> *create() 
      { 
345
	return new DOFVector<T>(feSpace, ""); 
346
      }
347

Thomas Witkowski's avatar
Thomas Witkowski committed
348
349
      void free(DOFVector<T> *vec) 
      { 
350
	delete vec; 
351
      }
352
353
354
355
356
357

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
358
    /// Empty constructor. No initialization!
359
360
361
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
362
    {}
363

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

Thomas Witkowski's avatar
Thomas Witkowski committed
367
    /// Initialization.
368
    void init(const FiniteElemSpace* f, std::string n);
369

Thomas Witkowski's avatar
Thomas Witkowski committed
370
    /// Copy Constructor
Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
    DOFVector(const DOFVector& rhs) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
373
      *this = rhs;   
374
375
376
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
377
    }
378

Thomas Witkowski's avatar
Thomas Witkowski committed
379
    /// Destructor
380
381
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
382
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
383
384
    typename std::vector<T>::iterator begin() 
    { 
385
      return vec.begin(); 
386
    }
387

Thomas Witkowski's avatar
Thomas Witkowski committed
388
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
    typename std::vector<T>::iterator end() 
    { 
391
      return vec.end(); 
392
    }
393
  
Praetorius, Simon's avatar
Praetorius, Simon committed
394
395
396
397
    /** \brief
     * Used by DOFAdmin to compress this DOFVector. Implementation of
     * DOFIndexedBase::compress()
     */
398
    virtual void compressDOFIndexed(int first, int last,
399
				    std::vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
402
403
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
404
      coarsenOperation = op; 
405
    }
406

Thomas Witkowski's avatar
Thomas Witkowski committed
407
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
    inline CoarsenOperation getCoarsenOperation() 
    { 
410
      return coarsenOperation; 
411
    }
412

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

416
    /// Interpolation after refinement. Implemented for DOFVector<double>
417
    inline void refineInterpol(RCNeighbourList&, int) {}
418

Thomas Witkowski's avatar
Thomas Witkowski committed
419
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
420
421
    std::vector<T>& getVector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
422
      return vec;
423
    }
424

Thomas Witkowski's avatar
Thomas Witkowski committed
425
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
426
427
    inline int getSize() const 
    { 
428
      return vec.size();
429
    } 
430

Thomas Witkowski's avatar
Thomas Witkowski committed
431
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
    inline int getUsedSize() const 
    { 
434
      return this->feSpace->getAdmin()->getUsedSize(); 
435
    }
436

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

Thomas Witkowski's avatar
Thomas Witkowski committed
445
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
    inline void resize(int n, T init) 
    { 
448
      FUNCNAME("DOFVector<T>::resize()");
449
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
450
      vec.resize(n, init);
451
    } 
452

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

Thomas Witkowski's avatar
Thomas Witkowski committed
462
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
    inline T& operator[](DegreeOfFreedom i) 
    {
465
      FUNCNAME("DOFVector<T>::operator[]");
466
467

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

470
      return vec[i];
471
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
472
 
Thomas Witkowski's avatar
Thomas Witkowski committed
473
    /// Calculates Integral of this DOFVector
474
475
    double Int(Quadrature* q = NULL) const;

Praetorius, Simon's avatar
Praetorius, Simon committed
476
477
478
479
    /** \brief
     * Calculates Integral of this DOFVector over parts of the domain
     * boundary, indicated by boundaryType. Implemented for DOFVector<double>
     **/
480
481
482
483
484
485
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
486
487
488
489
490
491

    /** \brief
     * Calculates Integral of this DOFVector times normal vector over parts 
     * of the domain boundary, indicated by boundaryType. Implemented for 
     * DOFVector<WorldVector<double> >
     **/
492
493
494
495
496
497
498
    double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = NULL) const
    {
      FUNCNAME("DOFVector::IntOnBoundaryNormal())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
499
    /// Calculates L1 norm of this DOFVector
500
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
501
 
Thomas Witkowski's avatar
Thomas Witkowski committed
502
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
503
504
    inline double L2Norm(Quadrature* q = NULL) const 
    {
505
      return sqrt(L2NormSquare());
506
    }
507

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

Thomas Witkowski's avatar
Thomas Witkowski committed
511
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
    inline double H1Norm(Quadrature* q = NULL) const 
    {
514
      return sqrt(H1NormSquare());
515
    }
516

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

Thomas Witkowski's avatar
Thomas Witkowski committed
520
    /// Calculates euclidian norm of this DOFVector
521
522
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
523
    /// Returns square of the euclidian norm.
524
525
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
526
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
527
528
    inline double l2norm() const 
    { 
529
      return nrm2();
530
    }
531

Thomas Witkowski's avatar
Thomas Witkowski committed
532
    /// Calculates the absolute sum of this DOFVector
533
534
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
535
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
536
537
    inline double l1norm() const 
    { 
538
      return asum();
539
    } 
540

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

Thomas Witkowski's avatar
Thomas Witkowski committed
550
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
    inline DOFVector<T>& operator=(T d) 
    {
553
554
      set(d); 
      return *this;
555
    } 
556

Thomas Witkowski's avatar
Thomas Witkowski committed
557
    /// Assignment operator between two vectors
558
    DOFVector<T>& operator=(const DOFVector<T>&); 
559

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

563
    /// Returns minimum of DOFVector.
564
565
    T min() const; 

566
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
567
568
    T max() const;

569
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
570
    T absMax() const;
571

572
573
574
575
    /// Returns the average value of the DOFVector.
    T average() const;

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
584
585
586
587
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
588
589
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

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

592

Praetorius, Simon's avatar
Praetorius, Simon committed
593
594
595
596
597
598
    /** eval DOFVector at given point p. If oldElInfo != NULL the search for the element, where p is inside,
     *  starts from oldElInfo. 
     *  implemented for: double, WorldVector< double >
     */
    inline const T& evalAtPoint(
      WorldVector<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const 
599
600
601
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
Praetorius, Simon's avatar
Praetorius, Simon committed
602
      return *value;
603
    }
604

Praetorius, Simon's avatar
Praetorius, Simon committed
605
606
607
608
    /** determine the DegreeOfFreedom that has coords with minimal euclidean distance to WorldVector p
     *  return true if DOF is found, and false otherwise
     */
    const bool getDOFidxAtPoint(WorldVector<double> &p, DegreeOfFreedom &idx, ElInfo *oldElInfo = NULL, bool useOldElInfo = false) const;
609

610
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
611
612
    void serialize(std::ostream &out) 
    {
613
614
615
      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));
616
    }
617

618
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
619
620
    void deserialize(std::istream &in) 
    {
621
622
623
624
      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));
625
    }
626

Praetorius, Simon's avatar
Praetorius, Simon committed
627
//     DOFVector<WorldVector<T> > *getGradient(DOFVector<WorldVector<T> >*) const;
628
629
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
630
631
632

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

Praetorius, Simon's avatar
Praetorius, Simon committed
633
//     DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
634
635
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
636
637

  protected: 
638

639
    /// Data container
640
    std::vector<T> vec; 
641
 
642
    /// Specifies what operation should be performed after coarsening
643
644
    CoarsenOperation coarsenOperation;

645
    /// Used by \ref interpol
646
647
    AbstractFunction<T, WorldVector<double> > *interFct;

648
    /// Used for mesh traversal
649
650
651
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
652

653
  template<>
654
655
656
657
658
659
660
661
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

  template<>
  double DOFVector<WorldVector<double> >::IntOnBoundaryNormal(
    BoundaryType boundaryType, Quadrature* q) const;

  template<>
Praetorius, Simon's avatar
Praetorius, Simon committed
662
663
  const double& DOFVector<double>::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
664
665

  template<>
Praetorius, Simon's avatar
Praetorius, Simon committed
666
667
  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
668

669
  template<>
670
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
671
672

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

Thomas Witkowski's avatar
Thomas Witkowski committed
675
676
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
677
    return v.min();
678
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
679

Thomas Witkowski's avatar
Thomas Witkowski committed
680
681
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
682
    return v.max();
683
  }
684
685
686
687
688
689
690
691
692
693


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Praetorius, Simon's avatar
Praetorius, Simon committed
694
695
696
697
    /** \brief
     * Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
     * as DOFContainer at DOFAdmin
     */
698
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
699
700
701
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
702
    }
703
  
704
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
705
706
    ~DOFVectorDOF() 
    {
707
      feSpace->getAdmin()->removeDOFContainer(this);
708
    }
709

Praetorius, Simon's avatar
Praetorius, Simon committed
710
711
712
713
    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVector<DegreeOfFreedom>::operator[]()
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
716
      return DOFVector<DegreeOfFreedom>::operator[](i);
717
    }
718

Thomas Witkowski's avatar
Thomas Witkowski committed
719
720
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
721
722
723
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

724
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
725
726
    int getSize() const 
    {
727
      return DOFVector<DegreeOfFreedom>::getSize();
728
    }
729

730
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
731
732
    void resize(int size) 
    {
733
734
735
736
737
738
739
740
741
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
742

743
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
  double norm(DOFVector<T> *vec) 
  {
746
    return vec->nrm2();
747
  }
748
749

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
  double L2Norm(DOFVector<T> *vec) 
  {
752
    return vec->L2Norm();
753
  }
754
755

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
756
757
  double H1Norm(DOFVector<T> *vec) 
  {
758
    return vec->H1Norm();
759
  }
760
761

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
762
  void print(DOFVector<T> *vec) 
763
  {
764
    vec->print();
765
  }
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811

  // 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
812
813
  inline void scal(T a, DOFVector<T>& y) 
  {
814
815
    y *= a;
  }
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831

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

832
833
834
835
836
837
838
839
840
  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);


841
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
842
843
  inline void set(DOFVector<T>& vec, T d) 
  {
844
    vec.set(d);
845
  }
846
847

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
848
849
  inline void setValue(DOFVector<T>& vec, T d) 
  {
850
    vec.set(d);