DOFVector.h 26.8 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
Thomas Witkowski's avatar
Thomas Witkowski committed
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
    double Int(Quadrature* q = NULL) const
    {
      return Int(-1, q);
    }

    /** \brief
     * Calculates Integral of this DOFVector. 
     *
     * \param[in]  meshlevel  Then mesh level on which the integral should be
     *                        calculated. Usually, only -1 for the leaf level
     *                        makes sence, but other values can be used, e.g.,
     *                        for debugging.
     * \param[in]  q          Quadrature object. If not specified, the function
     *                        creates a new quadrature object.
     */
    double Int(int meshLevel, Quadrature* q = NULL) const;

491

Praetorius, Simon's avatar
Praetorius, Simon committed
492
493
494
495
    /** \brief
     * Calculates Integral of this DOFVector over parts of the domain
     * boundary, indicated by boundaryType. Implemented for DOFVector<double>
     **/
496
497
498
499
500
501
    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
502
503
504
505
506
507

    /** \brief
     * Calculates Integral of this DOFVector times normal vector over parts 
     * of the domain boundary, indicated by boundaryType. Implemented for 
     * DOFVector<WorldVector<double> >
     **/
508
509
510
511
512
513
514
    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
515
    /// Calculates L1 norm of this DOFVector
516
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
517
 
Thomas Witkowski's avatar
Thomas Witkowski committed
518
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
519
520
    inline double L2Norm(Quadrature* q = NULL) const 
    {
521
      return sqrt(L2NormSquare());
522
    }
523

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

Thomas Witkowski's avatar
Thomas Witkowski committed
527
    /// Calculates H1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
528
529
    inline double H1Norm(Quadrature* q = NULL) const 
    {
530
      return sqrt(H1NormSquare());
531
    }
532

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

Thomas Witkowski's avatar
Thomas Witkowski committed
536
    /// Calculates euclidian norm of this DOFVector
537
538
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
539
    /// Returns square of the euclidian norm.
540
541
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
542
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
543
544
    inline double l2norm() const 
    { 
545
      return nrm2();
546
    }
547

Thomas Witkowski's avatar
Thomas Witkowski committed
548
    /// Calculates the absolute sum of this DOFVector
549
550
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
551
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
552
553
    inline double l1norm() const 
    { 
554
      return asum();
555
    } 
556

Thomas Witkowski's avatar
Thomas Witkowski committed
557
    /// Calculates doublewell of this DOFVector
558
    double DoubleWell(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
559
 
Thomas Witkowski's avatar
Thomas Witkowski committed
560
    /// Calculates the sum of this DOFVector
561
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
562
 
Thomas Witkowski's avatar
Thomas Witkowski committed
563
    /// Sets \ref vec[i] = val, i=0 , ... , size
564
565
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
566
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
567
568
    inline DOFVector<T>& operator=(T d) 
    {
569
570
      set(d); 
      return *this;
571
    } 
572

Thomas Witkowski's avatar
Thomas Witkowski committed
573
    /// Assignment operator between two vectors
574
    DOFVector<T>& operator=(const DOFVector<T>&); 
575

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

579
    /// Returns minimum of DOFVector.
580
581
    T min() const; 

582
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
583
584
    T max() const;

585
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
586
    T absMax() const;
587

588
589
590
591
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
592
593
    void print() const; 

594
    ///
595
596
    int calcMemoryUsage() const;

Praetorius, Simon's avatar
Praetorius, Simon committed
597
598
599
600
    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
601
602
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

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

605

Praetorius, Simon's avatar
Praetorius, Simon committed
606
607
608
609
610
611
    /** 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 
612
613
614
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
Praetorius, Simon's avatar
Praetorius, Simon committed
615
      return *value;
616
    }
617

Praetorius, Simon's avatar
Praetorius, Simon committed
618
619
620
621
    /** 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;
622

623
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
624
625
    void serialize(std::ostream &out) 
    {
626
627
628
      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));
629
    }
630

631
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
    void deserialize(std::istream &in) 
    {
634
635
636
637
      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));
638
    }
639

Praetorius, Simon's avatar
Praetorius, Simon committed
640
//     DOFVector<WorldVector<T> > *getGradient(DOFVector<WorldVector<T> >*) const;
641
642
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
643
644
645

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

Praetorius, Simon's avatar
Praetorius, Simon committed
646
//     DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
647
648
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
649
650

  protected: 
651

652
    /// Data container
653
    std::vector<T> vec; 
654
 
655
    /// Specifies what operation should be performed after coarsening
656
657
658
    CoarsenOperation coarsenOperation;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
659

660
  template<>
661
662
663
664
665
666
667
668
  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
669
670
  const double& DOFVector<double>::evalAtPoint(
    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
671
672

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

676
  template<>
677
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
678
679

  template<>
680
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
681

Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
684
    return v.min();
685
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
686

Thomas Witkowski's avatar
Thomas Witkowski committed
687
688
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
689
    return v.max();
690
  }
691
692
693
694
695
696
697
698
699
700


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Praetorius, Simon's avatar
Praetorius, Simon committed
701
702
703
704
    /** \brief
     * Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
     * as DOFContainer at DOFAdmin
     */
705
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
706
707
708
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
709
    }
