DOFVector.h 26.3 KB
Newer Older
1 2 3 4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6 7
// ==                                                                        ==
// ============================================================================
8 9 10 11 12 13 14 15 16 17 18 19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20 21 22 23 24 25

/** \file DOFVector.h */

#ifndef AMDIS_DOFVECTOR_H 
#define AMDIS_DOFVECTOR_H 
 
26 27 28 29
#include <vector> 
#include <memory> 
#include <list> 

30
#include "AMDiS_fwd.h"
31 32 33 34 35 36 37 38 39 40 41
#include "FixVec.h"
#include "Global.h" 
#include "Flag.h" 
#include "RCNeighbourList.h" 
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "DOFContainer.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
Thomas Witkowski's avatar
Thomas Witkowski committed
42
#include "BasisFunction.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
43
#include "FiniteElemSpace.h"
44
#include "SurfaceQuadrature.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
45 46 47
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/ParallelDofMapping.h"
#endif
48

49 50
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
51 52
  using namespace std;

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

58 59
    DOFVectorBase() 
      : feSpace(NULL),
60
	elementVector(3),
Thomas Witkowski's avatar
Thomas Witkowski committed
61 62
        boundaryManager(NULL),
        nBasFcts(0)
63 64
    {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
65
      dofMap = NULL;
66 67
#endif
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
68
    
Thomas Witkowski's avatar
Thomas Witkowski committed
69
    DOFVectorBase(const FiniteElemSpace *f, string n);
70

Thomas Witkowski's avatar
Thomas Witkowski committed
71
    virtual ~DOFVectorBase();
72

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

Thomas Witkowski's avatar
Thomas Witkowski committed
78 79
    /// Evaluates the DOF vector at a set of quadrature points defined on the 
    /// given element.
80 81 82 83
    void getVecAtQPs(const ElInfo *elInfo, 
		     const Quadrature *quad,
		     const FastQuadrature *quadFast,
		     mtl::dense_vector<T>& vecAtQPs) const;
84

85 86 87 88 89
    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
90

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

Thomas Witkowski's avatar
Thomas Witkowski committed
126
    /// Returns the FE space the DOF vector is defined on.
127
    inline const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
128
    {
129
      return feSpace;
130
    }
131

Thomas Witkowski's avatar
Thomas Witkowski committed
132 133
    /// 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
134 135 136 137 138 139 140 141 142
    void assemble(T factor, ElInfo *elInfo,			    
		  const BoundaryType *bound, 
		  Operator *op = NULL);

    void assemble2(T factor, 
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
		   const BoundaryType *bound, 
		   Operator *op = NULL);
143 144

    void addElementVector(T sign,
145
			  const ElementVector& elVec, 
146
			  const BoundaryType *bound,
147
			  ElInfo *elInfo,
148
			  bool add = true); 
Thomas Witkowski's avatar
Thomas Witkowski committed
149

Thomas Witkowski's avatar
Thomas Witkowski committed
150 151 152
    /// 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
153
    void finishAssembling();
154 155 156 157 158 159 160 161
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
162
    }
163

Thomas Witkowski's avatar
Thomas Witkowski committed
164
    inline vector<double*>::iterator getOperatorFactorBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
165
    {
166
      return operatorFactor.begin();
167
    }
168

Thomas Witkowski's avatar
Thomas Witkowski committed
169
    inline vector<double*>::iterator getOperatorFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
170
    {
171
      return operatorFactor.end();
172
    }
173

Thomas Witkowski's avatar
Thomas Witkowski committed
174
    inline vector<double*>::iterator getOperatorEstFactorBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
175
    {
176
      return operatorEstFactor.begin();
177
    }
178

Thomas Witkowski's avatar
Thomas Witkowski committed
179
    inline vector<double*>::iterator getOperatorEstFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
180
    {
181
      return operatorEstFactor.end();
182
    }
183 184


Thomas Witkowski's avatar
Thomas Witkowski committed
185
    inline vector<Operator*>::iterator getOperatorsBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
186
    {
187
      return operators.begin();
188
    }
189

Thomas Witkowski's avatar
Thomas Witkowski committed
190
    inline vector<Operator*>::iterator getOperatorsEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
191
    {
192
      return operators.end();
193
    }
194 195 196 197 198 199 200 201 202 203 204

    Flag getAssembleFlag();

    /** \brief
     * Evaluates \f[ u_h(x(\lambda)) = \sum_{i=0}^{m-1} vec[ind[i]] * 
     * \varphi^i(\lambda) \f] where \f$ \varphi^i \f$ is the i-th basis function,
     * \f$ x(\lambda) \f$ are the world coordinates of lambda and
     * \f$ m \f$ is the number of basis functions
     */
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

Thomas Witkowski's avatar
Thomas Witkowski committed
205
    inline vector<Operator*>& getOperators() 
Thomas Witkowski's avatar
Thomas Witkowski committed
206
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
207
      return operators; 
208
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
209

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

Thomas Witkowski's avatar
Thomas Witkowski committed
215
    inline vector<double*>& getOperatorEstFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
