DOFVector.h 26.6 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"
Thomas Witkowski's avatar
Thomas Witkowski committed
45
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
46
#include <mpi.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
47
48
#include "parallel/ParallelDofMapping.h"
#endif
49

50
51
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
  using namespace std;

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

59
60
    DOFVectorBase() 
      : feSpace(NULL),
61
	elementVector(3),
Thomas Witkowski's avatar
Thomas Witkowski committed
62
63
        boundaryManager(NULL),
        nBasFcts(0)
64
65
    {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
66
      dofMap = NULL;
67
68
#endif
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
69
    
Thomas Witkowski's avatar
Thomas Witkowski committed
70
    DOFVectorBase(const FiniteElemSpace *f, string n);
71

Thomas Witkowski's avatar
Thomas Witkowski committed
72
    virtual ~DOFVectorBase();
73

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

Thomas Witkowski's avatar
Thomas Witkowski committed
79
80
    /// 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

Thomas Witkowski's avatar
Thomas Witkowski committed
92
93
    /// 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;

Thomas Witkowski's avatar
Thomas Witkowski committed
105
106
    /// Evaluates the comp'th component of the derivative of a DOF vector at a
    /// set of quadrature points defined on the given element.
107
108
109
110
111
112
113
114
115
116
117
118
119
    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
120
121
    /// Evaluates the jacobian of a DOF vector at a set of quadrature points
    ///  defined on the given element.
122
123
124
125
    void getD2AtQPs(const ElInfo *elInfo,
		    const Quadrature *quad,
		    const FastQuadrature *quadFast,
		    mtl::dense_vector<typename D2Type<T>::type> &d2AtQPs) const;
126

Thomas Witkowski's avatar
Thomas Witkowski committed
127
    /// Returns the FE space the DOF vector is defined on.
128
    inline const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
129
    {
130
      return feSpace;
131
    }
132

Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
    /// 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
135
136
137
138
139
140
141
142
143
    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);
144
145

    void addElementVector(T sign,
146
			  const ElementVector& elVec, 
147
			  const BoundaryType *bound,
148
			  ElInfo *elInfo,
149
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
150

Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153
    /// 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
154
    void finishAssembling();
155
156
157
158
159
160
161
162
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
163
    }
164

Thomas Witkowski's avatar
Thomas Witkowski committed
165
    inline vector<double*>::iterator getOperatorFactorBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
166
    {
167
      return operatorFactor.begin();
168
    }
169

Thomas Witkowski's avatar
Thomas Witkowski committed
170
    inline vector<double*>::iterator getOperatorFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
171
    {
172
      return operatorFactor.end();
173
    }
174

Thomas Witkowski's avatar
Thomas Witkowski committed
175
    inline vector<double*>::iterator getOperatorEstFactorBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
176
    {
177
      return operatorEstFactor.begin();
178
    }
179

Thomas Witkowski's avatar
Thomas Witkowski committed
180
    inline vector<double*>::iterator getOperatorEstFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
181
    {
182
      return operatorEstFactor.end();
183
    }
184
185


Thomas Witkowski's avatar
Thomas Witkowski committed
186
    inline vector<Operator*>::iterator getOperatorsBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
187
    {
188
      return operators.begin();
189
    }
190

Thomas Witkowski's avatar
Thomas Witkowski committed
191
    inline vector<Operator*>::iterator getOperatorsEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
192
    {
193
      return operators.end();
194
    }
195
196
197

    Flag getAssembleFlag();

198
199
200
201
    /// 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
202
203
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

Thomas Witkowski's avatar
Thomas Witkowski committed
204
    inline vector<Operator*>& getOperators() 
Thomas Witkowski's avatar
Thomas Witkowski committed
205
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
206
      return operators; 
207
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
208

Thomas Witkowski's avatar
Thomas Witkowski committed
209
    inline vector<double*>& getOperatorFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
210
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
211
      return operatorFactor; 
212
    }
213

Thomas Witkowski's avatar
Thomas Witkowski committed
214
    inline vector<double*>& getOperatorEstFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
215
    { 
216
      return operatorEstFactor; 
217
    }
218

Thomas Witkowski's avatar
Thomas Witkowski committed
219
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
220
    inline string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
221
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
222
      return name; 
223
    } 
224

Thomas Witkowski's avatar
Thomas Witkowski committed
225
    inline void setName(string n)
226
    {
227
      name = n;
228
229
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
230
231
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
232
      return boundaryManager; 
233
    }
234

Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
237
      boundaryManager = bm;
238
    }
239

Thomas Witkowski's avatar
Thomas Witkowski committed
240
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
241
    void setDofMap(FeSpaceDofMap& m)
Thomas Witkowski's avatar
Thomas Witkowski committed
242
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
243
      dofMap = &m;
Thomas Witkowski's avatar
Thomas Witkowski committed
244
245
246
247
    }

    inline bool isRankDof(DegreeOfFreedom dof)
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
248
      TEST_EXIT_DBG(dofMap)("No rank dofs set!\n");
