DOFVector.h 27.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
23
24
25
26
27

/** \file DOFVector.h */

#ifndef AMDIS_DOFVECTOR_H 
#define AMDIS_DOFVECTOR_H 
 
28
29
#include <vector> 
#include <memory> 
30
#include <map> 
31

32
#include "AMDiS_fwd.h"
33
34
35
36
37
38
39
40
41
42
43
#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
44
#include "BasisFunction.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
45
#include "FiniteElemSpace.h"
46
#include "SurfaceQuadrature.h"
47
#include "Traits.h"
48

49
50
51
52
53
54
namespace AMDiS {

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

56
    DOFVectorBase() 
57
      : feSpace(NULL),
58
	elementVector(3),
59
        boundaryManager(NULL),
Thomas Witkowski's avatar
Thomas Witkowski committed
60
        nBasFcts(0)
61
    {}
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
    /// For the given element, this function returns an array of all DOFs of 
    /// this DOFVector that are defined on this element.
69
70
    virtual void getLocalVector(const Element *el, 
				mtl::dense_vector<T>& localVec) const;
71

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

79
80
81
82
83
    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
84

Thomas Witkowski's avatar
Thomas Witkowski committed
85
86
    /// Evaluates the gradient of a DOF vector at a set of quadrature points
    /// defined on the given element.
87
88
89
90
91
92
93
94
95
96
97
    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;

Thomas Witkowski's avatar
Thomas Witkowski committed
98
99
    /// Evaluates the comp'th component of the derivative of a DOF vector at a
    /// set of quadrature points defined on the given element.
100
101
102
103
104
105
106
107
108
109
110
111
112
    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;

Thomas Witkowski's avatar
Thomas Witkowski committed
113
114
    /// Evaluates the jacobian of a DOF vector at a set of quadrature points
    ///  defined on the given element.
115
116
117
118
    void getD2AtQPs(const ElInfo *elInfo,
		    const Quadrature *quad,
		    const FastQuadrature *quadFast,
		    mtl::dense_vector<typename D2Type<T>::type> &d2AtQPs) const;
119

Thomas Witkowski's avatar
Thomas Witkowski committed
120
    /// Returns the FE space the DOF vector is defined on.
121
    inline const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
122
    {
123
      return feSpace;
124
    }
125

Thomas Witkowski's avatar
Thomas Witkowski committed
126
127
    /// 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
128
129
    void assemble(T factor, ElInfo *elInfo,			    
		  const BoundaryType *bound, 
130
		  Operator *op = NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
131
132
133
134
135

    void assemble2(T factor, 
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
		   const BoundaryType *bound, 
136
		   Operator *op = NULL);
137
138

    void addElementVector(T sign,
139
			  const ElementVector& elVec, 
140
			  const BoundaryType *bound,
141
			  ElInfo *elInfo,
142
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
143

144
145
146
    ///
    void assembleOperator(Operator &op);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
147
148
149
    /// 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
150
    void finishAssembling();
151
152
 
    inline void addOperator(Operator* op, 
153
154
			    double *factor = NULL,
			    double *estFactor = NULL) 
155
156
157
158
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
159
    }
160

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

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

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

176
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
177
    {
178
      return operatorEstFactor.end();
179
    }
180
181


182
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
183
    {
184
      return operators.begin();
185
    }
186

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

    Flag getAssembleFlag();

194
195
196
197
    /// 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
198
199
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

200
    inline std::vector<Operator*>& getOperators() 
Thomas Witkowski's avatar
Thomas Witkowski committed
201
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
202
      return operators; 
203
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
204

205
    inline std::vector<double*>& getOperatorFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
206
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
207
      return operatorFactor; 
208
    }
209

210
    inline std::vector<double*>& getOperatorEstFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
211
    { 
212
      return operatorEstFactor; 
213
    }
214

Thomas Witkowski's avatar
Thomas Witkowski committed
215
    /// Returns \ref name
216
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
217
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
218
      return name; 
219
    } 
220

221
    inline void setName(std::string n)
222
    {
223
      name = n;
224
225
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
226
227
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
228
      return boundaryManager; 
229
    }
230

Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
233
      boundaryManager = bm;
234
    }
235

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
236
237
238
239
240
241
    inline void setDirichletDofValue(DegreeOfFreedom dof,
				     T value)
    {
      dirichletDofValues[dof] = value;
    }

