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

48
49
50
51
52
53
namespace AMDiS {

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
84
85
    /// Evaluates the gradient of a DOF vector at a set of quadrature points
    /// defined on the given element.
86
87
88
89
90
91
92
93
94
95
96
    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
97
98
    /// Evaluates the comp'th component of the derivative of a DOF vector at a
    /// set of quadrature points defined on the given element.
99
100
101
102
103
104
105
106
107
108
109
110
111
    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
112
113
    /// Evaluates the jacobian of a DOF vector at a set of quadrature points
    ///  defined on the given element.
114
115
116
117
    void getD2AtQPs(const ElInfo *elInfo,
		    const Quadrature *quad,
		    const FastQuadrature *quadFast,
		    mtl::dense_vector<typename D2Type<T>::type> &d2AtQPs) const;
118

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

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

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

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

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

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

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

167
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
168
    {
169
      return operatorEstFactor.begin();
170
    }
171

172
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
173
    {
174
      return operatorEstFactor.end();
175
    }
176
177


178
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
179
    {
180
      return operators.begin();
181
    }
182

183
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
184
    {
185
      return operators.end();
186
    }
187
188
189

    Flag getAssembleFlag();

190
191
192
193
    /// 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
194
195
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

196
    inline std::vector<Operator*>& getOperators() 
Thomas Witkowski's avatar
Thomas Witkowski committed
197
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
198
      return operators; 
199
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
200

201
    inline std::vector<double*>& getOperatorFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
202
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
203
      return operatorFactor; 
204
    }
205

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

Thomas Witkowski's avatar
Thomas Witkowski committed
211
    /// Returns \ref name
212
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
213
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
214
      return name; 
215
    } 
216

217
    inline void setName(std::string n)
218
    {
219
      name = n;
220
221
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
222
223
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
224
      return boundaryManager; 
225
    }
226

Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
229
      boundaryManager = bm;
230
    }
231

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
232
233
234
235
236
237
    inline void setDirichletDofValue(DegreeOfFreedom dof,
				     T value)
    {
      dirichletDofValues[dof] = value;
    }

238
    std::map<DegreeOfFreedom, T>& getDirichletValues()
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
239
240
241
242
    {
      return dirichletDofValues;
    }

243
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
244
    ///
245
246
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
247
    ///
248
    std::string name;
249

Thomas Witkowski's avatar
Thomas Witkowski committed
250
    ///
251
    ElementVector elementVector;
252

Thomas Witkowski's avatar
Thomas Witkowski committed
253
    ///
254
    std::vector<Operator*> operators;
255

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    ///
263
264
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
265
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
266
    int nBasFcts;
267

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
270

271
    std::map<DegreeOfFreedom, T> dirichletDofValues;
Thomas Witkowski's avatar
Thomas Witkowski committed
272
  };
273