249

250
      return dofMap->isRankDof(dof);
Thomas Witkowski's avatar
Thomas Witkowski committed
251
    }
252
253
254
255
256
257
258
259
260
261
262

    inline void setDirichletDofValue(DegreeOfFreedom dof,
				     T value)
    {
      dirichletDofValues[dof] = value;
    }

    map<DegreeOfFreedom, T>& getDirichletValues()
    {
      return dirichletDofValues;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
#endif

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

Thomas Witkowski's avatar
Thomas Witkowski committed
269
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
270
    string name;
271

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

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

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

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

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

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

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

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
294
    FeSpaceDofMap *dofMap;
295

296
297
    map<DegreeOfFreedom, T> dirichletDofValues;
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
298
  };
299

300

Thomas Witkowski's avatar
Thomas Witkowski committed
301
  /// Specifies which operation should be done after coarsening
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
330
331
332
  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)
333
      {}
334
335
336

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
337
      {}
338
339
340
341
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
342
343
344
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
345

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

Thomas Witkowski's avatar
Thomas Witkowski committed
367
    /// Constructs a DOFVector with name n belonging to FiniteElemSpace f
Thomas Witkowski's avatar
Thomas Witkowski committed
368
    DOFVector(const FiniteElemSpace* f, string n); 
369

Thomas Witkowski's avatar
Thomas Witkowski committed
370
    /// Initialization.
Thomas Witkowski's avatar
Thomas Witkowski committed
371
    void init(const FiniteElemSpace* f, string n);
372

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

Thomas Witkowski's avatar
Thomas Witkowski committed
382
    /// Destructor
383
384
    virtual ~DOFVector();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
391
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
392
    typename vector<T>::iterator end() 
Thomas Witkowski's avatar
Thomas Witkowski committed
393
    { 
394
      return vec.end(); 
395
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
396
397
398
399
400
401
402
403
404
405
406
407
408

    /// Returns const_iterator to the begin of \ref vec
    typename vector<T>::const_iterator begin() const
    {
      return vec.begin();
    }

    /// Returns const_iterator to the end of \ref vec
    typename vector<T>::const_iterator end() const
    {
      return vec.end();
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
409
410
    /// Used by DOFAdmin to compress this DOFVector. Implementation of
    /// \ref DOFIndexedBase::compress()
411
    virtual void compressDOFIndexed(int first, int last,
Thomas Witkowski's avatar
Thomas Witkowski committed
412
				    vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
413
414

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
417
      coarsenOperation = op; 
418
    }
419

Thomas Witkowski's avatar
Thomas Witkowski committed
420
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
    inline CoarsenOperation getCoarsenOperation() 
    { 
423
      return coarsenOperation; 
424
    }
425

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

429
    /// Interpolation after refinement. Implemented for DOFVector<double>
430
    inline void refineInterpol(RCNeighbourList&, int) {}
431

Thomas Witkowski's avatar
Thomas Witkowski committed
432
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
433
    vector<T>& getVector() 
Thomas Witkowski's avatar
Thomas Witkowski committed
434
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
435
      return vec;
436
    }
437

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
475
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
    inline T& operator[](DegreeOfFreedom i) 
    {
478
      FUNCNAME("DOFVector<T>::operator[]");
479
480

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

483
      return vec[i];
484
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
485
 
Thomas Witkowski's avatar
Thomas Witkowski committed
486
    /// Calculates Integral of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
    double Int(Quadrature* q = NULL) const
    {
      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.
     */
    double Int(int meshLevel, Quadrature* q = NULL) const;

504

Thomas Witkowski's avatar
Thomas Witkowski committed
505
506
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
507
508
509
510
511
512
    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
513

Thomas Witkowski's avatar
Thomas Witkowski committed
514
515
516
    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
517
518
519
520
521
522
523
    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
524
    /// Calculates L1 norm of this DOFVector
525
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
526
 
Thomas Witkowski's avatar
Thomas Witkowski committed
527
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
528
529
    inline double L2Norm(Quadrature* q = NULL) const 
    {
530
      return sqrt(L2NormSquare());
531
    }
532

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

Thomas Witkowski's avatar
Thomas Witkowski committed
536
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
537
538
    inline double H1Norm(Quadrature* q = NULL) const 
    {
539
      return sqrt(H1NormSquare());
540
    }
541

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

Thomas Witkowski's avatar
Thomas Witkowski committed
545
    /// Calculates euclidian norm of this DOFVector
546
547
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
548
    /// Returns square of the euclidian norm.
549
550
    double squareNrm2() const;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
557
    /// Calculates the absolute sum of this DOFVector
558
559
    T asum() const; 

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
582
    /// Assignment operator between two vectors
583
    DOFVector<T>& operator=(const DOFVector<T>&); 
584

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

588
    /// Returns minimum of DOFVector.
589
590
    T min() const; 

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

594
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
595
    T absMax() const;
596

597
598
599
600
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
601
602
    void print() const; 

603
    ///
604
605
    int calcMemoryUsage() const;

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

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

612

Thomas Witkowski's avatar
Thomas Witkowski committed
613
614
615
    /// 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 >
616
617
    T evalAtPoint(WorldVector<double> &p, 
		  ElInfo *oldElInfo = NULL) const 
618
619
620
621
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
    }
622

Thomas Witkowski's avatar
Thomas Witkowski committed
623
624
625
626
627
628
629
    /// 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;
630

631
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
632
    void serialize(ostream &out) 
Thomas Witkowski's avatar
Thomas Witkowski committed
633
    {
634
635
636
      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));
637
    }