242
    std::map<DegreeOfFreedom, T>& getDirichletValues()
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
243
244
245
246
    {
      return dirichletDofValues;
    }

247
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
248
    ///
249
250
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
251
    ///
252
    std::string name;
253

Thomas Witkowski's avatar
Thomas Witkowski committed
254
    ///
255
    ElementVector elementVector;
256

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

Thomas Witkowski's avatar
Thomas Witkowski committed
260
    ///
261
    std::vector<double*> operatorFactor;
262

Thomas Witkowski's avatar
Thomas Witkowski committed
263
    ///
264
    std::vector<double*> operatorEstFactor;
265

Thomas Witkowski's avatar
Thomas Witkowski committed
266
    ///
267
268
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
269
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
270
    int nBasFcts;
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
273
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
274

275
    std::map<DegreeOfFreedom, T> dirichletDofValues;
Thomas Witkowski's avatar
Thomas Witkowski committed
276
  };
277

278

Thomas Witkowski's avatar
Thomas Witkowski committed
279
  /// Specifies which operation should be done after coarsening
280
281
282
  typedef enum{
    NO_OPERATION = 0,   
    COARSE_RESTRICT = 1,
283
284
285
    COARSE_INTERPOL = 2, 
    REFINE_INTERPOL = 4
  } RefineCoarsenOperation;
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
 

  /** \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
  {  
303
304
305
306
307
308
  public:
    typedef typename DOFVectorBase<T>::value_type	value_type;
    typedef typename DOFVectorBase<T>::size_type		size_type;
    typedef typename DOFVectorBase<T>::reference		reference;
    typedef typename DOFVectorBase<T>::const_reference	const_reference;
  
309
310
311
312
313
314
315
316
317
  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)
318
      {}
319
320
321

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
322
      {}
323
324
325
326
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
327
328
329
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
330

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

Thomas Witkowski's avatar
Thomas Witkowski committed
352
    /// Constructs a DOFVector with name n belonging to FiniteElemSpace f
353
    DOFVector(const FiniteElemSpace* f, std::string n, bool addToSynch = false); 
354

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    /// Initialization.
356
    void init(const FiniteElemSpace* f, std::string n, bool addToSynch = false);
357

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

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

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

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

    /// Returns const_iterator to the begin of \ref vec
383
    typename std::vector<T>::const_iterator begin() const
Praetorius, Simon's avatar
Praetorius, Simon committed
384
385
386
387
388
    {
      return vec.begin();
    }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
405
    /// Returns \ref coarsenOperation
406
    inline RefineCoarsenOperation getCoarsenOperation() 
Thomas Witkowski's avatar
Thomas Witkowski committed
407
    { 
408
      return coarsenOperation; 
409
    }
410

411
412
413
414
415
416
417
418
419
420
421
422
    /// Sets \ref refineOperation to op
    inline void setRefineOperation(RefineCoarsenOperation op) 
    { 
      refineOperation = op; 
    }

    /// Returns \ref refineOperation
    inline RefineCoarsenOperation getRefineOperation() 
    { 
      return refineOperation; 
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
423
    /// Restriction after coarsening. Implemented for DOFVector<double>
424
    inline void coarseRestrict(RCNeighbourList&, int) {}
425

426
    /// Interpolation after refinement. Implemented for DOFVector<double>
427
    inline void refineInterpol(RCNeighbourList&, int) {}
428

Thomas Witkowski's avatar
Thomas Witkowski committed
429
    /// Returns \ref vec
430
    std::vector<T>& getVector() 
Thomas Witkowski's avatar
Thomas Witkowski committed
431
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
432
      return vec;
433
    }
434

Thomas Witkowski's avatar
Thomas Witkowski committed
435
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
    inline int getSize() const 
    { 
438
      return vec.size();
439
    } 
440

Thomas Witkowski's avatar
Thomas Witkowski committed
441
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
442
443
    inline int getUsedSize() const 
    { 
444
      return this->feSpace->getAdmin()->getUsedSize(); 
445
    }
446

Thomas Witkowski's avatar
Thomas Witkowski committed
447
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
448
449
    inline void resize(int n) 
    { 
450
      FUNCNAME_DBG("DOFVector<T>::resize()");
451
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
452
      vec.resize(n);
453
    } 
454

Thomas Witkowski's avatar
Thomas Witkowski committed
455
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
    inline void resize(int n, T init) 
    { 
458
      FUNCNAME_DBG("DOFVector<T>::resize()");
459
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
460
      vec.resize(n, init);
461
    } 
462

Thomas Witkowski's avatar
Thomas Witkowski committed
463
    /// Returns \ref vec[i]
464
    inline const_reference operator[](DegreeOfFreedom i) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
465
    {
466
      FUNCNAME_DBG("DOFVector<T>::operator[]");
467
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
468
	("Illegal vector index %d.\n", i);
469
      return vec[i];
470
    } 
471

Thomas Witkowski's avatar
Thomas Witkowski committed
472
    /// Returns \ref vec[i]
473
    inline reference operator[](DegreeOfFreedom i) 
Thomas Witkowski's avatar
Thomas Witkowski committed
474
    {
475
      FUNCNAME_DBG("DOFVector<T>::operator[]");
476
477

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

480
      return vec[i];
481
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
482
 
Thomas Witkowski's avatar
Thomas Witkowski committed
483
    /// Calculates Integral of this DOFVector
484
    double Int(Quadrature* q = NULL) const
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
487
488
489
490
491
492
493
494
495
496
497
498
    {
      return Int(-1, q);
    }

    /** \brief
     * Calculates Integral of this DOFVector. 
     *
     * \param[in]  meshlevel  Then mesh level on which the integral should be
     *                        calculated. Usually, only -1 for the leaf level
     *                        makes sence, but other values can be used, e.g.,
     *                        for debugging.
     * \param[in]  q          Quadrature object. If not specified, the function
     *                        creates a new quadrature object.
     */