274

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

  /** \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)
308
      {}
309
310
311

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
312
      {}
313
314
315
316
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
317
318
319
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
320

Thomas Witkowski's avatar
Thomas Witkowski committed
321
322
      DOFVector<T> *create() 
      { 
323
	return new DOFVector<T>(feSpace, ""); 
324
      }
325

Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
      void free(DOFVector<T> *vec) 
      { 
328
	delete vec; 
329
      }
330
331
332
333
334
335

    private:
      FiniteElemSpace *feSpace;
    };

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    /// Empty constructor. No initialization!
337
338
339
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
340
    {}
341

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

Thomas Witkowski's avatar
Thomas Witkowski committed
345
    /// Initialization.
346
    void init(const FiniteElemSpace* f, std::string n);
347

Thomas Witkowski's avatar
Thomas Witkowski committed
348
    /// Copy Constructor
349
    DOFVector(const DOFVector& rhs) : DOFVectorBase<T>()
Thomas Witkowski's avatar
Thomas Witkowski committed
350
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
351
      *this = rhs;   
352
353
354
      this->name = rhs.name + "copy";
      if (this->feSpace && this->feSpace->getAdmin())
	(dynamic_cast<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
355
    }
356

Thomas Witkowski's avatar
Thomas Witkowski committed
357
    /// Destructor
358
359
    virtual ~DOFVector();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
366
    /// Returns iterator to the end of \ref vec
367
    typename std::vector<T>::iterator end() 
Thomas Witkowski's avatar
Thomas Witkowski committed
368
    { 
369
      return vec.end(); 
370
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
371
372

    /// Returns const_iterator to the begin of \ref vec
373
    typename std::vector<T>::const_iterator begin() const
Praetorius, Simon's avatar
Praetorius, Simon committed
374
375
376
377
378
    {
      return vec.begin();
    }

    /// Returns const_iterator to the end of \ref vec
379
    typename std::vector<T>::const_iterator end() const
Praetorius, Simon's avatar
Praetorius, Simon committed
380
381
382
383
    {
      return vec.end();
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
384
385
    /// Used by DOFAdmin to compress this DOFVector. Implementation of
    /// \ref DOFIndexedBase::compress()
386
    virtual void compressDOFIndexed(int first, int last,
387
				    std::vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
388
389

    /// Sets \ref coarsenOperation to op
390
    inline void setCoarsenOperation(RefineCoarsenOperation op) 
Thomas Witkowski's avatar
Thomas Witkowski committed
391
    { 
392
      coarsenOperation = op; 
393
    }
394

Thomas Witkowski's avatar
Thomas Witkowski committed
395
    /// Returns \ref coarsenOperation
396
    inline RefineCoarsenOperation getCoarsenOperation() 
Thomas Witkowski's avatar
Thomas Witkowski committed
397
    { 
398
      return coarsenOperation; 
399
    }
400

401
402
403
404
405
406
407
408
409
410
411
412
    /// 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
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
420
    std::vector<T>& getVector() 
Thomas Witkowski's avatar
Thomas Witkowski committed
421
    { 
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_DBG("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_DBG("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_DBG("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_DBG("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
    double Int(Quadrature* q = nullptr) const
Thomas Witkowski's avatar
Thomas Witkowski committed
475
476
477
478
479
480
481
482
483
484
485
486
487
488
    {
      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.
     */
489
    double Int(int meshLevel, Quadrature* q = nullptr) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
490

491

Thomas Witkowski's avatar
Thomas Witkowski committed
492
493
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
494
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = nullptr) const
495
496
497
498
499
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
500

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

Thomas Witkowski's avatar
Thomas Witkowski committed
511
    /// Calculates L1 norm of this DOFVector
512
    double L1Norm(Quadrature* q = nullptr) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
513
 
Thomas Witkowski's avatar
Thomas Witkowski committed
514
    /// Calculates L2 norm of this DOFVector
515
    inline double L2Norm(Quadrature* q = nullptr) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
516
    {
517
      return sqrt(L2NormSquare());
518
    }
519

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

Thomas Witkowski's avatar
Thomas Witkowski committed
523
    /// Calculates H1 norm of this DOFVector
524
    inline double H1Norm(Quadrature* q = nullptr) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
525
    {
526
      return sqrt(H1NormSquare());
527
    }
528

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

Thomas Witkowski's avatar
Thomas Witkowski committed
532
    /// Calculates euclidian norm of this DOFVector
533
534
    double nrm2() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
535
    /// Returns square of the euclidian norm.
536
537
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
538
    /// Calculates l2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
539
540
    inline double l2norm() const 
    { 
541
      return nrm2();
542
    }
543

Thomas Witkowski's avatar
Thomas Witkowski committed
544
    /// Calculates the absolute sum of this DOFVector
545
546
    T asum() const; 

Thomas Witkowski's avatar
Thomas Witkowski committed
547
    /// Calculates the l1 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
548
549
    inline double l1norm() const 
    { 
550
      return asum();
551
    } 
552

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
569
    /// Assignment operator between two vectors
570
    DOFVector<T>& operator=(const DOFVector<T>&); 
571

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

575
    /// Returns minimum of DOFVector.
576
577
    T min() const; 

578
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
579
580
    T max() const;

581
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
582
    T absMax() const;
583

584
585
586
587
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
588
589
    void print() const; 

590
    ///
591
    size_t calcMemoryUsage() const;
