DOFVector.h 27.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
23
24
25
26
27

/** \file DOFVector.h */

#ifndef AMDIS_DOFVECTOR_H 
#define AMDIS_DOFVECTOR_H 
 
28
29
#include <vector> 
#include <memory> 
30
#include <map> 
31

32
#include "AMDiS_fwd.h"
33
34
35
36
37
38
39
40
41
42
43
#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
44
#include "BasisFunction.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
45
#include "FiniteElemSpace.h"
46
#include "SurfaceQuadrature.h"
47
#include "Traits.h"
48

49
50
51
52
53
54
namespace AMDiS {

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

56
    DOFVectorBase() 
57
      : feSpace(nullptr),
58
	elementVector(3),
59
        boundaryManager(nullptr),
Thomas Witkowski's avatar
Thomas Witkowski committed
60
        nBasFcts(0)
61
    {}
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    
63
    DOFVectorBase(const FiniteElemSpace *f, std::string n);
64

Thomas Witkowski's avatar
Thomas Witkowski committed
65
    virtual ~DOFVectorBase();
66

Thomas Witkowski's avatar
Thomas Witkowski committed
67
68
    /// For the given element, this function returns an array of all DOFs of 
    /// this DOFVector that are defined on this element.
69
70
    virtual void getLocalVector(const Element *el, 
				mtl::dense_vector<T>& localVec) const;
71

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

79
80
81
82
83
    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
84

Thomas Witkowski's avatar
Thomas Witkowski committed
85
86
    /// Evaluates the gradient of a DOF vector at a set of quadrature points
    /// defined on the given element.
87
88
89
90
91
92
93
94
95
96
97
    void getGrdAtQPs( const ElInfo *elInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;

    void getGrdAtQPs( const ElInfo *smallElInfo,
		      const ElInfo *largeElInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
98
99
    /// Evaluates the comp'th component of the derivative of a DOF vector at a
    /// set of quadrature points defined on the given element.
100
101
102
103
104
105
106
107
108
109
110
111
112
    void getDerivativeAtQPs( const ElInfo *elInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      int comp,
		      mtl::dense_vector<T> &derivativeAtQPs) const;

    void getDerivativeAtQPs( const ElInfo *smallElInfo,
		      const ElInfo *largeElInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      int comp,
		      mtl::dense_vector<T> &derivativeAtQPs) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
113
114
    /// Evaluates the jacobian of a DOF vector at a set of quadrature points
    ///  defined on the given element.
115
116
117
118
    void getD2AtQPs(const ElInfo *elInfo,
		    const Quadrature *quad,
		    const FastQuadrature *quadFast,
		    mtl::dense_vector<typename D2Type<T>::type> &d2AtQPs) const;
119

Thomas Witkowski's avatar
Thomas Witkowski committed
120
    /// Returns the FE space the DOF vector is defined on.
121
    inline const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
122
    {
123
      return feSpace;
124
    }
125

Thomas Witkowski's avatar
Thomas Witkowski committed
126
127
    /// 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
128
129
    void assemble(T factor, ElInfo *elInfo,			    
		  const BoundaryType *bound, 
130
		  Operator *op = nullptr);
Thomas Witkowski's avatar
Thomas Witkowski committed
131
132
133
134
135

    void assemble2(T factor, 
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
		   const BoundaryType *bound, 
136
		   Operator *op = nullptr);
137
138

    void addElementVector(T sign,
139
			  const ElementVector& elVec, 
140
			  const BoundaryType *bound,
141
			  ElInfo *elInfo,
142
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
143

144
145
146
    ///
    void assembleOperator(Operator &op);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
147
148
149
    /// 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
150
    void finishAssembling();
151
152
 
    inline void addOperator(Operator* op, 
153
154
			    double *factor = nullptr,
			    double *estFactor = nullptr) 
155
156
157
158
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
159
    }
160

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

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

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

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


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

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

    Flag getAssembleFlag();

194
195
196
197
    /// 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
198
199
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
236
237
238
239
240
241
    inline void setDirichletDofValue(DegreeOfFreedom dof,
				     T value)
    {
      dirichletDofValues[dof] = value;
    }

242
    std::map<DegreeOfFreedom, T>& getDirichletValues()
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
243
244
245
246
    {
      return dirichletDofValues;
    }

247
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
248
    ///
249
250
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
251
    ///
252
    std::string name;
253

Thomas Witkowski's avatar
Thomas Witkowski committed
254
    ///
255
    ElementVector elementVector;
256

Thomas Witkowski's avatar
Thomas Witkowski committed
257
    ///
258
    std::vector<Operator*> operators;
259

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
266
    ///
267
268
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
269
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
270
    int nBasFcts;
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
273
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
274

275
    std::map<DegreeOfFreedom, T> dirichletDofValues;
Thomas Witkowski's avatar
Thomas Witkowski committed
276
  };
277

278

Thomas Witkowski's avatar
Thomas Witkowski committed
279
  /// Specifies which operation should be done after coarsening
280
281
282
  typedef enum{
    NO_OPERATION = 0,   
    COARSE_RESTRICT = 1,
283
284
285
    COARSE_INTERPOL = 2, 
    REFINE_INTERPOL = 4
  } RefineCoarsenOperation;
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
 

  /** \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
  {  
303
304
305
306
307
308
  public:
    typedef typename DOFVectorBase<T>::value_type	value_type;
    typedef typename DOFVectorBase<T>::size_type		size_type;
    typedef typename DOFVectorBase<T>::reference		reference;
    typedef typename DOFVectorBase<T>::const_reference	const_reference;
  
309
310
311
312
313
314
315
316
317
  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)
318
      {}
319
320
321

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

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

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

Thomas Witkowski's avatar
Thomas Witkowski committed
352
    /// Constructs a DOFVector with name n belonging to FiniteElemSpace f
353
    DOFVector(const FiniteElemSpace* f, std::string n, bool addToSynch = true); 
354

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    /// Initialization.
356
    void init(const FiniteElemSpace* f, std::string n, bool addToSynch = true);
357

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

Thomas Witkowski's avatar
Thomas Witkowski committed
367
    /// Destructor
368
369
    virtual ~DOFVector();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
376
    /// Returns iterator to the end of \ref vec
377
    typename std::vector<T>::iterator end() 
Thomas Witkowski's avatar
Thomas Witkowski committed
378
    { 
379
      return vec.end(); 
380
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
381
382

    /// Returns const_iterator to the begin of \ref vec
383
    typename std::vector<T>::const_iterator begin() const
Praetorius, Simon's avatar
Praetorius, Simon committed
384
385
386
387
388
    {
      return vec.begin();
    }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
405
    /// Returns \ref coarsenOperation
406
    inline RefineCoarsenOperation getCoarsenOperation() 
Thomas Witkowski's avatar
Thomas Witkowski committed
407
    { 
408
      return coarsenOperation; 
409
    }
410

411
412
413
414
415
416
417
418
419
420
421
422
    /// Sets \ref refineOperation to op
    inline void setRefineOperation(RefineCoarsenOperation op) 
    { 
      refineOperation = op; 
    }

    /// Returns \ref refineOperation
    inline RefineCoarsenOperation getRefineOperation() 
    { 
      return refineOperation; 
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
423
    /// Restriction after coarsening. Implemented for DOFVector<double>
424
    inline void coarseRestrict(RCNeighbourList&, int) {}
425

426
    /// Interpolation after refinement. Implemented for DOFVector<double>
427
    inline void refineInterpol(RCNeighbourList&, int) {}
428

Thomas Witkowski's avatar
Thomas Witkowski committed
429
    /// Returns \ref vec
430
    std::vector<T>& getVector() 
Thomas Witkowski's avatar
Thomas Witkowski committed
431
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
432
      return vec;
433
    }
434

Thomas Witkowski's avatar
Thomas Witkowski committed
435
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
    inline int getSize() const 
    { 
438
      return vec.size();
439
    } 
440

Thomas Witkowski's avatar
Thomas Witkowski committed
441
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
442
443
    inline int getUsedSize() const 
    { 
444
      return this->feSpace->getAdmin()->getUsedSize(); 
445
    }
446

Thomas Witkowski's avatar
Thomas Witkowski committed
447
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
448
449
    inline void resize(int n) 
    { 
450
      FUNCNAME_DBG("DOFVector<T>::resize()");
451
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
452
      vec.resize(n);
453
    } 
454

Thomas Witkowski's avatar
Thomas Witkowski committed
455
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
    inline void resize(int n, T init) 
    { 
458
      FUNCNAME_DBG("DOFVector<T>::resize()");
459
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
460
      vec.resize(n, init);
461
    } 
462

Thomas Witkowski's avatar
Thomas Witkowski committed
463
    /// Returns \ref vec[i]
464
    inline const_reference operator[](DegreeOfFreedom i) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
465
    {
466
      FUNCNAME_DBG("DOFVector<T>::operator[]");
467
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
468
	("Illegal vector index %d.\n", i);
469
      return vec[i];
470
    } 
471

Thomas Witkowski's avatar
Thomas Witkowski committed
472
    /// Returns \ref vec[i]
473
    inline reference operator[](DegreeOfFreedom i) 
Thomas Witkowski's avatar
Thomas Witkowski committed
474
    {
475
      FUNCNAME_DBG("DOFVector<T>::operator[]");
476
477

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

480
      return vec[i];
481
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
482
 
Thomas Witkowski's avatar
Thomas Witkowski committed
483
    /// Calculates Integral of this DOFVector
484
    double Int(Quadrature* q = nullptr) const
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
487
488
489
490
491
492
493
494
495
496
497
498
    {
      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.
     */
499
    double Int(int meshLevel, Quadrature* q = nullptr) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
500

501

Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
504
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = nullptr) const
505
506
507
508
509
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
510

Thomas Witkowski's avatar
Thomas Witkowski committed
511
512
513
    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
514
    double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = nullptr) const
