DOFVector.h 26.1 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
30
31
#include <vector> 
#include <memory> 
#include <list> 

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

48
49
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
50
51
  using namespace std;

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

57
58
    DOFVectorBase() 
      : feSpace(NULL),
59
	elementVector(3),
Thomas Witkowski's avatar
Thomas Witkowski committed
60
61
        boundaryManager(NULL),
        nBasFcts(0)
62
    {}
Thomas Witkowski's avatar
Thomas Witkowski committed
63
    
Thomas Witkowski's avatar
Thomas Witkowski committed
64
    DOFVectorBase(const FiniteElemSpace *f, string n);
65

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
86
87
    /// Evaluates the gradient of a DOF vector at a set of quadrature points
    /// defined on the given element.
88
89
90
91
92
93
94
95
96
97
98
    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
99
100
    /// Evaluates the comp'th component of the derivative of a DOF vector at a
    /// set of quadrature points defined on the given element.
101
102
103
104
105
106
107
108
109
110
111
112
113
    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
114
115
    /// Evaluates the jacobian of a DOF vector at a set of quadrature points
    ///  defined on the given element.
116
117
118
119
    void getD2AtQPs(const ElInfo *elInfo,
		    const Quadrature *quad,
		    const FastQuadrature *quadFast,
		    mtl::dense_vector<typename D2Type<T>::type> &d2AtQPs) const;
120

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

Thomas Witkowski's avatar
Thomas Witkowski committed
127
128
    /// 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
129
130
131
132
133
134
135
136
137
    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);
138
139

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

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

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

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

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

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


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

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

    Flag getAssembleFlag();

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
219
    inline void setName(string n)
220
    {
221
      name = n;
222
223
    }

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
249
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
250
    string name;
251

Thomas Witkowski's avatar
Thomas Witkowski committed
252
    ///
253
    ElementVector elementVector;
254

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
264
    ///
265
266
    BoundaryManager *boundaryManager;

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

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
273
    map<DegreeOfFreedom, T> dirichletDofValues;
Thomas Witkowski's avatar
Thomas Witkowski committed
274
  };
275

276

Thomas Witkowski's avatar
Thomas Witkowski committed
277
  /// Specifies which operation should be done after coarsening
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
  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)
309
      {}
310
311
312

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
313
      {}
314
315
316
317
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
318
319
320
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
321

Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
      DOFVector<T> *create() 
      { 
324
	return new DOFVector<T>(feSpace, ""); 
325
      }
326

Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
      void free(DOFVector<T> *vec) 
      { 
329
	delete vec; 
330
      }
331
332
333
334
335
336

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
337
    /// Empty constructor. No initialization!
338
339
340
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
341
    {}
342

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

Thomas Witkowski's avatar
Thomas Witkowski committed
346
    /// Initialization.
Thomas Witkowski's avatar
Thomas Witkowski committed
347
    void init(const FiniteElemSpace* f, string n);
348

Thomas Witkowski's avatar
Thomas Witkowski committed
349
    /// Copy Constructor
350
    DOFVector(const DOFVector& rhs) : DOFVectorBase<T>()
Thomas Witkowski's avatar
Thomas Witkowski committed
351
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
352
      *this = rhs;   
353
354
355
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
356
    }
357

Thomas Witkowski's avatar
Thomas Witkowski committed
358
    /// Destructor
359
360
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    /// Returns iterator to the begin of \ref vec
362
    typename std::vector<T>::iterator begin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
363
    { 
364
      return vec.begin(); 
365
    }
366

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

    /// Returns const_iterator to the begin of \ref vec
374
    typename std::vector<T>::const_iterator begin() const
Praetorius, Simon's avatar
Praetorius, Simon committed
375
376
377
378
379
    {
      return vec.begin();
    }

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

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

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

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

405
    /// Interpolation after refinement. Implemented for DOFVector<double>
406
    inline void refineInterpol(RCNeighbourList&, int) {}
407

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

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

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

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

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

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

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

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

459
      return vec[i];
460
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
461
 
Thomas Witkowski's avatar
Thomas Witkowski committed
462
    /// Calculates Integral of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
    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;

480

Thomas Witkowski's avatar
Thomas Witkowski committed
481
482
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
483
484
485
486
487
488
    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
489

Thomas Witkowski's avatar
Thomas Witkowski committed
490
491
492
    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
493
494
495
496
497
498
499
    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
500
    /// Calculates L1 norm of this DOFVector
501
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
502
 
Thomas Witkowski's avatar
Thomas Witkowski committed
503
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
    inline double L2Norm(Quadrature* q = NULL) const 
    {
506
      return sqrt(L2NormSquare());
507
    }
508

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

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

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

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

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

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() const 
    { 
530
      return nrm2();
531
    }
532

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

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

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

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

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

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

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

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

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

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

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

579
    ///
580
581
    int calcMemoryUsage() const;

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

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

588

Thomas Witkowski's avatar
Thomas Witkowski committed
589
590
591
    /// 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 >
592
593
    T evalAtPoint(WorldVector<double> &p, 
		  ElInfo *oldElInfo = NULL) const 
594
595
596
597
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
    }
598

Thomas Witkowski's avatar
Thomas Witkowski committed
599
600
601
    /// Determine the DegreeOfFreedom that has coords with minimal euclidean 
    /// distance to WorldVector p. return true if DOF is found, and false 
    /// otherwise.