499
    double Int(int meshLevel, Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
500

501

Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
504
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
505
506
507
508
509
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
510

Thomas Witkowski's avatar
Thomas Witkowski committed
511
512
513
    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
514
    double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = NULL) const
515
516
517
518
519
520
    {
      FUNCNAME("DOFVector::IntOnBoundaryNormal())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
521
    /// Calculates L1 norm of this DOFVector
522
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
523
 
Thomas Witkowski's avatar
Thomas Witkowski committed
524
    /// Calculates L2 norm of this DOFVector
525
    inline double L2Norm(Quadrature* q = NULL) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
526
    {
Wieland Marth's avatar
Wieland Marth committed
527
      return std::sqrt(L2NormSquare());
528
    }
529

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

Thomas Witkowski's avatar
Thomas Witkowski committed
533
    /// Calculates H1 norm of this DOFVector
534
    inline double H1Norm(Quadrature* q = NULL) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
535
    {
Wieland Marth's avatar
Wieland Marth committed
536
      return std::sqrt(H1NormSquare());
537
    }
538

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

Thomas Witkowski's avatar
Thomas Witkowski committed
542
    /// Calculates euclidian norm of this DOFVector
543
544
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
545
    /// Returns square of the euclidian norm.
546
547
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
548
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
549
550
    inline double l2norm() const 
    { 
551
      return nrm2();
552
    }
553

Thomas Witkowski's avatar
Thomas Witkowski committed
554
    /// Calculates the absolute sum of this DOFVector
555
556
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
557
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
558
559
    inline double l1norm() const 
    { 
560
      return asum();
561
    } 
562

Thomas Witkowski's avatar
Thomas Witkowski committed
563
    /// Calculates doublewell of this DOFVector
564
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
565
 
Thomas Witkowski's avatar
Thomas Witkowski committed
566
    /// Calculates the sum of this DOFVector
567
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
568
 
Thomas Witkowski's avatar
Thomas Witkowski committed
569
    /// Sets \ref vec[i] = val, i=0 , ... , size
570
571
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
572
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
573
574
    inline DOFVector<T>& operator=(T d) 
    {
575
576
      set(d); 
      return *this;
577
    } 
578

Thomas Witkowski's avatar
Thomas Witkowski committed
579
    /// Assignment operator between two vectors
580
    DOFVector<T>& operator=(const DOFVector<T>&); 
581

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

585
    /// Returns minimum of DOFVector.
586
587
    T min() const; 

588
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
589
590
    T max() const;

591
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
592
    T absMax() const;
593

594
595
596
597
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
598
599
    void print() const; 

600
    ///
601
    size_t calcMemoryUsage() const;
602

Thomas Witkowski's avatar
Thomas Witkowski committed
603
604
    /// Computes the coefficients of the interpolant of the function fct and
    /// stores these in the DOFVector
605
606
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

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

609

610
    /// Eval DOFVector at given point p. If oldElInfo != NULL the search for 
Thomas Witkowski's avatar
Thomas Witkowski committed
611
612
    /// the element, where p is inside, starts from oldElInfo. implemented for:
    /// double, WorldVector< double >
613
    T evalAtPoint(WorldVector<double> &p, 
614
		  ElInfo *oldElInfo = NULL) const 
615
616
617
618
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
    }
619

Thomas Witkowski's avatar
Thomas Witkowski committed
620
621
622
    /// Determine the DegreeOfFreedom that has coords with minimal euclidean 
    /// distance to WorldVector p. return true if DOF is found, and false 
    /// otherwise.
623
    bool getDofIdxAtPoint(WorldVector<double> &p, 
Thomas Witkowski's avatar
Thomas Witkowski committed
624
				DegreeOfFreedom &idx, 
625
				ElInfo *oldElInfo = NULL, 
Thomas Witkowski's avatar
Thomas Witkowski committed
626
				bool useOldElInfo = false) const;
627

628
629
    template<typename S>
    void serialize(std::ostream &out, boost::mpl::identity<S>) 
Thomas Witkowski's avatar
Thomas Witkowski committed
630
    {
631
      unsigned int size = static_cast<unsigned int>(vec.size());
632
633
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
634
    }
635
636
637
638
639
640
641
642
643
644
645
646
    
    // no serialization for bool vectors
    void serialize(std::ostream &out, boost::mpl::identity<bool>) {}
    
    /// Writes the data of the DOFVector to an output stream.
    void serialize(std::ostream &out) 
    {
      serialize(out, boost::mpl::identity<T>());
    }
    
    template<typename S>
    void deserialize(std::istream &in, boost::mpl::identity<S>) 
Thomas Witkowski's avatar
Thomas Witkowski committed
647
    {
648
649
650
651
      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));
