DOFVector.h 26.8 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * 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.
18
 *
19
 ******************************************************************************/
20
21


22
23
24

/** \file DOFVector.h */

25
26
27
28
29
30
#ifndef AMDIS_DOFVECTOR_H
#define AMDIS_DOFVECTOR_H

#include <vector>
#include <memory>
#include <map>
31

32
#include "AMDiS_fwd.h"
33
#include "FixVec.h"
34
35
36
#include "Global.h"
#include "Flag.h"
#include "RCNeighbourList.h"
37
38
39
40
41
42
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "DOFContainer.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
43
#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
namespace AMDiS {

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

56
    DOFVectorBase()
57
      : feSpace(NULL),
58
	elementVector(3),
59
        boundaryManager(NULL),
Thomas Witkowski's avatar
Thomas Witkowski committed
60
        nBasFcts(0)
61
    {}
62

63
    DOFVectorBase(const FiniteElemSpace *f, std::string n);
64

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

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

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

79
    void getVecAtQPs(const ElInfo *smallElInfo,
80
81
82
83
		     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.
128
129
    void assemble(T factor, ElInfo *elInfo,
		  const BoundaryType *bound,
130
		  Operator *op = NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
131

132
    void assemble2(T factor,
Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
135
		   const BoundaryType *bound,
136
		   Operator *op = NULL);
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
    ///
    void assembleOperator(Operator &op);
146

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
			    double *factor = NULL,
154
			    double *estFactor = NULL)
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
    /// Evaluates \f[ u_h(x(\lambda)) = \sum_{i=0}^{m-1} vec[ind[i]] *
195
196
197
    /// \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
201
202
    inline std::vector<Operator*>& getOperators()
    {
      return operators;
203
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
204

205
206
207
    inline std::vector<double*>& getOperatorFactor()
    {
      return operatorFactor;
208
    }
209

210
211
212
    inline std::vector<double*>& getOperatorEstFactor()
    {
      return operatorEstFactor;
213
    }
214

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

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

226
227
228
    inline BoundaryManager* getBoundaryManager() const
    {
      return boundaryManager;
229
    }
230

231
    inline void setBoundaryManager(BoundaryManager *bm)
Thomas Witkowski's avatar
Thomas Witkowski committed
232
    {
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

279

280
  /** \ingroup DOFAdministration
281
   * \brief
282
283
284
285
286
287
   * 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
288
289
290
291
   * 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.
   */
292
  template<typename T>
293
  class DOFVector : public DOFVectorBase<T>, public Serializable
294
  {
295
296
297
298
299
  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;
300

301
302
303
304
305
306
307
308
309
  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)
310
      {}
311
312
313

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

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

323
324
325
      DOFVector<T> *create()
      {
	return new DOFVector<T>(feSpace, "");
326
      }
327

328
329
330
      void free(DOFVector<T> *vec)
      {
	delete vec;
331
      }
332
333
334
335
336
337

    private:
      FiniteElemSpace *feSpace;
    };

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

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

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

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

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

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

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

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

    /// Returns const_iterator to the end of \ref vec
380
    typename std::vector<T>::const_iterator end() const
Praetorius, Simon's avatar
Praetorius, Simon committed
381
382
383
    {
      return vec.end();
    }
384

Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
    /// Used by DOFAdmin to compress this DOFVector. Implementation of
    /// \ref DOFIndexedBase::compress()
387
    virtual void compressDOFIndexed(int first, int last,
388
				    std::vector<DegreeOfFreedom> &newDof);
389

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    /// Restriction after coarsening. Implemented for DOFVector<double>
391
    void coarseRestrict(RCNeighbourList&, int) {}
392

393
    /// Interpolation after refinement. Implemented for DOFVector<double>
394
395
396
397
    void refineInterpol(RCNeighbourList& r, int i)
    {
        refineInterpolImpl(r, i, id<T>());
    }
398

Thomas Witkowski's avatar
Thomas Witkowski committed
399
    /// Returns \ref vec
400
401
    std::vector<T>& getVector()
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
402
      return vec;
403
    }
404

Thomas Witkowski's avatar
Thomas Witkowski committed
405
    /// Returns size of \ref vec
406
407
    inline int getSize() const
    {
408
      return vec.size();
409
    }
410

Thomas Witkowski's avatar
Thomas Witkowski committed
411
    /// Returns used size of the vector.
412
413
414
    inline int getUsedSize() const
    {
      return this->feSpace->getAdmin()->getUsedSize();
415
    }
416

Thomas Witkowski's avatar
Thomas Witkowski committed
417
    /// Resizes \ref vec to n
418
419
    inline void resize(int n)
    {
420
      FUNCNAME_DBG("DOFVector<T>::resize()");
421
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n");
422
      vec.resize(n);
423
    }
424

Thomas Witkowski's avatar
Thomas Witkowski committed
425
    /// Resizes \ref vec to n and inits new values with init
426
427
    inline void resize(int n, T init)
    {
428
      FUNCNAME_DBG("DOFVector<T>::resize()");
429
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n");
430
      vec.resize(n, init);
431
    }
432

Thomas Witkowski's avatar
Thomas Witkowski committed
433
    /// Returns \ref vec[i]
434
    inline const_reference operator[](DegreeOfFreedom i) const
Thomas Witkowski's avatar
Thomas Witkowski committed
435
    {
436
      FUNCNAME_DBG("DOFVector<T>::operator[]");
437
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
438
	("Illegal vector index %d.\n", i);
439
      return vec[i];
440
    }
441

Thomas Witkowski's avatar
Thomas Witkowski committed
442
    /// Returns \ref vec[i]
443
    inline reference operator[](DegreeOfFreedom i)
Thomas Witkowski's avatar
Thomas Witkowski committed
444
    {
445
      FUNCNAME_DBG("DOFVector<T>::operator[]");
446

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

450
      return vec[i];
451
    }
452

Thomas Witkowski's avatar
Thomas Witkowski committed
453
    /// Calculates Integral of this DOFVector
454
    double Int(Quadrature* q = NULL) const
Thomas Witkowski's avatar
Thomas Witkowski committed
455
456
457
458
459
    {
      return Int(-1, q);
    }

    /** \brief
460
     * Calculates Integral of this DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
463
464
465
466
467
468
     *
     * \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.
     */
469
    double Int(int meshLevel, Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
470

471

Thomas Witkowski's avatar
Thomas Witkowski committed
472
473
    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
474
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
475
476
477
478
479
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
480

481
482
    /// Calculates Integral of this DOFVector times normal vector over parts
    /// of the domain boundary, indicated by boundaryType. Implemented for
Thomas Witkowski's avatar
Thomas Witkowski committed
483
    /// DOFVector<WorldVector<double> >
484
    double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = NULL) const
485
486
487
488
489
490
    {
      FUNCNAME("DOFVector::IntOnBoundaryNormal())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
494
    /// Calculates L2 norm of this DOFVector
495
    inline double L2Norm(Quadrature* q = NULL) const
Thomas Witkowski's avatar
Thomas Witkowski committed
496
    {
Wieland Marth's avatar
Wieland Marth committed
497
      return std::sqrt(L2NormSquare());
498
    }
499

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

Thomas Witkowski's avatar
Thomas Witkowski committed
503
    /// Calculates H1 norm of this DOFVector
504
    inline double H1Norm(Quadrature* q = NULL) const
Thomas Witkowski's avatar
Thomas Witkowski committed
505
    {
Wieland Marth's avatar
Wieland Marth committed
506
      return std::sqrt(H1NormSquare());
507
    }
508

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

Thomas Witkowski's avatar
Thomas Witkowski committed
512
    /// Calculates euclidian norm of this DOFVector
513
    double nrm2() const;
514

Thomas Witkowski's avatar
Thomas Witkowski committed
515
    /// Returns square of the euclidian norm.
516
517
    double squareNrm2() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
518
    /// Calculates l2 norm of this DOFVector
519
520
    inline double l2norm() const
    {
521
      return nrm2();
522
    }
523

Thomas Witkowski's avatar
Thomas Witkowski committed
524
    /// Calculates the absolute sum of this DOFVector
525
    T asum() const;
526

Thomas Witkowski's avatar
Thomas Witkowski committed
527
    /// Calculates the l1 norm of this DOFVector
528
529
    inline double l1norm() const
    {
530
      return asum();
531
    }
532

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

Thomas Witkowski's avatar
Thomas Witkowski committed
536
    /// Calculates the sum of this DOFVector
537
538
    T sum() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
539
    /// Sets \ref vec[i] = val, i=0 , ... , size
540
    void set(T val);
541

Thomas Witkowski's avatar
Thomas Witkowski committed
542
    /// Assignment operator for setting complete vector to a certain value d
543
    inline DOFVector<T>& operator=(T d)
Thomas Witkowski's avatar
Thomas Witkowski committed
544
    {
545
      set(d);
546
      return *this;
547
    }
548

Thomas Witkowski's avatar
Thomas Witkowski committed
549
    /// Assignment operator between two vectors
550
    DOFVector<T>& operator=(const DOFVector<T>&);
551

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

555
    /// Returns minimum of DOFVector.
556
    T min() const;
557

558
    /// Returns maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
559
560
    T max() const;

561
    /// Returns absolute maximum of DOFVector.
Thomas Witkowski's avatar
Thomas Witkowski committed
562
    T absMax() const;
563

564
565
566
567
    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
568
    void print() const;
569

570
    ///
571
    size_t calcMemoryUsage() const;
572

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

Praetorius, Simon's avatar
Praetorius, Simon committed
577
    void interpol(const DOFVector<T> *v, double factor = 1.0);
578

579

580
    /// Eval DOFVector at given point p. If oldElInfo != NULL the search for
Thomas Witkowski's avatar
Thomas Witkowski committed
581
582
    /// the element, where p is inside, starts from oldElInfo. implemented for:
    /// double, WorldVector< double >
583
584
    T evalAtPoint(WorldVector<double> const& p,
		  ElInfo *oldElInfo = NULL) const
585
    {
586
        return evalAtPointImpl(p, oldElInfo, id<T>());
587
    }
588

589
590
591
592
    T operator()(WorldVector<double> const& p) const
    {
      return evalAtPoint(p);
    }
593

594
595
    /// Determine the DegreeOfFreedom that has coords with minimal euclidean
    /// distance to WorldVector p. return true if DOF is found, and false
Thomas Witkowski's avatar
Thomas Witkowski committed
596
    /// otherwise.
597
598
599
    bool getDofIdxAtPoint(WorldVector<double> &p,
				DegreeOfFreedom &idx,
				ElInfo *oldElInfo = NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
600
				bool useOldElInfo = false) const;
601

602
    template<typename S>
603
    void serialize(std::ostream &out, boost::mpl::identity<S>)
Thomas Witkowski's avatar
Thomas Witkowski committed
604
    {
605
      unsigned int size = static_cast<unsigned int>(vec.size());
606
607
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
608
    }
609

610
611
    // no serialization for bool vectors
    void serialize(std::ostream &out, boost::mpl::identity<bool>) {}
612

613
    /// Writes the data of the DOFVector to an output stream.
614
    void serialize(std::ostream &out)
615
616
617
    {
      serialize(out, boost::mpl::identity<T>());
    }
618

619
    template<typename S>
620
    void deserialize(std::istream &in, boost::mpl::identity<S>)
Thomas Witkowski's avatar
Thomas Witkowski committed
621
    {
622
623
624
625
      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));
626
    }
627

628
629
630
631
    // no de-serialization for bool vectors
    void deserialize(std::istream &in, boost::mpl::identity<bool>) {}

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

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

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

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

645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
  protected:

    template <class S>
    void refineInterpolImpl(RCNeighbourList& r, int i, id<S>) { }

    // some specializations for concrete types
    void refineInterpolImpl(RCNeighbourList& r, int i, id<double>);
    void refineInterpolImpl(RCNeighbourList& r, int i, id<WorldVector<double> >);

    template <class S>
    S evalAtPointImpl(WorldVector<double> const&, ElInfo*, id<S>) const
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
    }

