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


20
21
22
23
24
25

/** \file DOFVector.h */

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

30
#include "AMDiS_fwd.h"
31
32
33
34
35
36
37
38
39
40
41
#include "FixVec.h"
#include "Global.h" 
#include "Flag.h" 
#include "RCNeighbourList.h" 
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "DOFContainer.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
Thomas Witkowski's avatar
Thomas Witkowski committed
42
#include "BasisFunction.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
43
#include "FiniteElemSpace.h"
44
#include "SurfaceQuadrature.h"
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
70
71
    /** \brief 
     * For the given element, this function returns an array of all DOFs of this
     * DOFVector that are defined on this element.
     */
72
73
    virtual void getLocalVector(const Element *el, 
				mtl::dense_vector<T>& localVec) const;
74

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

84
85
86
87
88
    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
89

90
91
    const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo, 
				      const Quadrature *quad,
92
				      const FastQuadrature *quadFast,
93
				      WorldVector<T> *grdAtQPs) const;
94

Thomas Witkowski's avatar
Thomas Witkowski committed
95
96
97
98
99
100
    const WorldVector<T> *getGrdAtQPs(const ElInfo *smallElInfo, 
				      const ElInfo *largeElInfo,
				      const Quadrature *quad,
				      const FastQuadrature *quadFast,
				      WorldVector<T> *grdAtQPs) const;

101
102
    const WorldMatrix<T> *getD2AtQPs(const ElInfo *elInfo, 
				     const Quadrature *quad,
103
				     const FastQuadrature *quadFast,
104
				     WorldMatrix<T> *d2AtQPs) const;
105

Thomas Witkowski's avatar
Thomas Witkowski committed
106
    /// Returns the FE space the DOF vector is defined on.
107
    inline const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
108
    {
109
      return feSpace;
110
    }
111

Thomas Witkowski's avatar
Thomas Witkowski committed
112
113
114
115
116
117
118
119
120
121
122
123
124
    /** \brief
     * Assembles the element vector for the given ellement and adds the
     * element matrix to the current DOF vector.
     */ 
    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);
125
126

    void addElementVector(T sign,
127
			  const ElementVector& elVec, 
128
			  const BoundaryType *bound,
129
			  ElInfo *elInfo,
130
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
131
132
133
134
135
136
137

    /* \brief
     * 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.
     */
    void finishAssembling();
138
139
140
141
142
143
144
145
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
146
    }
147

Thomas Witkowski's avatar
Thomas Witkowski committed
148
149
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
150
      return operatorFactor.begin();
151
    }
152

Thomas Witkowski's avatar
Thomas Witkowski committed
153
154
    inline std::vector<double*>::iterator getOperatorFactorEnd() 
    {
155
      return operatorFactor.end();
156
    }
157

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

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


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

Thomas Witkowski's avatar
Thomas Witkowski committed
174
175
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
176
      return operators.end();
177
    }
178
179
180
181
182
183
184
185
186
187
188

    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
189
190
    inline std::vector<Operator*>& getOperators() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
191
      return operators; 
192
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
193

Thomas Witkowski's avatar
Thomas Witkowski committed
194
195
    inline std::vector<double*>& getOperatorFactor() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
196
      return operatorFactor; 
197
    }
198

Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
    inline std::vector<double*>& getOperatorEstFactor() 
    { 
201
      return operatorEstFactor; 
202
    }
203

Thomas Witkowski's avatar
Thomas Witkowski committed
204
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
205
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
206
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
207
      return name; 
208
    } 
209

210
211
    inline void setName(std::string n)
    {
212
      name = n;
213
214
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
215
216
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
217
      return boundaryManager; 
218
    }
219

Thomas Witkowski's avatar
Thomas Witkowski committed
220
221
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
222
      boundaryManager = bm;
223
    }
224

Thomas Witkowski's avatar
Thomas Witkowski committed
225
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
226
    inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
227
    {
228
      rankDofs = &dofmap;
Thomas Witkowski's avatar
Thomas Witkowski committed
229
230
231
232
    }

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

235
      return (*rankDofs)[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
238
    }
#endif

239
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
240
    ///
241
242
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
243
    ///
