DOFVector.h 26.2 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

Praetorius, Simon's avatar
Praetorius, Simon 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

Praetorius, Simon's avatar
Praetorius, Simon 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

Praetorius, Simon's avatar
Praetorius, Simon committed
90
91
92
93
    /** \brief
     * Evaluates the gradient of a DOF vector at a set of quadrature points defined on the
     * given element.
     */
94
95
96
97
98
99
100
101
102
103
104
    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;

Praetorius, Simon's avatar
Praetorius, Simon committed
105
106
107
108
    /** \brief
     * Evaluates the comp'th component of the derivative of a DOF vector at a set of quadrature points defined on the
     * given element.
     */
109
110
111
112
113
114
115
116
117
118
119
120
121
    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;

Praetorius, Simon's avatar
Praetorius, Simon committed
122
123
124
125
    /** \brief
     * Evaluates the jacobian of a DOF vector at a set of quadrature points defined on the
     * given element.
     */
126
127
128
129
    void getD2AtQPs(const ElInfo *elInfo,
		    const Quadrature *quad,
		    const FastQuadrature *quadFast,
		    mtl::dense_vector<typename D2Type<T>::type> &d2AtQPs) const;
130

Thomas Witkowski's avatar
Thomas Witkowski committed
131
    /// Returns the FE space the DOF vector is defined on.
132
    inline const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
133
    {
134
      return feSpace;
135
    }
136

Praetorius, Simon's avatar
Praetorius, Simon committed
137
138
139
140
    /** \brief
     * 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
141
142
143
144
145
146
147
148
149
    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);
150
151

    void addElementVector(T sign,
152
			  const ElementVector& elVec, 
153
			  const BoundaryType *bound,
154
			  ElInfo *elInfo,
155
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
156

Praetorius, Simon's avatar
Praetorius, Simon committed
157
158
159
160
161
    /* \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.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
162
    void finishAssembling();
163
164
165
166
167
168
169
170
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
171
    }
172

Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
175
      return operatorFactor.begin();
176
    }
177

Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
    inline std::vector<double*>::iterator getOperatorFactorEnd() 
    {
180
      return operatorFactor.end();
181
    }
182

Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
185
      return operatorEstFactor.begin();
186
    }
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188
189
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() 
    {
190
      return operatorEstFactor.end();
191
    }
192
193


Thomas Witkowski's avatar
Thomas Witkowski committed
194
195
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
196
      return operators.begin();
197
    }
198

Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
201
      return operators.end();
202
    }
203
204
205
206
207
208
209
210
211
212
213

    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
214
215
    inline std::vector<Operator*>& getOperators() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
216
      return operators; 
217
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
218

Thomas Witkowski's avatar
Thomas Witkowski committed
219
220
    inline std::vector<double*>& getOperatorFactor() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
221
      return operatorFactor; 
222
    }
223

Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
    inline std::vector<double*>& getOperatorEstFactor() 
    { 
226
      return operatorEstFactor; 
227
    }
228

Thomas Witkowski's avatar
Thomas Witkowski committed
229
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
230
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
231
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
232
      return name; 
233
    } 
234

235
236
    inline void setName(std::string n)
    {
237
      name = n;
238
239
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
240
241
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
242
      return boundaryManager; 
243
    }
244

Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
247
      boundaryManager = bm;
248
    }
249

Thomas Witkowski's avatar
Thomas Witkowski committed
250
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
251
    inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
252
    {
253
      rankDofs = &dofmap;
Thomas Witkowski's avatar
Thomas Witkowski committed
254
255
256
257
    }

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

260
      return (*rankDofs)[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
263
    }
#endif

264
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
265
    ///
266
267
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    ///
269
    std::string name;
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    ///
272
    ElementVector elementVector;
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
    ///
275
    std::vector<Operator*> operators;
276

Thomas Witkowski's avatar
Thomas Witkowski committed
277
    ///
278
    std::vector<double*> operatorFactor;
279

Thomas Witkowski's avatar
Thomas Witkowski committed
280
    ///
281
    std::vector<double*> operatorEstFactor;
282

Thomas Witkowski's avatar
Thomas Witkowski committed
283
    ///
284
285
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
286
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
287
    int nBasFcts;
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
290
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
291
292

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
293
    std::map<DegreeOfFreedom, bool> *rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
294
#endif
295

Thomas Witkowski's avatar
Thomas Witkowski committed
296
  };
297

Thomas Witkowski's avatar
Thomas Witkowski committed
298
  /// Specifies which operation should be done after coarsening
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
  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)
330
      {}
331
332
333

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
334
      {}
335
336
337
338
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
339
340
341
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
342

Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
      DOFVector<T> *create() 
      { 
345
	return new DOFVector<T>(feSpace, ""); 
346
      }
347

Thomas Witkowski's avatar
Thomas Witkowski committed
348
349
      void free(DOFVector<T> *vec) 
      { 
350
	delete vec; 
351
      }
352
353
354
355
356
357

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
358
    /// Empty constructor. No initialization!
359
360
361
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
362
    {}
363

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

Thomas Witkowski's avatar
Thomas Witkowski committed
367
    /// Initialization.
368
    void init(const FiniteElemSpace* f, std::string n);
369

Thomas Witkowski's avatar
Thomas Witkowski committed
370
    /// Copy Constructor
Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
    DOFVector(const DOFVector& rhs) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
373
      *this = rhs;   
374
375
376
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
377
    }
378

Thomas Witkowski's avatar
Thomas Witkowski committed
379
    /// Destructor
380
381
    virtual ~DOFVector();

Thomas Witkowski's avatar
Thomas Witkowski committed
382
    /// Returns iterator to the begin of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
383
384
    typename std::vector<T>::iterator begin() 
    { 
385
      return vec.begin(); 
386
    }
387

Thomas Witkowski's avatar
Thomas Witkowski committed
388
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
    typename std::vector<T>::iterator end() 
    { 
391
      return vec.end(); 
392
    }
393
  
Praetorius, Simon's avatar
Praetorius, Simon committed
394
395
396
397
    /** \brief
     * Used by DOFAdmin to compress this DOFVector. Implementation of
     * DOFIndexedBase::compress()
     */