    // some specializations for concrete types
    double evalAtPointImpl(WorldVector<double> const&, ElInfo*, id<double>) const;
    WorldVector<double> evalAtPointImpl(WorldVector<double> const&, ElInfo*, id<WorldVector<double> >) const;
664

665
    /// Data container
666
667
    std::vector<T> vec;
  };
668

Thomas Witkowski's avatar
Thomas Witkowski committed
669

670
  template<>
671
672
673
674
675
676
677
678
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

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

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

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

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


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

708
    /// Deregisters itself at DOFAdmin.
709
    ~DOFVectorDOF()
Thomas Witkowski's avatar
Thomas Witkowski committed
710
    {
711
      feSpace->getAdmin()->removeDOFContainer(this);
712
    }
713

714
    /// Implements DOFContainer::operator[]() by calling
Thomas Witkowski's avatar
Thomas Witkowski committed
715
    /// DOFVector<DegreeOfFreedom>::operator[]()
716
    DegreeOfFreedom& operator[](DegreeOfFreedom i)
Thomas Witkowski's avatar
Thomas Witkowski committed
717
    {
718
      return DOFVector<DegreeOfFreedom>::operator[](i);
719
    }
720

721
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const
Thomas Witkowski's avatar
Thomas Witkowski committed
722
    {
723
724
725
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

726
    /// Implements DOFIndexedBase::getSize()
727
    int getSize() const
Thomas Witkowski's avatar
Thomas Witkowski committed
728
    {
729
      return DOFVector<DegreeOfFreedom>::getSize();
730
    }
731

732
    /// Implements DOFIndexedBase::resize()
733
    void resize(int size)
Thomas Witkowski's avatar
Thomas Witkowski committed
734
    {
735
736
737
738
739
740
741
742
743
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
744

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

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

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

  template<typename T>
764
  void print(DOFVector<T> *vec)
765
  {
766
    vec->print();
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
795
796
797
798
799

  // 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>
800
  void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
801
802
803

  // matrix vector product
  template<typename T>
804
805
  void mv(MatrixTranspose transpose,
	  const DOFMatrix &a,
806
807
	  const DOFVector<T> &x,
	  DOFVector<T> &result,
808
	  bool add = false);
809
810
811
812
813

  template<typename T>
  void xpay(double a,const DOFVector<T>& x,DOFVector<T>& y);

  template<typename T>
814
  inline void scal(T a, DOFVector<T>& y)
Thomas Witkowski's avatar
Thomas Witkowski committed
815
  {
816
817
    y *= a;
  }
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832