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

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

598
599
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
600
601
602

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

603
604
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
605
606

  protected: 
607

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

614
    /// Used by \ref interpol
615
616
    AbstractFunction<T, WorldVector<double> > *interFct;

617
    /// Used for mesh traversal
618
619
620
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
621

622
  template<>
623
624
625
626
627
628
629
630
  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
631
632
  const double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
					      ElInfo *oldElInfo) const;
633
634

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

639
  template<>
640
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
641
642

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

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

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


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

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
708

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

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
728
  void print(DOFVector<T> *vec) 
729
  {
730
    vec->print();
731
  }
732
733
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

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

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

798
799
800
801
802
803
804
805
806
  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);


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
825
826
827
828
829
830
831
832
833
834
  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());
  }

835
836
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
837
838
839
840
841
  
  template<typename T>
  std::vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     std::vector<DOFVector<double>*> *res);
  
842
  /** \brief
Thomas Witkowski's avatar
Thomas Witkowski committed
843
844
845
   * 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.
846
847
848
849
850
   */
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct);

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
863
864
865
  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct);

866
867
868
869
  double integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<double, WorldVector<double> > *fct);


870
871
872
873
874
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_