216
    { 
217
      return operatorEstFactor; 
218
    }
219

Thomas Witkowski's avatar
Thomas Witkowski committed
220
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
221
    inline string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
222
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
223
      return name; 
224
    } 
225

Thomas Witkowski's avatar
Thomas Witkowski committed
226
    inline void setName(string n)
227
    {
228
      name = n;
229 230
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
231 232
    inline BoundaryManager* getBoundaryManager() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
233
      return boundaryManager; 
234
    }
235

Thomas Witkowski's avatar
Thomas Witkowski committed
236 237
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
238
      boundaryManager = bm;
239
    }
240

Thomas Witkowski's avatar
Thomas Witkowski committed
241
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
242
    void setDofMap(FeSpaceDofMap& m)
Thomas Witkowski's avatar
Thomas Witkowski committed
243
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
244
      dofMap = &m;
Thomas Witkowski's avatar
Thomas Witkowski committed
245 246 247 248
    }

    inline bool isRankDof(DegreeOfFreedom dof)
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
249
      TEST_EXIT_DBG(dofMap)("No rank dofs set!\n");
250

251
      return dofMap->isRankDof(dof, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
252 253 254
    }
#endif

255
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
256
    ///
257 258
    const FiniteElemSpace *feSpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
259
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
260
    string name;
261

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    ///
263
    ElementVector elementVector;
264

Thomas Witkowski's avatar
Thomas Witkowski committed
265
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
266
    vector<Operator*> operators;
267

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    vector<double*> operatorFactor;
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
272
    vector<double*> operatorEstFactor;
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
    ///
275 276
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
277
    /// Number of basis functions of the used finite element space.
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    int nBasFcts;
279

Thomas Witkowski's avatar
Thomas Witkowski committed
280
    /// Dimension of the mesh this DOFVectorBase belongs to
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    int dim;
Thomas Witkowski's avatar
Thomas Witkowski committed
282 283

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
284
    FeSpaceDofMap *dofMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
285
#endif
286

Thomas Witkowski's avatar
Thomas Witkowski committed
287
  };
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
  /// Specifies which operation should be done after coarsening
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320
  typedef enum{
    NO_OPERATION = 0,   
    COARSE_RESTRICT = 1,
    COARSE_INTERPOL = 2 
  } CoarsenOperation;
 

  /** \ingroup DOFAdministration 
   * \brief
   * The DOFs described above are just integers that can be used as indices into 
   * vectors and matrices. During refinement and coarsening of the mesh, the 
   * number of used DOFs, the meaning of one integer index, and even the total 
   * range of DOFs change. To be able to handle these changes automatically for 
   * all vectors, which are indexed by the DOFs, special data structures are 
   * used which contain such vector data. Lists of these structures are kept in 
   * DOFAdmin, so that all vectors in the lists can be resized together with the
   * range of DOFs. During refinement and coarsening of elements, values can be
   * interpolated automatically to new DOFs, and restricted from old DOFs.
   */
  template<typename T> 
  class DOFVector : public DOFVectorBase<T>, public Serializable
  {  
  public:
    /** \ingroup DOFAdministration
     * \brief
     * Enables the access of DOFVector<T>::Iterator. Alias for DOFIterator<T>
     */
    class Iterator : public DOFIterator<T> {
    public:
      Iterator(DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(c, type)
321
      {}
322 323 324

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
325
      {}
326 327 328 329
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
330 331 332
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}
333

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

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

    private:
      FiniteElemSpace *feSpace;
    };

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

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    /// Constructs a DOFVector with name n belonging to FiniteElemSpace f
Thomas Witkowski's avatar
Thomas Witkowski committed
356
    DOFVector(const FiniteElemSpace* f, string n); 
357

Thomas Witkowski's avatar
Thomas Witkowski committed
358
    /// Initialization.
Thomas Witkowski's avatar
Thomas Witkowski committed
359
    void init(const FiniteElemSpace* f, string n);
360

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

Thomas Witkowski's avatar
Thomas Witkowski committed
370
    /// Destructor
371 372
    virtual ~DOFVector();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
379
    /// Returns iterator to the end of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
380
    typename vector<T>::iterator end() 
Thomas Witkowski's avatar
Thomas Witkowski committed
381
    { 
382
      return vec.end(); 
383
    }
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,
Thomas Witkowski's avatar
Thomas Witkowski committed
388
				    vector<DegreeOfFreedom> &newDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
389 390

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

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

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

405
    /// Interpolation after refinement. Implemented for DOFVector<double>
406
    inline void refineInterpol(RCNeighbourList&, int) {}
407

Thomas Witkowski's avatar
Thomas Witkowski committed
408
    /// Returns \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
409
    vector<T>& getVector() 
Thomas Witkowski's avatar
Thomas Witkowski committed
410
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
411
      return vec;
412
    }
413

Thomas Witkowski's avatar
Thomas Witkowski committed
414
    /// Returns size of \ref vec
Thomas Witkowski's avatar
Thomas Witkowski committed
415 416
    inline int getSize() const 
    { 
417
      return vec.size();
418
    } 