515
516
517
518
519
520
    {
      FUNCNAME("DOFVector::IntOnBoundaryNormal())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
521
    /// Calculates L1 norm of this DOFVector
522
    double L1Norm(Quadrature* q = nullptr) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
523
 
Thomas Witkowski's avatar
Thomas Witkowski committed
524
    /// Calculates L2 norm of this DOFVector
525
    inline double L2Norm(Quadrature* q = nullptr) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
526
    {
527
      return sqrt(L2NormSquare());
528
    }
529

Thomas Witkowski's avatar
Thomas Witkowski committed
530
    /// Calculates square of L2 norm of this DOFVector
531
    double L2NormSquare(Quadrature* q = nullptr) const;
532

Thomas Witkowski's avatar
Thomas Witkowski committed
533
    /// Calculates H1 norm of this DOFVector
534
    inline double H1Norm(Quadrature* q = nullptr) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
535
    {
536
      return sqrt(H1NormSquare());
537
    }
538

Thomas Witkowski's avatar
Thomas Witkowski committed
539
    /// Calculates square of H1 norm of this DOFVector
540
    double H1NormSquare(Quadrature* q = nullptr) const;  
541

Thomas Witkowski's avatar
Thomas Witkowski committed
542
    /// Calculates euclidian norm of this DOFVector
543
544
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
545
    /// Returns square of the euclidian norm.
546
547
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
548
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
549
550
    inline double l2norm() const 
    { 
551
      return nrm2();
552
    }
553

Thomas Witkowski's avatar
Thomas Witkowski committed
554
    /// Calculates the absolute sum of this DOFVector
555
556
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
557
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
558
559
    inline double l1norm() const 
    { 
560
      return asum();
561
    } 
562

Thomas Witkowski's avatar
Thomas Witkowski committed
563
    /// Calculates doublewell of this DOFVector
564
    double DoubleWell(Quadrature* q = nullptr) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
565
 
Thomas Witkowski's avatar
Thomas Witkowski committed
566
    /// Calculates the sum of this DOFVector
567
    T sum() const; 
Thomas Witkowski's avatar
Thomas Witkowski committed
568
 
Thomas Witkowski's avatar
Thomas Witkowski committed
569
    /// Sets \ref vec[i] = val, i=0 , ... , size
570
571
    void set(T val); 

Thomas Witkowski's avatar
Thomas Witkowski committed
572
    /// Assignment operator for setting complete vector to a certain value d
Thomas Witkowski's avatar
Thomas Witkowski committed
573
574
    inline DOFVector<T>& operator=(T d) 
    {
575
576
      set(d); 
      return *this;
577
    } 
578

Thomas Witkowski's avatar
Thomas Witkowski committed
579
    /// Assignment operator between two vectors
580
    DOFVector<T>& operator=(const DOFVector<T>&); 
581

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

585
    /// Returns minimum of DOFVector.
586
587
    T min() const; 

588
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
589
590
    T max() const;

591
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
592
    T absMax() const;
593

594
595
596
597
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
598
599
    void print() const; 

600
    ///
601
    size_t calcMemoryUsage() const;
602

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

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

609

610
    /// Eval DOFVector at given point p. If oldElInfo != nullptr the search for 
Thomas Witkowski's avatar
Thomas Witkowski committed
611
612
    /// the element, where p is inside, starts from oldElInfo. implemented for:
    /// double, WorldVector< double >
613
    T evalAtPoint(WorldVector<double> &p, 
614
		  ElInfo *oldElInfo = nullptr) const 
615
616
617
618
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
    }
619

Thomas Witkowski's avatar
Thomas Witkowski committed
620
621
622
    /// Determine the DegreeOfFreedom that has coords with minimal euclidean 
    /// distance to WorldVector p. return true if DOF is found, and false 
    /// otherwise.
623
    bool getDofIdxAtPoint(WorldVector<double> &p, 
Thomas Witkowski's avatar
Thomas Witkowski committed
624
				DegreeOfFreedom &idx, 
625
				ElInfo *oldElInfo = nullptr, 
Thomas Witkowski's avatar
Thomas Witkowski committed
626
				bool useOldElInfo = false) const;
627

628
629
    template<typename S>
    void serialize(std::ostream &out, boost::mpl::identity<S>) 
Thomas Witkowski's avatar
Thomas Witkowski committed
630
    {
631
      unsigned int size = static_cast<unsigned int>(vec.size());
632
633
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
634
    }
635
636
637
638
639
640
641
642
643
644
645
646
    
    // no serialization for bool vectors
    void serialize(std::ostream &out, boost::mpl::identity<bool>) {}
    
    /// Writes the data of the DOFVector to an output stream.
    void serialize(std::ostream &out) 
    {
      serialize(out, boost::mpl::identity<T>());
    }
    
    template<typename S>
    void deserialize(std::istream &in, boost::mpl::identity<S>) 
Thomas Witkowski's avatar
Thomas Witkowski committed
647
    {
648
649
650
651
      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));
