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
Blub    
Thomas Witkowski committed
240
241
242
243
244
245
246
247
248
249
250
251
    inline void setDirichletDofValue(DegreeOfFreedom dof,
				     T value)
    {
      dirichletDofValues[dof] = value;
    }

    map<DegreeOfFreedom, T>& getDirichletValues()
    {
      return dirichletDofValues;
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
252
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
253
    void setDofMap(FeSpaceDofMap& m)
Thomas Witkowski's avatar
Thomas Witkowski committed
254
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
255
      dofMap = &m;
Thomas Witkowski's avatar
Thomas Witkowski committed
256
257
258
259
    }

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

262
      return dofMap->isRankDof(dof);
Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
265
    }
#endif

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
294
295
    map<DegreeOfFreedom, T> dirichletDofValues;

Thomas Witkowski's avatar
Thomas Witkowski committed
296
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
297
    FeSpaceDofMap *dofMap;
298
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
299
  };
300

301

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

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

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

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

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

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

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

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

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

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

    /// 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
410
411
    /// Used by DOFAdmin to compress this DOFVector. Implementation of
    /// \ref DOFIndexedBase::compress()
412
    virtual void compressDOFIndexed(int first, int last,
Thomas Witkowski's avatar
Thomas Witkowski committed
413
				    vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
414
415

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

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

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

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

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

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

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

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

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

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

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

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

484
      return vec[i];
485
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
486
 
Thomas Witkowski's avatar
Thomas Witkowski committed
487
    /// Calculates Integral of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
    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;

505

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

613

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

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

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

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

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

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

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

  protected: 
658

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

Thomas Witkowski's avatar
Thomas Witkowski committed
666

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

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

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

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

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

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

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

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


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

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
752

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

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

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

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

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

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

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

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


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

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

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

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

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

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

  
  /// 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
912
		  AbstractFunction<TOut, T> *fct = NULL)
913
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
914
    return fct == NULL ? vec.Int() : integrate_Vec(vec, fct);
915
  }
916
  
917
  /** \brief
918
919
   * Computes the integral of a function that includes two different DOFVectors:
   * \f$ \int f(v(\vec{x}), w(\vec{x}))\,\text{d}\vec{x}\f$