398
    virtual void compressDOFIndexed(int first, int last,
399
				    std::vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401

    /// Sets \ref coarsenOperation to op
Thomas Witkowski's avatar
Thomas Witkowski committed
402
403
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
404
      coarsenOperation = op; 
405
    }
406

Thomas Witkowski's avatar
Thomas Witkowski committed
407
    /// Returns \ref coarsenOperation
Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
    inline CoarsenOperation getCoarsenOperation() 
    { 
410
      return coarsenOperation; 
411
    }
412

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

416
    /// Interpolation after refinement. Implemented for DOFVector<double>
417
    inline void refineInterpol(RCNeighbourList&, int) {}
418

Thomas Witkowski's avatar
Thomas Witkowski committed
419
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
420
421
    std::vector<T>& getVector() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
422
      return vec;
423
    }
424

Thomas Witkowski's avatar
Thomas Witkowski committed
425
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
426
427
    inline int getSize() const 
    { 
428
      return vec.size();
429
    } 
430

Thomas Witkowski's avatar
Thomas Witkowski committed
431
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
    inline int getUsedSize() const 
    { 
434
      return this->feSpace->getAdmin()->getUsedSize(); 
435
    }
436

Thomas Witkowski's avatar
Thomas Witkowski committed
437
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
438
439
    inline void resize(int n) 
    { 
440
      FUNCNAME("DOFVector<T>::resize()");
441
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
442
      vec.resize(n);
443
    } 
444

Thomas Witkowski's avatar
Thomas Witkowski committed
445
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
    inline void resize(int n, T init) 
    { 
448
      FUNCNAME("DOFVector<T>::resize()");
449
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
450
      vec.resize(n, init);
451
    } 