652
    }
653
654
655
656
657
658
659
660
661
    
    // no de-serialization for bool vectors
    void deserialize(std::istream &in, boost::mpl::identity<bool>) {}

    /// Reads data of an DOFVector from an input stream.
    void deserialize(std::istream &in) 
    {
      deserialize(in, boost::mpl::identity<T>());
    }
662

663
664
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
665
666
667

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

668
669
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
670
671

  protected: 
672

673
    /// Data container
674
    std::vector<T> vec; 
675
 
676
    /// Specifies what operation should be performed after coarsening
677
678
    RefineCoarsenOperation coarsenOperation;
    RefineCoarsenOperation refineOperation;
679
680
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
681

682
  template<>
683
684
685
686
687
688
689
690
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

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

  template<>
691
692
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
					 ElInfo *oldElInfo) const;
693
694

  template<>
695
696
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const;
697

698
  template<>
699
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
700
701

  template<>
702
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
703

Thomas Witkowski's avatar
Thomas Witkowski committed
704
705
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
706
    return v.min();
707
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
708

Thomas Witkowski's avatar
Thomas Witkowski committed
709
710
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
711
    return v.max();
712
  }
713
714
715
716
717
718
719
720
721
722


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
Thomas Witkowski's avatar
Thomas Witkowski committed
723
724
    /// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
    /// as DOFContainer at DOFAdmin