638

639
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
640
    void deserialize(istream &in) 
Thomas Witkowski's avatar
Thomas Witkowski committed
641
    {
642
643
644
645
      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));
646
    }
647

648
649
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
650
651
652

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

653
654
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
655
656

  protected: 
657

658
    /// Data container
Thomas Witkowski's avatar
Thomas Witkowski committed
659
    vector<T> vec; 
660
 
661
    /// Specifies what operation should be performed after coarsening
662
663
664
    CoarsenOperation coarsenOperation;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
665

666
  template<>
667
668
669
670
671
672
673
674
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

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

  template<>
675
676
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
					 ElInfo *oldElInfo) const;
677
678

  template<>
679
680
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const;
681

682
  template<>
683
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
684
685

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

Thomas Witkowski's avatar
Thomas Witkowski committed
688
689
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
690
    return v.min();
691
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
692

Thomas Witkowski's avatar
Thomas Witkowski committed
693
694
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
695
    return v.max();
696
  }
697
698
699
700
701
702
703
704
705
706


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Thomas Witkowski's avatar
Thomas Witkowski committed
707
708
709
    /// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
    /// as DOFContainer at DOFAdmin
    DOFVectorDOF(const FiniteElemSpace* feSpace_, string name_)
710
711
712
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
713
    }
714
  
715
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
    ~DOFVectorDOF() 
    {
718
      feSpace->getAdmin()->removeDOFContainer(this);
719
    }
720

Thomas Witkowski's avatar
Thomas Witkowski committed
721
722
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
723
724
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
725
      return DOFVector<DegreeOfFreedom>::operator[](i);
726
    }
727

Thomas Witkowski's avatar
Thomas Witkowski committed
728
729
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
730
731
732
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

733
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
734
735
    int getSize() const 
    {
736
      return DOFVector<DegreeOfFreedom>::getSize();
737
    }
738

739
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
740
741
    void resize(int size) 
    {
742
743
744
745
746
747
748
749
750
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
751

752
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
753
754
  double norm(DOFVector<T> *vec) 
  {
755
    return vec->nrm2();
756
  }
757
758

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
  double L2Norm(DOFVector<T> *vec) 
  {
761
    return vec->L2Norm();
762
  }
763
764

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
765
766
  double H1Norm(DOFVector<T> *vec) 
  {
767
    return vec->H1Norm();
768
  }
769
770

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
771
  void print(DOFVector<T> *vec) 
772
  {
773
    vec->print();
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

  // 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>
807
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
808
809
810
811
812
813
814
815
816
817
818
819
820

  // 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
821
822
  inline void scal(T a, DOFVector<T>& y) 
  {
823
824
    y *= a;
  }
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840

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

841
842
843
844
845
846
847
848
849
  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);


850
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
851
852
  inline void set(DOFVector<T>& vec, T d) 
  {
853
    vec.set(d);
854
  }
855
856

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
857
858
  inline void setValue(DOFVector<T>& vec, T d) 
  {
859
    vec.set(d);
860
  }
861
862

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
863
864
  inline int size(DOFVector<T> *vec) 
  {
865
    return vec->getUsedSize();
866
  }
867

868
869
870
871
872
873
  template<typename T>
  inline int size(const DOFVector<T>& vec) 
  {
    return vec.getUsedSize();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
874
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
875
  inline void checkFeSpace(const FiniteElemSpace* feSpace, const vector<T>& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
876
877
878
879
880
881
882
883
  {
    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());
  }

884
885
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
886
887
  
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
888
889
  vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     vector<DOFVector<double>*> *res);
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910

  
  /// Computes the integral: \f$ \int f(\vec{x})\,\text{d}\vec{x}\f$
  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct);
  
  template<typename TOut>
  TOut integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
  {
    return integrate_Coords(feSpace, fct);
  }
  
  /// Computes the integral: \f$ \int f(v(\vec{x}))\,\text{d}\vec{x}\f$
  template<typename TOut, typename T>
  TOut integrate_Vec(const DOFVector<T> &vec,
		  AbstractFunction<TOut, T> *fct);
  
  template<typename TOut, typename T>
  TOut integrate(const DOFVector<T> &vec,
Praetorius, Simon's avatar
Praetorius, Simon committed
911
		  AbstractFunction<TOut, T> *fct = NULL)
912
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
913
    return fct == NULL ? vec.Int() : integrate_Vec(vec, fct);
914
  }
915