244
    std::string name;
245

Thomas Witkowski's avatar
Thomas Witkowski committed
246
    ///
247
    ElementVector elementVector;
248

Thomas Witkowski's avatar
Thomas Witkowski committed
249
    ///
250
    std::vector<Operator*> operators;
251

Thomas Witkowski's avatar
Thomas Witkowski committed
252
    ///
253
    std::vector<double*> operatorFactor;
254

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

Thomas Witkowski's avatar
Thomas Witkowski committed
258
    ///
259
260
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
261
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
262
    int nBasFcts;
263

Thomas Witkowski's avatar
Thomas Witkowski committed
264
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
265
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
266
267

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
268
    std::map<DegreeOfFreedom, bool> *rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
269
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
270
  };
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
  /// Specifies which operation should be done after coarsening
273
274
275
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
  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)
304
      {}
305
306
307

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
308
      {}
309
310
311
312
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
313
314
315
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
316

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
341
    /// Initialization.
342
    void init(const FiniteElemSpace* f, std::string n);
343

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

Thomas Witkowski's avatar
Thomas Witkowski committed
353
    /// Destructor
354
355
    virtual ~DOFVector();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
362
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
363
364
    typename std::vector<T>::iterator end() 
    { 
365
      return vec.end(); 
366
    }
367
368
369
370
371
372
  
    /** \brief
     * Used by DOFAdmin to compress this DOFVector. Implementation of
     * DOFIndexedBase::compress()
     */
    virtual void compressDOFIndexed(int first, int last,
373
				    std::vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
374
375

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
376
377
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
378
      coarsenOperation = op; 
379
    }
380

Thomas Witkowski's avatar
Thomas Witkowski committed
381
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
    inline CoarsenOperation getCoarsenOperation() 
    { 
384
      return coarsenOperation; 
385
    }
386

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

390
    /// Interpolation after refinement. Implemented for DOFVector<double>
391
    inline void refineInterpol(RCNeighbourList&, int) {}
392

Thomas Witkowski's avatar
Thomas Witkowski committed
393
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
394
395
    std::vector<T>& getVector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
396
      return vec;
397
    }
398

Thomas Witkowski's avatar
Thomas Witkowski committed
399
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
    inline int getSize() const 
    { 
402
      return vec.size();
403
    } 
404

Thomas Witkowski's avatar
Thomas Witkowski committed
405
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
406
407
    inline int getUsedSize() const 
    { 
408
      return this->feSpace->getAdmin()->getUsedSize(); 
409
    }
410

Thomas Witkowski's avatar
Thomas Witkowski committed
411
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
    inline void resize(int n) 
    { 
414
      FUNCNAME("DOFVector<T>::resize()");
415
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
416
      vec.resize(n);
417
    } 
418

Thomas Witkowski's avatar
Thomas Witkowski committed
419
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
420
421
    inline void resize(int n, T init) 
    { 
422
      FUNCNAME("DOFVector<T>::resize()");
423
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
424
      vec.resize(n, init);
425
    } 
426

Thomas Witkowski's avatar
Thomas Witkowski committed
427
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
428
429
    inline const T& operator[](DegreeOfFreedom i) const 
    {
430
      FUNCNAME("DOFVector<T>::operator[]");
431
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
432
	("Illegal vector index %d.\n", i);
433
      return vec[i];
434
    } 
435

Thomas Witkowski's avatar
Thomas Witkowski committed
436
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
437
438
    inline T& operator[](DegreeOfFreedom i) 
    {
439
      FUNCNAME("DOFVector<T>::operator[]");
440
441

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

444
      return vec[i];
445
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
446
 
Thomas Witkowski's avatar
Thomas Witkowski committed
447
    /// Calculates Integral of this DOFVector
448
449
    double Int(Quadrature* q = NULL) const;

450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
    /** \brief
     * Calculates Integral of this DOFVector over parts of the domain
     * boundary, indicated by boundaryType. Implemented for DOFVector<double>
     **/
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }

    /** \brief
     * Calculates Integral of this DOFVector times normal vector over parts 
     * of the domain boundary, indicated by boundaryType. Implemented for 
     * DOFVector<WorldVector<double> >
     **/
    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
