DOFVector.h 23.9 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
48
49
50
51
namespace AMDiS {

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

53
54
    DOFVectorBase() 
      : feSpace(NULL),
55
	elementVector(3),
Thomas Witkowski's avatar
Thomas Witkowski committed
56
57
        boundaryManager(NULL),
        nBasFcts(0)
58
59
60
61
62
    {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      rankDofs = NULL;
#endif
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
63
    
64
    DOFVectorBase(const FiniteElemSpace *f, std::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

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

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

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

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


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

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

    Flag getAssembleFlag();

    /** \brief
     * 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
     */
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
236
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
237
    inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
238
    {
239
      rankDofs = &dofmap;
Thomas Witkowski's avatar
Thomas Witkowski committed
240
241
242
243
    }

    inline bool isRankDof(DegreeOfFreedom dof)
    {
244
      TEST_EXIT_DBG(rankDofs)("No rank dofs set!\n");
245

246
      return (*rankDofs)[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
247
248
249
    }
#endif

250
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
251
    ///
252
253
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
254
    ///
255
    std::string name;
256

Thomas Witkowski's avatar
Thomas Witkowski committed
257
    ///
258
    ElementVector elementVector;
259

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
266
    ///
267
    std::vector<double*> operatorEstFactor;
268

Thomas Witkowski's avatar
Thomas Witkowski committed
269
    ///
270
271
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
272
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
273
    int nBasFcts;
274

Thomas Witkowski's avatar
Thomas Witkowski committed
275
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
276
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
277
278

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
279
    std::map<DegreeOfFreedom, bool> *rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
280
#endif
281

Thomas Witkowski's avatar
Thomas Witkowski committed
282
  };
283

Thomas Witkowski's avatar
Thomas Witkowski committed
284
  /// Specifies which operation should be done after coarsening
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
  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)
316
      {}
317
318
319

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

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

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

Thomas Witkowski's avatar
Thomas Witkowski committed
350
    /// Constructs a DOFVector with name n belonging to FiniteElemSpace f
351
    DOFVector(const FiniteElemSpace* f, std::string n); 
352

Thomas Witkowski's avatar
Thomas Witkowski committed
353
    /// Initialization.
354
    void init(const FiniteElemSpace* f, std::string n);
355

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

Thomas Witkowski's avatar
Thomas Witkowski committed
365
    /// Destructor
366
367
    virtual ~DOFVector();

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
391
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
    inline CoarsenOperation getCoarsenOperation() 
    { 
394
      return coarsenOperation; 
395
    }
396

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

400
    /// Interpolation after refinement. Implemented for DOFVector<double>
401
    inline void refineInterpol(RCNeighbourList&, int) {}
402

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
446
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
    inline T& operator[](DegreeOfFreedom i) 
    {
449
      FUNCNAME("DOFVector<T>::operator[]");
450
451

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

454
      return vec[i];
455
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
456
 
Thomas Witkowski's avatar
Thomas Witkowski committed
457
    /// Calculates Integral of this DOFVector
458
459
    double Int(Quadrature* q = NULL) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
460
461
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
462
463
464
465
466
467
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
470
471
    
    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
472
473
474
475
476
477
478
    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
479
    /// Calculates L1 norm of this DOFVector
480
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
481
 
Thomas Witkowski's avatar
Thomas Witkowski committed
482
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
483
484
    inline double L2Norm(Quadrature* q = NULL) const 
    {
485
      return sqrt(L2NormSquare());
486
    }
487

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

Thomas Witkowski's avatar
Thomas Witkowski committed
491
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
492
493
    inline double H1Norm(Quadrature* q = NULL) const 
    {
494
      return sqrt(H1NormSquare());
495
    }
496

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

Thomas Witkowski's avatar
Thomas Witkowski committed
500
    /// Calculates euclidian norm of this DOFVector
501
502
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
503
    /// Returns square of the euclidian norm.
504
505
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
506
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
507
508
    inline double l2norm() const 
    { 
509
      return nrm2();
510
    }
511

Thomas Witkowski's avatar
Thomas Witkowski committed
512
    /// Calculates the absolute sum of this DOFVector
513
514
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
515
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
516
517
    inline double l1norm() const 
    { 
518
      return asum();
519
    } 
520

Thomas Witkowski's avatar
Thomas Witkowski committed
521
    /// Calculates doublewell of this DOFVector
522
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
523
 
Thomas Witkowski's avatar
Thomas Witkowski committed
524
    /// Calculates the sum of this DOFVector
525
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
526
 
Thomas Witkowski's avatar
Thomas Witkowski committed
527
    /// Sets \ref vec[i] = val, i=0 , ... , size
528
529
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
530
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
531
532
    inline DOFVector<T>& operator=(T d) 
    {
533
534
      set(d); 
      return *this;
535
    } 
536

Thomas Witkowski's avatar
Thomas Witkowski committed
537
    /// Assignment operator between two vectors
538
    DOFVector<T>& operator=(const DOFVector<T>&); 
539

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

543
    /// Returns minimum of DOFVector.
544
545
    T min() const; 

546
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
547
548
    T max() const;

549
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
550
    T absMax() const;
551

552
553
554
555
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Used by interpol while mesh traversal.
556
    void interpolFct(ElInfo* elinfo);
557

558
    /// Prints \ref vec to stdout.
559
560
    void print() const; 

561
    ///
562
563
    int calcMemoryUsage() const;

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

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

570

Thomas Witkowski's avatar
Thomas Witkowski committed
571
572
573
574
575
    /// 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>
    inline const T evalAtPoint(WorldVector<double> &p, 
			       ElInfo *oldElInfo = NULL) const 
576
577
578
579
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
    }
580

581
582
583
    const bool getDOFidxAtPoint(WorldVector<double> &p, DegreeOfFreedom &idx, ElInfo *oldElInfo = NULL, bool useOldElInfo = false);


584
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
585
586
    void serialize(std::ostream &out) 
    {
587
588
589
      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));
590
    }