452

Thomas Witkowski's avatar
Thomas Witkowski committed
453
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
454
455
    inline const T& operator[](DegreeOfFreedom i) const 
    {
456
      FUNCNAME("DOFVector<T>::operator[]");
457
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
458
	("Illegal vector index %d.\n", i);
459
      return vec[i];
460
    } 
461

Thomas Witkowski's avatar
Thomas Witkowski committed
462
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
    inline T& operator[](DegreeOfFreedom i) 
    {
465
      FUNCNAME("DOFVector<T>::operator[]");
466
467

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

470
      return vec[i];
471
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
472
 
Thomas Witkowski's avatar
Thomas Witkowski committed
473
    /// Calculates Integral of this DOFVector
474
475
    double Int(Quadrature* q = NULL) const;

Praetorius, Simon's avatar
Praetorius, Simon committed
476
477
478
479
    /** \brief
     * Calculates Integral of this DOFVector over parts of the domain
     * boundary, indicated by boundaryType. Implemented for DOFVector<double>
     **/
480
481
482
483
484
485
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
486
487
488
489
490
491

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

589

Praetorius, Simon's avatar
Praetorius, Simon committed
590
591
592
593
594
595
    /** 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, T* value = NULL) const 
596
597
598
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
Praetorius, Simon's avatar
Praetorius, Simon committed
599
      return *value;
600
    }
601

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

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

615
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
616
617
    void deserialize(std::istream &in) 
    {
618
619
620
621
      unsigned int size;
      in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int));
      vec.resize(size);
      in.read(reinterpret_cast<char*>(&(vec[0])), size * sizeof(T));
622
    }
623

Praetorius, Simon's avatar
Praetorius, Simon committed
624
//     DOFVector<WorldVector<T> > *getGradient(DOFVector<WorldVector<T> >*) const;
625
626
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
627
628
629

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

Praetorius, Simon's avatar
Praetorius, Simon committed
630
//     DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
631
632
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
633
634

  protected: 
635

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

Thomas Witkowski's avatar
Thomas Witkowski committed
643

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

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

  template<>
Praetorius, Simon's avatar
Praetorius, Simon committed
653
654
  const double& DOFVector<double>::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
655
656

  template<>
Praetorius, Simon's avatar
Praetorius, Simon committed
657
658
  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
659

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

  template<>
664
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
665

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

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


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

Praetorius, Simon's avatar
Praetorius, Simon committed
701
702
703
704
    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVector<DegreeOfFreedom>::operator[]()
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
705
706
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
707
      return DOFVector<DegreeOfFreedom>::operator[](i);
708
    }
709

Thomas Witkowski's avatar
Thomas Witkowski committed
710
711
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
712
713
714
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

715
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
    int getSize() const 
    {
718
      return DOFVector<DegreeOfFreedom>::getSize();
719
    }
720

721
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
722
723
    void resize(int size) 
    {
724
725
726
727
728
729
730
731
732
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
733

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
  double H1Norm(DOFVector<T> *vec) 
  {
749
    return vec->H1Norm();
750
  }
751
752

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

  // 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>
789
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
790
791
792
793
794
795
796
797
798
799
800
801
802

  // 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
803
804
  inline void scal(T a, DOFVector<T>& y) 
  {
805
806
    y *= a;
  }
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822

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

823
824
825
826
827
828
829
830
831
  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);


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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
839
840
  inline void setValue(DOFVector<T>& vec, T d) 
  {
841
    vec.set(d);
842
  }
843
844

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
845
846
  inline int size(DOFVector<T> *vec) 
  {
847
    return vec->getUsedSize();
848
  }
849

850
851
852
853
854
855
  template<typename T>
  inline int size(const DOFVector<T>& vec) 
  {
    return vec.getUsedSize();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
856
857
858
859
860
861
862
863
864
865
  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());
  }

866
867
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
868
869
870
871
  
  template<typename T>
  std::vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     std::vector<DOFVector<double>*> *res);
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896

  
  /// Computes the integral: \f$ \int f(\vec{x})\,\text{d}\vec{x}\f$
  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct);
  
  template<typename TOut>
  TOut integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
  {
    return integrate_Coords(feSpace, fct);
  }
  
  /// Computes the integral: \f$ \int f(v(\vec{x}))\,\text{d}\vec{x}\f$
  template<typename TOut, typename T>
  TOut integrate_Vec(const DOFVector<T> &vec,
		  AbstractFunction<TOut, T> *fct);
  
  template<typename TOut, typename T>
  TOut integrate(const DOFVector<T> &vec,
		  AbstractFunction<TOut, T> *fct)
  {
    return integrate_Vec(vec, fct);
  }
897
  
898
  /** \brief
899
900
901
   * Computes the integral of a function that includes two different DOFVectors:
   * \f$ \int f(v(\vec{x}), w(\vec{x}))\,\text{d}\vec{x}\f$
   * This function works also for the case that the DOFVectors are defined on two different
Praetorius, Simon's avatar
Praetorius, Simon committed
902
   * meshes.
903
   */
904
905
906
907
908
909
910
911
912
913
914
915
916
  template<typename TOut, typename T1, typename T2>
  TOut integrate_Vec2(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct);

  template<typename TOut, typename T1, typename T2>
  TOut integrate(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct)
  {
    return integrate_Vec2(vec1, vec2, fct);
  }
  
Praetorius, Simon's avatar
Praetorius, Simon committed
917
  /// Computes the integral of a function with two DOFVectors defined on the same mesh.
918
919
920
921
  template<typename TOut, typename T1, typename T2>
  TOut int_Vec2_SingleMesh(const DOFVector<T1> &vec1,
		       const DOFVector<T2> &vec2,
		       BinaryAbstractFunction<TOut, T1, T2> *fct);
922

Praetorius, Simon's avatar
Praetorius, Simon committed
923
  /// Computes the integral of a function with two DOFVectors defined on different meshes.
924
925
926
927
928
929
930
931
932
933
934
935
936
937
  template<typename TOut, typename T1, typename T2>
  TOut int_Vec2_DualMesh(const DOFVector<T1> &vec1,
		     const DOFVector<T2> &vec2,
		     BinaryAbstractFunction<TOut, T1, T2> *fct);
  
  /// Computes the integral: \f$ \int v(\vec{x}) f(\vec{x})\,\text{d}\vec{x}\f$
  template<typename T1, typename T2>
  typename ProductType<T1,T2>::type integrate_VecTimesCoords(const DOFVector<T1> &vec,
		   AbstractFunction<T2, WorldVector<double> > *fct);
  
  /// Computes the integral: \f$ \int f(v(\vec{x}), \vec{x})\,\text{d}\vec{x}\f$
  template<typename TOut, typename T>
  TOut integrate_VecAndCoords(const DOFVector<T> &vec,
			   BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct);
938
  
939
940
941
942
943
944
945
946
  template<typename TOut, typename T>
  TOut integrate_VecAndCoords(const DOFVector<T> &vec,
			   BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct)
  {
    return integrate(vec, fct);
  }
  
  /// Computes the integral: \f$ \int f(\{v_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$
947
948
949
  double integrateGeneral(const std::vector<DOFVector<double>*> &vecs,
			  AbstractFunction<double, std::vector<double> > *fct);
  
950
  /// Computes the integral: \f$ \int f(\{v_i(\vec{x})\}_i,\{\nabla w_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$
951
952
  double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
				  const std::vector<DOFVector<double>*> &grds,
Praetorius, Simon's avatar
Praetorius, Simon committed
953
				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > > *fct);
954

955
956
957
958
959
}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_