473
    /// Calculates L1 norm of this DOFVector
474
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
475
 
Thomas Witkowski's avatar
Thomas Witkowski committed
476
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
477
478
    inline double L2Norm(Quadrature* q = NULL) const 
    {
479
      return sqrt(L2NormSquare());
480
    }
481

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

Thomas Witkowski's avatar
Thomas Witkowski committed
485
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
486
487
    inline double H1Norm(Quadrature* q = NULL) const 
    {
488
      return sqrt(H1NormSquare());
489
    }
490

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

Thomas Witkowski's avatar
Thomas Witkowski committed
494
    /// Calculates euclidian norm of this DOFVector
495
496
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
497
    /// Returns square of the euclidian norm.
498
499
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
500
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
501
502
    inline double l2norm() const 
    { 
503
      return nrm2();
504
    }
505

Thomas Witkowski's avatar
Thomas Witkowski committed
506
    /// Calculates the absolute sum of this DOFVector
507
508
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
509
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
510
511
    inline double l1norm() const 
    { 
512
      return asum();
513
    } 
514

Thomas Witkowski's avatar
Thomas Witkowski committed
515
    /// Calculates doublewell of this DOFVector
516
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
517
 
Thomas Witkowski's avatar
Thomas Witkowski committed
518
    /// Calculates the sum of this DOFVector
519
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
520
 
Thomas Witkowski's avatar
Thomas Witkowski committed
521
    /// Sets \ref vec[i] = val, i=0 , ... , size
522
523
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
524
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
525
526
    inline DOFVector<T>& operator=(T d) 
    {
527
528
      set(d); 
      return *this;
529
    } 
530

Thomas Witkowski's avatar
Thomas Witkowski committed
531
    /// Assignment operator between two vectors
532
    DOFVector<T>& operator=(const DOFVector<T>&); 
533

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

537
    /// Returns minimum of DOFVector.
538
539
    T min() const; 

540
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
541
542
    T max() const;

543
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
544
    T absMax() const;
545

546
547
548
549
    /// Returns the average value of the DOFVector.
    T average() const;

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

552
    /// Prints \ref vec to stdout.
553
554
    void print() const; 

555
    ///
556
557
    int calcMemoryUsage() const;

558
559
560
561
562
563
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

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

566

567
568
569
570
    /** 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 >
     */
571
572
    inline const T& evalAtPoint(
      WorldVector<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const 
573
574
575
576
577
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
      return *value;
    }
578

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

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

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

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

600
    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
601
602

  protected: 
603

604
    /// Data container
605
    std::vector<T> vec; 
606
 
607
    /// Specifies what operation should be performed after coarsening
608
609
    CoarsenOperation coarsenOperation;

610
    /// Used by \ref interpol
611
612
    AbstractFunction<T, WorldVector<double> > *interFct;

613
    /// Used for mesh traversal
614
615
616
    static DOFVector<T> *traverseVector;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
617

618
  template<>
619
620
621
622
623
624
625
626
627
628
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

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

  template<>
  const double& DOFVector<double>::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
629
630

  template<>
631
632
  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
633

634
  template<>
635
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
636
637

  template<>
638
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
639

Thomas Witkowski's avatar
Thomas Witkowski committed
640
641
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
642
    return v.min();
643
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
644

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


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

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

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

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

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

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
707

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

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

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

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

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

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

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


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

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

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

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

834
835
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855

  /** \brief
   * 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.
   */
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct);

  /// Computes the integral of a function with two DOFVectors defined on the same mesh.
  double intSingleMesh(const DOFVector<double> &vec1,
		       const DOFVector<double> &vec2,
		       BinaryAbstractFunction<double, double, double> *fct);

  /// Computes the integral of a function with two DOFVectors defined on different meshes.
  double intDualMesh(const DOFVector<double> &vec1,
		     const DOFVector<double> &vec2,
		     BinaryAbstractFunction<double, double, double> *fct);

Thomas Witkowski's avatar
Thomas Witkowski committed
856
857
858
  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct);

859
860
861
862
  double integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<double, WorldVector<double> > *fct);


863
864
865
866
867
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_