591

592
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
593
594
    void deserialize(std::istream &in) 
    {
595
596
597
598
      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));
599
    }
600

601
602
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
603
604
605

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

606
607
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
608
609

  protected: 
610

611
    /// Data container
612
    std::vector<T> vec; 
613
 
614
    /// Specifies what operation should be performed after coarsening
615
616
    CoarsenOperation coarsenOperation;

617
    /// Used by \ref interpol
618
619
    AbstractFunction<T, WorldVector<double> > *interFct;

620
    /// Used for mesh traversal
621
622
623
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
624

625
  template<>
626
627
628
629
630
631
632
633
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

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

  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
634
635
  const double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
					      ElInfo *oldElInfo) const;
636
637

  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
638
639
640
  const WorldVector<double> 
    DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p,
 						 ElInfo *oldElInfo) const;
641

642
  template<>
643
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
644
645

  template<>
646
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
647

Thomas Witkowski's avatar
Thomas Witkowski committed
648
649
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
650
    return v.min();
651
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
652

Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
655
    return v.max();
656
  }
657
658
659
660
661
662
663
664
665
666


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Thomas Witkowski's avatar
Thomas Witkowski committed
667
668
    /// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
    /// as DOFContainer at DOFAdmin
669
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
670
671
672
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
673
    }
674
  
675
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
    ~DOFVectorDOF() 
    {
678
      feSpace->getAdmin()->removeDOFContainer(this);
679
    }
680

Thomas Witkowski's avatar
Thomas Witkowski committed
681
682
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
683
684
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
685
      return DOFVector<DegreeOfFreedom>::operator[](i);
686
    }
687

Thomas Witkowski's avatar
Thomas Witkowski committed
688
689
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
690
691
692
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

693
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
694
695
    int getSize() const 
    {
696
      return DOFVector<DegreeOfFreedom>::getSize();
697
    }
698

699
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
700
701
    void resize(int size) 
    {
702
703
704
705
706
707
708
709
710
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
711

712
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
713
714
  double norm(DOFVector<T> *vec) 
  {
715
    return vec->nrm2();
716
  }
717
718

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
719
720
  double L2Norm(DOFVector<T> *vec) 
  {
721
    return vec->L2Norm();
722
  }
723
724

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
731
  void print(DOFVector<T> *vec) 
732
  {
733
    vec->print();
734
  }
735
736
737
738
739
740
741
742
743
744
745
746
747
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>
  void axpy(double a,const DOFVector<T>& x, DOFVector<T>& y);

  // 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
781
782
  inline void scal(T a, DOFVector<T>& y) 
  {
783
784
    y *= a;
  }
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800

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

801
802
803
804
805
806
807
808
809
  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);


810
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
811
812
  inline void set(DOFVector<T>& vec, T d) 
  {
813
    vec.set(d);
814
  }
815
816

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
817
818
  inline void setValue(DOFVector<T>& vec, T d) 
  {
819
    vec.set(d);
820
  }
821
822

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
823
824
  inline int size(DOFVector<T> *vec) 
  {
825
    return vec->getUsedSize();
826
  }
827

Thomas Witkowski's avatar
Thomas Witkowski committed
828
829
830
831
832
833
834
835
836
837
  template<typename T>
  inline void checkFeSpace(const FiniteElemSpace* feSpace, const std::vector<T>& vec)
  {
    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());
  }

838
839
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
840
841
842
843
844
  
  template<typename T>
  std::vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     std::vector<DOFVector<double>*> *res);
  
845
  /** \brief
Thomas Witkowski's avatar
Thomas Witkowski committed
846
847
848
   * Computes the integral of a function that includes two different DOFVectors.
   * This function works also for the case that the DOFVectors are defined on
   * two different meshes.
849
850
851
852
853
   */
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct);

Thomas Witkowski's avatar
Thomas Witkowski committed
854
855
  /// Computes the integral of a function with two DOFVectors defined on the 
  /// same mesh.
856
857
858
859
  double intSingleMesh(const DOFVector<double> &vec1,
		       const DOFVector<double> &vec2,
		       BinaryAbstractFunction<double, double, double> *fct);

Thomas Witkowski's avatar
Thomas Witkowski committed
860
861
  /// Computes the integral of a function with two DOFVectors defined on 
  /// different meshes.
862
863
864
865
  double intDualMesh(const DOFVector<double> &vec1,
		     const DOFVector<double> &vec2,
		     BinaryAbstractFunction<double, double, double> *fct);

Thomas Witkowski's avatar
Thomas Witkowski committed
866
867
868
  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct);

869
870
  double integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<double, WorldVector<double> > *fct);
871
872
873
874
875
876
877
  
  double integrateGeneral(const std::vector<DOFVector<double>*> &vecs,
			  AbstractFunction<double, std::vector<double> > *fct);
  
  double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
				  const std::vector<DOFVector<double>*> &grds,
				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > *fct);
878
879
880
881
882
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_