725
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
726
727
728
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
729
    }
730
  
731
    /// Deregisters itself at DOFAdmin.
Thomas Witkowski's avatar
Thomas Witkowski committed
732
733
    ~DOFVectorDOF() 
    {
734
      feSpace->getAdmin()->removeDOFContainer(this);
735
    }
736

Thomas Witkowski's avatar
Thomas Witkowski committed
737
738
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
741
      return DOFVector<DegreeOfFreedom>::operator[](i);
742
    }
743

Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
746
747
748
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

749
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
    int getSize() const 
    {
752
      return DOFVector<DegreeOfFreedom>::getSize();
753
    }
754

755
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
756
757
    void resize(int size) 
    {
758
759
760
761
762
763
764
765
766
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
767

768
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
769
770
  double norm(DOFVector<T> *vec) 
  {
771
    return vec->nrm2();
772
  }
773
774

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
  double L2Norm(DOFVector<T> *vec) 
  {
777
    return vec->L2Norm();
778
  }
779
780

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
781
782
  double H1Norm(DOFVector<T> *vec) 
  {
783
    return vec->H1Norm();
784
  }
785
786

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
787
  void print(DOFVector<T> *vec) 
788
  {
789
    vec->print();
790
  }
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822

  // 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>
823
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
824
825
826
827
828
829
830
831
832
833
834
835
836

  // 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
837
838
  inline void scal(T a, DOFVector<T>& y) 
  {
839
840
    y *= a;
  }
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856

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

857
858
859
860
861
862
863
864
865
  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);


866
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
867
868
  inline void set(DOFVector<T>& vec, T d) 
  {
869
    vec.set(d);
870
  }
871
872

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
873
874
  inline void setValue(DOFVector<T>& vec, T d) 
  {
875
    vec.set(d);
876
  }
877

878
879
880
881
882
883
884
885
886
887
888
//   template<typename T>
//   inline int size(DOFVector<T> *vec) 
//   {
//     return vec->getUsedSize();
//   }
// 
//   template<typename T>
//   inline int size(const DOFVector<T>& vec) 
//   {
//     return vec.getUsedSize();
//   }
889