DOFVector.h 26.3 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20
21
22
23
24
25

/** \file DOFVector.h */

#ifndef AMDIS_DOFVECTOR_H 
#define AMDIS_DOFVECTOR_H 
 
26
27
28
29
#include <vector> 
#include <memory> 
#include <list> 

30
#include "AMDiS_fwd.h"
31
32
33
34
35
36
37
38
39
40
41
#include "FixVec.h"
#include "Global.h" 
#include "Flag.h" 
#include "RCNeighbourList.h" 
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "DOFContainer.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
Thomas Witkowski's avatar
Thomas Witkowski committed
42
#include "BasisFunction.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
43
#include "FiniteElemSpace.h"
44
#include "SurfaceQuadrature.h"
45

46
47
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
48
49
  using namespace std;

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
64
    virtual ~DOFVectorBase();
65

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

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

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

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

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

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

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

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

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

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

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

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


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

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

    Flag getAssembleFlag();

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

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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
247
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
248
    string name;
249

Thomas Witkowski's avatar
Thomas Witkowski committed
250
    ///
251
    ElementVector elementVector;
252

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    ///
263
264
    BoundaryManager *boundaryManager;

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

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

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

274

Thomas Witkowski's avatar
Thomas Witkowski committed
275
  /// Specifies which operation should be done after coarsening
276
277
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
  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)
307
      {}
308
309
310

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

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

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
356
    /// Destructor
357
358
    virtual ~DOFVector();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

478

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    /// Prints \ref vec to stdout.
575
576
    void print() const; 

577
    ///
578
579
    int calcMemoryUsage() const;

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

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

586

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

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

605
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
606
    void serialize(ostream &out) 
Thomas Witkowski's avatar
Thomas Witkowski committed
607
    {
608
609
610
      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));
611
    }
612

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

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

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

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

  protected: 
631

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

Thomas Witkowski's avatar
Thomas Witkowski committed
639

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

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

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

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

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

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

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

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


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

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
725

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

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
745
  void print(DOFVector<T> *vec) 
746
  {
747
    vec->print();
748
  }
749
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

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

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

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

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


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

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

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

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

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

859
860
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
861
862
  
  template<typename T>
863
864
  std::vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     std::vector<DOFVector<double>*> *res);
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885

  
  /// 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
886
		  AbstractFunction<TOut, T> *fct = NULL)
887
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
888
    return fct == NULL ? vec.Int() : integrate_Vec(vec, fct);
889
  }
890
  
891
  /** \brief
892
893
   * 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$
Thomas Witkowski's avatar
Thomas Witkowski committed
894
895
   * This function works also for the case that the DOFVectors are defined on 
   * two different meshes.
896
   */
897
898
899
900
901
902
903
904
905
906
907
908
909
  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
910
911
  /// Computes the integral of a function with two DOFVectors defined on the 
  /// same mesh.
912
913
914
915
  template<typename TOut, typename T1, typename T2>
  TOut int_Vec2_SingleMesh(const DOFVector<T1> &vec1,
		       const DOFVector<T2> &vec2,
		       BinaryAbstractFunction<TOut, T1, T2> *fct);