592

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

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

599

600
    /// Eval DOFVector at given point p. If oldElInfo != nullptr the search for 
Thomas Witkowski's avatar
Thomas Witkowski committed
601
602
    /// the element, where p is inside, starts from oldElInfo. implemented for:
    /// double, WorldVector< double >
603
    T evalAtPoint(WorldVector<double> &p, 
604
		  ElInfo *oldElInfo = nullptr) const 
605
606
607
608
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
    }
609

Thomas Witkowski's avatar
Thomas Witkowski committed
610
611
612
    /// Determine the DegreeOfFreedom that has coords with minimal euclidean 
    /// distance to WorldVector p. return true if DOF is found, and false 
    /// otherwise.
613
    bool getDofIdxAtPoint(WorldVector<double> &p, 
Thomas Witkowski's avatar
Thomas Witkowski committed
614
				DegreeOfFreedom &idx, 
615
				ElInfo *oldElInfo = nullptr, 
Thomas Witkowski's avatar
Thomas Witkowski committed
616
				bool useOldElInfo = false) const;
617

618
    /// Writes the data of the DOFVector to an output stream.
619
    void serialize(std::ostream &out) 
Thomas Witkowski's avatar
Thomas Witkowski committed
620
    {
621
      unsigned int size = static_cast<unsigned int>(vec.size());
622
623
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
624
    }
625

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

635
636
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
637
638
639

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

640
641
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
642
643

  protected: 
644

645
    /// Data container
646
    std::vector<T> vec; 
647
 
648
    /// Specifies what operation should be performed after coarsening
649
650
    RefineCoarsenOperation coarsenOperation;
    RefineCoarsenOperation refineOperation;
651
652
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
653

654
  template<>
655
656
657
658
659
660
661
662
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

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

  template<>
663
664
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
					 ElInfo *oldElInfo) const;
665
666

  template<>
667
668
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const;
669

670
  template<>
671
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
672
673

  template<>
674
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
675

Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
  inline double min(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
678
    return v.min();
679
  } 
Thomas Witkowski's avatar
Thomas Witkowski committed
680

Thomas Witkowski's avatar
Thomas Witkowski committed
681
682
  inline double max(const DOFVector<double>& v) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
683
    return v.max();
684
  }
685
686
687
688
689
690
691
692
693
694


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

Thomas Witkowski's avatar
Thomas Witkowski committed
709
710
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
711
712
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
713
      return DOFVector<DegreeOfFreedom>::operator[](i);
714
    }
715

Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
718
719
720
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

721
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
722
723
    int getSize() const 
    {
724
      return DOFVector<DegreeOfFreedom>::getSize();
725
    }
726

727
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
728
729
    void resize(int size) 
    {
730
731
732
733
734
735
736
737
738
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
739

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

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

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
753
754
  double H1Norm(DOFVector<T> *vec) 
  {
755
    return vec->H1Norm();
756
  }
757
758

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
759
  void print(DOFVector<T> *vec) 
760
  {
761
    vec->print();
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
789
790
791
792
793
794

  // 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>
795
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
796
797
798
799
800
801
802
803
804
805
806
807
808

  // 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
809
810
  inline void scal(T a, DOFVector<T>& y) 
  {
811
812
    y *= a;
  }
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828

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

829
830
831
832
833
834
835
836
837
  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);


838
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
839
840
  inline void set(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 void setValue(DOFVector<T>& vec, T d) 
  {
847
    vec.set(d);
848
  }
849
850

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
851
852
  inline int size(DOFVector<T> *vec) 
  {
853
    return vec->getUsedSize();
854
  }
855

856
857
858
859
860
861
  template<typename T>
  inline int size(const DOFVector<T>& vec) 
  {
    return vec.getUsedSize();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
862
  template<typename T>
863
  inline void checkFeSpace(const FiniteElemSpace* feSpace, const std::vector<T>& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
864
  {
865
    FUNCNAME_DBG("checkFeSpace()");
866
867
    TEST_EXIT_DBG(feSpace)("feSpace is nullptr\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is nullptr\n");
Thomas Witkowski's avatar
Thomas Witkowski committed