419

Thomas Witkowski's avatar
Thomas Witkowski committed
420
    /// Returns used size of the vector.
Thomas Witkowski's avatar
Thomas Witkowski committed
421 422
    inline int getUsedSize() const 
    { 
423
      return this->feSpace->getAdmin()->getUsedSize(); 
424
    }
425

Thomas Witkowski's avatar
Thomas Witkowski committed
426
    /// Resizes \ref vec to n
Thomas Witkowski's avatar
Thomas Witkowski committed
427 428
    inline void resize(int n) 
    { 
429
      FUNCNAME("DOFVector<T>::resize()");
430
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
431
      vec.resize(n);
432
    } 
433

Thomas Witkowski's avatar
Thomas Witkowski committed
434
    /// Resizes \ref vec to n and inits new values with init
Thomas Witkowski's avatar
Thomas Witkowski committed
435 436
    inline void resize(int n, T init) 
    { 
437
      FUNCNAME("DOFVector<T>::resize()");
438
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); 
439
      vec.resize(n, init);
440
    } 
441

Thomas Witkowski's avatar
Thomas Witkowski committed
442
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
443 444
    inline const T& operator[](DegreeOfFreedom i) const 
    {
445
      FUNCNAME("DOFVector<T>::operator[]");
446
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
447
	("Illegal vector index %d.\n", i);
448
      return vec[i];
449
    } 
450

Thomas Witkowski's avatar
Thomas Witkowski committed
451
    /// Returns \ref vec[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
452 453
    inline T& operator[](DegreeOfFreedom i) 
    {
454
      FUNCNAME("DOFVector<T>::operator[]");
455 456

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

459
      return vec[i];
460
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
461
 
Thomas Witkowski's avatar
Thomas Witkowski committed
462
    /// Calculates Integral of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479
    double Int(Quadrature* q = NULL) const
    {
      return Int(-1, q);
    }

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

480

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
500
    /// Calculates L1 norm of this DOFVector
501
    double L1Norm(Quadrature* q = NULL) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
502
 
Thomas Witkowski's avatar
Thomas Witkowski committed
503
    /// Calculates L2 norm of this DOFVector
Thomas Witkowski's avatar
Thomas Witkowski committed
504 505
    inline double L2Norm(Quadrature* q = NULL) const 
    {
506
      return sqrt(L2NormSquare());
507
    }
508

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

588

Thomas Witkowski's avatar
Thomas Witkowski committed
589 590 591 592 593 594
    /// Eval DOFVector at given point p. If oldElInfo != NULL the search for 
    /// the element, where p is inside, starts from oldElInfo. implemented for:
    /// double, WorldVector< double >
    inline const T& evalAtPoint(WorldVector<double> &p, 
				ElInfo *oldElInfo = NULL, 
				T* value = NULL) const 
595 596 597
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
Praetorius, Simon's avatar
Praetorius, Simon committed
598
      return *value;
599
    }
600

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

609
    /// Writes the data of the DOFVector to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
610
    void serialize(ostream &out) 
Thomas Witkowski's avatar
Thomas Witkowski committed
611
    {
612 613 614
      unsigned int size = vec.size();
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
615
    }
616

617
    /// Reads data of an DOFVector from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
618
    void deserialize(istream &in) 
Thomas Witkowski's avatar
Thomas Witkowski committed
619
    {
620 621 622 623
      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));
624
    }
625

626 627
    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
628 629 630

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

631 632
    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
633 634

  protected: 
635

636
    /// Data container
Thomas Witkowski's avatar
Thomas Witkowski committed
637
    vector<T> vec; 
638
 
639
    /// Specifies what operation should be performed after coarsening
640 641 642
    CoarsenOperation coarsenOperation;
  }; 

Thomas Witkowski's avatar
Thomas Witkowski committed
643

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

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

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

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

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

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

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

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


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

Thomas Witkowski's avatar
Thomas Witkowski committed
699 700
    /// Implements DOFContainer::operator[]() by calling 
    /// DOFVector<DegreeOfFreedom>::operator[]()
Thomas Witkowski's avatar
Thomas Witkowski committed
701 702
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
703
      return DOFVector<DegreeOfFreedom>::operator[](i);
704
    }
705

Thomas Witkowski's avatar
Thomas Witkowski committed
706 707
    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
708 709 710
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

711
    /// Implements DOFIndexedBase::getSize()
Thomas Witkowski's avatar
Thomas Witkowski committed
712 713
    int getSize() const 
    {
714
      return DOFVector<DegreeOfFreedom>::getSize();
715
    }
716

717
    /// Implements DOFIndexedBase::resize()
Thomas Witkowski's avatar
Thomas Witkowski committed
718 719
    void resize(int size) 
    {
720 721 722 723 724 725 726 727 728
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
729

730
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
731 732
  double norm(DOFVector<T> *vec) 
  {
733
    return vec->nrm2();
734
  }
735 736

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
737 738
  double L2Norm(DOFVector<T> *vec) 
  {
739
    return vec->L2Norm();
740
  }
741 742

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