652
    }
653
654
655
656
657
658
659
660
661
    
    // no de-serialization for bool vectors
    void deserialize(std::istream &in, boost::mpl::identity<bool>) {}

    /// Reads data of an DOFVector from an input stream.
    void deserialize(std::istream &in) 
    {
      deserialize(in, boost::mpl::identity<T>());
    }
662

663
664
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
665
666
667

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

668
669
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
670
671

  protected: 
672

673
    /// Data container
674
    std::vector<T> vec; 
675
 
676
    /// Specifies what operation should be performed after coarsening
677
678
    RefineCoarsenOperation coarsenOperation;
    RefineCoarsenOperation refineOperation;
679
680
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
681

682
  template<>
683
684
685
686
687
688
689
690
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

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

  template<>
691
692
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
					 ElInfo *oldElInfo) const;
693
694

  template<>
695
696
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const;
697

698
  template<>
699
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
700
701

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

Thomas Witkowski's avatar
Thomas Witkowski committed
704
705
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
706
    return v.min();
707
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
708

Thomas Witkowski's avatar
Thomas Witkowski committed
709
710
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
711
    return v.max();
712
  }
713
714
715
716
717
718
719
720
721
722


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Thomas Witkowski's avatar
Thomas Witkowski committed
723
724
    /// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
    /// as DOFContainer at DOFAdmin
725
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
726
727
728
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
729
    }
730
  
731
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
732
733
    ~DOFVectorDOF() 
    {
734
      feSpace->getAdmin()->removeDOFContainer(this);
735
    }
736

Thomas Witkowski's avatar
Thomas Witkowski committed
737
738
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
741
      return DOFVector<DegreeOfFreedom>::operator[](i);
742
    }
743

Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
746
747
748
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

749
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
    int getSize() const 
    {
752
      return DOFVector<DegreeOfFreedom>::getSize();
753
    }
754

755
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
756
757
    void resize(int size) 
    {
758
759
760
761
762
763
764
765
766
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
767

768
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
769
770
  double norm(DOFVector<T> *vec) 
  {
771
    return vec->nrm2();
772
  }
773
774

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
  double L2Norm(DOFVector<T> *vec) 
  {
777
    return vec->L2Norm();