602
    bool getDofIdxAtPoint(WorldVector<double> &p, 
Thomas Witkowski's avatar
Thomas Witkowski committed
603
604
605
				DegreeOfFreedom &idx, 
				ElInfo *oldElInfo = NULL, 
				bool useOldElInfo = false) const;
606

607
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
608
    void serialize(ostream &out) 
Thomas Witkowski's avatar
Thomas Witkowski committed
609
    {
610
      unsigned int size = static_cast<unsigned int>(vec.size());
611
612
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
613
    }
614

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

624
625
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
626
627
628

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

629
630
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
631
632

  protected: 
633

634
    /// Data container
635
    std::vector<T> vec; 
636
 
637
    /// Specifies what operation should be performed after coarsening
638
639
640
    CoarsenOperation coarsenOperation;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
641

642
  template<>
643
644
645
646
647
648
649
650
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

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

  template<>
651
652
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
					 ElInfo *oldElInfo) const;
653
654

  template<>
655
656
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const;
657

658
  template<>
659
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
660
661

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

Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
666
    return v.min();
667
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
668

Thomas Witkowski's avatar
Thomas Witkowski committed
669
670
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
671
    return v.max();
672
  }
673
674
675
676
677
678
679
680
681
682


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Thomas Witkowski's avatar
Thomas Witkowski committed
683
684
685
    /// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
    /// as DOFContainer at DOFAdmin
    DOFVectorDOF(const FiniteElemSpace* feSpace_, string name_)
686
687
688
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
689
    }
690
  
691
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
    ~DOFVectorDOF() 
    {
694
      feSpace->getAdmin()->removeDOFContainer(this);
695
    }
696

Thomas Witkowski's avatar
Thomas Witkowski committed
697
698
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
701
      return DOFVector<DegreeOfFreedom>::operator[](i);
702
    }
703

Thomas Witkowski's avatar
Thomas Witkowski committed
704
705
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
706
707
708
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

709
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
710
711
    int getSize() const 
    {
712
      return DOFVector<DegreeOfFreedom>::getSize();
713
    }
714

715
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
    void resize(int size) 
    {
718
719
720
721
722
723
724
725
726
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
727

728
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
729
730
  double norm(DOFVector<T> *vec) 
  {
731
    return vec->nrm2();
732
  }
733
734

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
  double L2Norm(DOFVector<T> *vec) 
  {
737
    return vec->L2Norm();
738
  }
739
740

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
741
742
  double H1Norm(DOFVector<T> *vec) 
  {
743
    return vec->H1Norm();
744
  }
745
746

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
747
  void print(DOFVector<T> *vec) 
748
  {
749
    vec->print();
750
  }
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782

  // 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>
783
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
784
785
786
787
788
789
790
791
792
793
794
795
796

  // 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
797
798
  inline void scal(T a, DOFVector<T>& y) 
  {
799
800
    y *= a;
  }
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816

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

817
818
819
820
821
822
823
824
825
  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);


826
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
827
828
  inline void set(DOFVector<T>& vec, T d) 
  {
829
    vec.set(d);
830
  }
831
832

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
833
834
  inline void setValue(DOFVector<T>& vec, T d) 
  {
835
    vec.set(d);
836
  }
837
838

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
839
840
  inline int size(DOFVector<T> *vec) 
  {
841
    return vec->getUsedSize();
842
  }
843

844
845
846
847
848
849
  template<typename T>
  inline int size(const DOFVector<T>& vec) 
  {
    return vec.getUsedSize();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
850
  template<typename T>
851
  inline void checkFeSpace(const FiniteElemSpace* feSpace, const std::vector<T>& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
852
  {
853
    FUNCNAME_DBG("checkFeSpace()");
Thomas Witkowski's avatar
Thomas Witkowski committed
854
855
856
857
858
859
860
    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());
  }

861
862
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
863
864
  
  template<typename T>
865
866
  std::vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     std::vector<DOFVector<double>*> *res);
867
868

  
869
  /// Computes the integral: \f$ \int f(\vec{x})\,d\vec{x}\f$
870
871
872
873
874
875
876
877
878
879
880
  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);
  }
  
881
  /// Computes the integral: \f$ \int f(v(\vec{x}))\,d\vec{x}\f$
882
883
884
885
886
887
  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
888
		  AbstractFunction<TOut, T> *fct = NULL)
889
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
890
    return fct == NULL ? vec.Int() : integrate_Vec(vec, fct);
891
  }
892
  
893
  /** \brief
894
   * Computes the integral of a function that includes two different DOFVectors:
895
   * \f$ \int f(v(\vec{x}), w(\vec{x}))\,d\vec{x}\f$
Thomas Witkowski's avatar
Thomas Witkowski committed
896
897
   * This function works also for the case that the DOFVectors are defined on 
   * two different meshes.
898
   */
899
900
901
902
903
904
905
906
907
908
909
910
911
  template<typename TOut, typename T1, typename T2>
  TOut integrate_Vec2(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct);

  template<typename TOut, typename T1, typename T2>
  TOut integrate(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct)
  {
    return integrate_Vec2(vec1, vec2, fct);
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
912
913
  /// Computes the integral of a function with two DOFVectors defined on the 
  /// same mesh.
914
915
916
917
  template<typename TOut, typename T1, typename T2>
  TOut int_Vec2_SingleMesh(const DOFVector<T1> &vec1,
		       const DOFVector<T2> &vec2,
		       BinaryAbstractFunction<TOut, T1, T2> *fct);