710
  
711
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
712
713
    ~DOFVectorDOF() 
    {
714
      feSpace->getAdmin()->removeDOFContainer(this);
715
    }
716

Praetorius, Simon's avatar
Praetorius, Simon committed
717
718
719
720
    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVector<DegreeOfFreedom>::operator[]()
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
721
722
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
723
      return DOFVector<DegreeOfFreedom>::operator[](i);
724
    }
725

Thomas Witkowski's avatar
Thomas Witkowski committed
726
727
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
728
729
730
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

731
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
732
733
    int getSize() const 
    {
734
      return DOFVector<DegreeOfFreedom>::getSize();
735
    }
736

737
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
738
739
    void resize(int size) 
    {
740
741
742
743
744
745
746
747
748
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
749

750
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
751
752
  double norm(DOFVector<T> *vec) 
  {
753
    return vec->nrm2();
754
  }
755
756

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
757
758
  double L2Norm(DOFVector<T> *vec) 
  {
759
    return vec->L2Norm();
760
  }
761
762

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
763
764
  double H1Norm(DOFVector<T> *vec) 
  {
765
    return vec->H1Norm();
766
  }
767
768

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
769
  void print(DOFVector<T> *vec) 
770
  {
771
    vec->print();
772
  }
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804

  // 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>
805
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
806
807
808
809
810
811
812
813
814
815
816
817
818

  // 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
819
820
  inline void scal(T a, DOFVector<T>& y) 
  {
821
822
    y *= a;
  }
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838

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

839
840
841
842
843
844
845
846
847
  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);


848
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
849
850
  inline void set(DOFVector<T>& vec, T d) 
  {
851
    vec.set(d);
852
  }
853
854

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
855
856
  inline void setValue(DOFVector<T>& vec, T d) 
  {
857
    vec.set(d);
858
  }
859
860

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
861
862
  inline int size(DOFVector<T> *vec) 
  {
863
    return vec->getUsedSize();
864
  }
865

866
867
868
869
870
871
  template<typename T>
  inline int size(const DOFVector<T>& vec) 
  {
    return vec.getUsedSize();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
872
873
874
875
876
877
878
879
880
881
  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());
  }

882
883
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
884
885
886
887
  
  template<typename T>
  std::vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     std::vector<DOFVector<double>*> *res);
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912

  
  /// 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);
  }
913
  
914
  /** \brief
915
916
917
   * 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
918
   * meshes.
919
   */
920
921
922
923
924
925
926
927
928
929
930
931
932
  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
933
  /// Computes the integral of a function with two DOFVectors defined on the same mesh.
934
935
936
937
  template<typename TOut, typename T1, typename T2>
  TOut int_Vec2_SingleMesh(const DOFVector<T1> &vec1,
		       const DOFVector<T2> &vec2,
		       BinaryAbstractFunction<TOut, T1, T2> *fct);
938

Praetorius, Simon's avatar
Praetorius, Simon committed
939
  /// Computes the integral of a function with two DOFVectors defined on different meshes.
940
941
942
943
944
945
946
947
948
949
950
951
952
953
  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);
954
  
955
956
957
958
959
960
961
962
  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$
963
964
965
  double integrateGeneral(const std::vector<DOFVector<double>*> &vecs,
			  AbstractFunction<double, std::vector<double> > *fct);
  
966
  /// Computes the integral: \f$ \int f(\{v_i(\vec{x})\}_i,\{\nabla w_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$
967
968
  double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
				  const std::vector<DOFVector<double>*> &grds,
Praetorius, Simon's avatar
Praetorius, Simon committed
969
				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > > *fct);
970

971
972
<