DOFVector.hh 31.9 KB
Newer Older
1
#include <list>
Thomas Witkowski's avatar
Thomas Witkowski committed
2
#include <algorithm>
3
#include <math.h>
4

5
6
7
8
9
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/utility/category.hpp>
#include <boost/numeric/linear_algebra/identity.hpp>

10
11
12
13
14
15
16
17
18
19
20
21
#include "FixVec.h"
#include "Boundary.h"
#include "DOFAdmin.h"
#include "ElInfo.h"
#include "Error.h"
#include "FiniteElemSpace.h"
#include "Global.h"
#include "Mesh.h"
#include "Quadrature.h"
#include "AbstractFunction.h"
#include "BoundaryManager.h"
#include "Assembler.h"
22
#include "OpenMP.h"
23
#include "Operator.h"
24
#include "Parameters.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
25
#include "Traverse.h"
26

27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
// Defining the interface for MTL4
namespace mtl {

  // Let MTL4 know that DOFVector it is a column vector
  namespace traits {
    template <typename T>
    struct category< AMDiS::DOFVector<T> > 
    {
      typedef tag::dense_col_vector type;
    };
  }

  namespace ashape {
    template <typename T>
    struct ashape< AMDiS::DOFVector<T> > 
    {
      typedef cvec<typename ashape<T>::type> type;
    };
  }

  // Modelling Collection and MutableCollection
  template <typename T>
  struct Collection< AMDiS::DOFVector<T> > 
  {
    typedef T               value_type;
    typedef const T&        const_reference;
    typedef std::size_t     size_type;
  };

  template <typename T>
  struct MutableCollection< AMDiS::DOFVector<T> > 
    : public Collection< AMDiS::DOFVector<T> > 
  {
    typedef T&              reference;
  };


} // namespace mtl



68
69
70
namespace AMDiS {

  template<typename T>
71
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
72
73
    : feSpace(f),
      name(n),
74
      elementVector(f->getBasisFcts()->getNumber()),
Thomas Witkowski's avatar
Thomas Witkowski committed
75
      boundaryManager(NULL)
76
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
77
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
78
    dim = feSpace->getMesh()->getDim();
79

80
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
81
    localVectors.resize(omp_get_overall_max_threads());
82
    grdPhis.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
83
    grdTmp.resize(omp_get_overall_max_threads());
84
    D2Phis.resize(omp_get_overall_max_threads());
85

86
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
87
88
      localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
      localVectors[i] = new T[this->nBasFcts];
89
90
91
      grdPhis[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      grdTmp[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      D2Phis[i] = new DimMat<double>(dim, NO_INIT);
92
93
    }
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
94
  
95
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
96
97
98
  DOFVectorBase<T>::~DOFVectorBase()
  {
    for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
99
100
      delete [] localIndices[i];
      delete [] localVectors[i];
101
102
103
      delete grdPhis[i];
      delete grdTmp[i];
      delete D2Phis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
104
105
106
107
    }
  }
  
  template<typename T>
108
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
109
    : DOFVectorBase<T>(f, n),
110
      coarsenOperation(NO_OPERATION)
111
  {
112
    vec.resize(0);
113
114
115
116
    init(f, n);
  } 

  template<typename T>
117
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
118
  {
119
120
121
122
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
123
    this->boundaryManager = new BoundaryManager(f);
124
125
126
127
128
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
129
130
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
131

Thomas Witkowski's avatar
Thomas Witkowski committed
132
    if (this->boundaryManager)
133
      delete this->boundaryManager;
134
135

    vec.clear();
136
137
138
139
140
  }

  template<typename T>
  DOFVector<T> * DOFVector<T>::traverseVector = NULL;

Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
143
144
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
145
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
146
147
148
149
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

150
151
152
153
    Vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
						feSpace->getAdmin(),
						&indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
154

155
    for (DegreeOfFreedom i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
156
157
158
      BoundaryCondition *condition = 
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;

Thomas Witkowski's avatar
Thomas Witkowski committed
159
      if (!(condition && condition->isDirichlet())) {
160
	DegreeOfFreedom irow = indices[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
161
162
163
164
165
166

	if (add) {
	  (*this)[irow] += factor * elVec[i];
	} else {  
	  (*this)[irow] = factor * elVec[i];
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
169
170
      }
    }
  }

171
172
173
174
175
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

176
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
177

178
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
181
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
182
183
184
185
186
187
188
189
190
191
      nrm += (*vecIterator) * (*vecIterator);

    return(sqrt(nrm));
  }

  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

192
    checkFeSpace(this->feSpace, vec);
193
194
    
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
197
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
198
199
200
201
202
203
204
205
      nrm += (*vecIterator) * (*vecIterator);

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
206
    FUNCNAME("DOFVector<T>::asum()");
207

208
    checkFeSpace(this->feSpace, vec);
209

210
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
211
212
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
213
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
214
215
216
217
218
219
220
221
      nrm += abs(*vecIterator);

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
222
    FUNCNAME("DOFVector<T>::sum()");
223

224
    checkFeSpace(this->feSpace, vec);
225

226
    double nrm = 0.0;    
Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
229
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
230
231
232
233
234
235
236
237
      nrm += *vecIterator;

    return(nrm);
  }

  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
238
    FUNCNAME("DOFVector<T>::set()");
239

240
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
241

242
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
243
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
244
245
246
247
248
249
250
      *vecIterator = alpha ; 
  }


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
251
    FUNCNAME("DOFVector<T>::copy()");
Thomas Witkowski's avatar
Thomas Witkowski committed
252

253
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
254

255
256
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >= 
		  this->feSpace->getAdmin()->getUsedSize())
257
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
258
       this->feSpace->getAdmin()->getUsedSize());
259
    
260
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
263
264
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), 
		       USED_DOFS);
    for (vecIterator.reset(), xIterator.reset(); !vecIterator.end();
	 ++vecIterator, ++xIterator)      
265
      *vecIterator = *xIterator; 
266
267
268
269
270
271
  }


  template<typename T>
  T DOFVector<T>::min() const
  {
272
273
    FUNCNAME("DOFVector<T>::min()");
    
274
    checkFeSpace(this->feSpace, vec);
275
276

    T m;    
277
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
279
      m = std::min(m, *vecIterator);
280

281
282
283
284
285
286
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
287
288
    FUNCNAME("DOFVector<T>::max()");
    
289
    checkFeSpace(this->feSpace, vec);
290

291
    T m;    
292
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
293
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
294
      m = std::max(m, *vecIterator);
295

296
297
298
    return m;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
301
302
303
304
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

305
306
307
  template<typename T>
  void DOFVector<T>::print() const
  {
308
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
309

310
    const DOFAdmin *admin = NULL;
311
    const char *format;
312

313
314
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
315

316
    MSG("Vec `%s':\n", this->name.c_str());
317
    int j = 0;
318
319
320
321
322
323
324
325
326
    if (admin) {
      if (admin->getUsedSize() > 100)
	format = "%s(%3d,%10.5le)";
      else if (admin->getUsedSize() > 10)
	format = "%s(%2d,%10.5le)";
      else
	format = "%s(%1d,%10.5le)";

      Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
327
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
328
	if ((j % 3) == 0) {
329
330
	  if (j) 
	    Msg::print("\n");
331
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
332
	} else {
333
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
334
	}
335
	j++;
336
      }
337
      Msg::print("\n");
338
    } else {
339
340
	MSG("no DOFAdmin, print whole vector.\n");
    
341
342
343
344
345
346
	for (int i = 0; i < static_cast<int>( vec.size()); i++) {
	  if ((j % 3) == 0) {
	    if (j) 
	      Msg::print("\n");
	    MSG("(%d,%10.5e)",i,vec[i]);
	  } else {
347
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
348
	  }
349
350
351
352
353
354
355
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }

356
357
358
359
360
361
362
363
364
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
365
366
367
368
369
370
371
372
373

  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda, 
			     DegreeOfFreedom* dof_indices)
  {
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFESpace()->getBasisFcts());
    int numberOfBasFcts = phi->getNumber();
    T val = 0.0;

374
    for (int i = 0; i < numberOfBasFcts; i++)
375
376
377
378
379
380
381
382
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
383
    FUNCNAME("interpol");
384

385
386
387
    TEST_EXIT_DBG(fct)("no function to interpolate\n");

    interFct = fct;
388

389
390
391
392
393
    if (!this->getFESpace()) {
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
394

395
396
397
398
399
    if (!(this->getFESpace()->getAdmin())) {
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
	  this->getFESpace()->getName().c_str());
      return;
    }
400
  
401
402
403
404
405
    if (!(this->getFESpace()->getBasisFcts())) {
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
406

407
408
409
410
411
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
412
413
414

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
419
420
421
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this->getFESpace()->getMesh(), -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
422
423
424
  }

  template<typename T>
425
  void DOFVector<T>::interpolFct(ElInfo* elinfo)
426
427
428
  {
    const BasisFunction *basFct = traverseVector->getFESpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFESpace()->getAdmin();
429
430
431
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
432
433
434

    DegreeOfFreedom *myLocalIndices = this->localIndices[omp_get_thread_num()];
    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), admin, myLocalIndices);
435
436
    int nBasFcts = basFct->getNumber();
    for (int i = 0; i < nBasFcts; i++)
437
      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
438
439
440
441
442
  }

  template<typename T>
  double DOFVector<T>::Int(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
443
    FUNCNAME("DOFVector::Int()");
444
  
445
    Mesh* mesh = this->feSpace->getMesh();
446

Thomas Witkowski's avatar
Thomas Witkowski committed
447
    if (!q) {
448
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
451

452
    FastQuadrature *quadFast = 
453
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
454

455
456
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
Thomas Witkowski's avatar
Thomas Witkowski committed
457
    std::vector<T> uh_vec(nPoints);
458
459
460
461
462
463
464
465
466
467
468
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getVecAtQPs(elInfo, NULL, quadFast, &(uh_vec[0]));
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;
469

470
471
      elInfo = stack.traverseNext(elInfo);
    }
472

473
    return result;  
474
475
476
477
478
  }

  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
479
    FUNCNAME("DOFVector::L1Norm()");
480
  
481
482
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
483
    if (!q) {
484
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
487

488
489
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
490

491
492
493
494
495
496
497
498
499
500
501
502
503
504
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    std::vector<T> uh_vec(this->feSpace->getBasisFcts()->getNumber());
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getVecAtQPs(elInfo, NULL, quadFast, &(uh_vec[0]));
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
505

506
507
      elInfo = stack.traverseNext(elInfo);
    }
508

509
    return result;  
510
511
512
513
514
  }

  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
515
    FUNCNAME("DOFVector::L2NormSquare()");
516

517
518
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
519
    if (!q) {
520
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
523

524
525
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
526

527
528
529
530
531
532
533
534
535
536
537
538
539
540
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    std::vector<T> uh_vec(this->feSpace->getBasisFcts()->getNumber());
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getVecAtQPs(elInfo, NULL, quadFast, &(uh_vec[0]));
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;
541

542
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
543
    }
544

545
    return result;  
546
547
  }

548
  template<typename T>  
549
550
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
551
    FUNCNAME("DOFVector::H1NormSquare()");
552

553
554
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
555
    if (!q) {
556
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
557
558
559
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
    std::vector<WorldVector<T> > grduh_vec(this->feSpace->getBasisFcts()->getNumber());
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getGrdAtQPs(elInfo, NULL, quadFast, &(grduh_vec[0]));

      for (int iq = 0; iq < nPoints; iq++) {
	double norm2 = 0.0;
	for (int j = 0; j < dimOfWorld; j++)
	  norm2 += sqr(grduh_vec[iq][j]);
	normT += quadFast->getWeight(iq) * norm2;
      }
582

583
      result += det * normT;
584

585
586
      elInfo = stack.traverseNext(elInfo);
    }
587

588
    return result;  
589
590
591
592
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
593
					std::vector<DegreeOfFreedom> &newDOF)
594
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
595
596
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
597
	vec[newDOF[i]] = vec[i];
598
599
600
601
602
603
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
604
605
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
606
	 op != this->operators.end(); ++op)
607
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
608

609
610
611
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
612
613
614
615
616
  template<typename T>
  void DOFVectorBase<T>::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = this->operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
617
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
618
619
620
      (*it)->finishAssembling();
  }

621
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
622
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
623
  {
624
    this->feSpace = rhs.feSpace;
625
    this->nBasFcts = rhs.nBasFcts;
626
    vec = rhs.vec;
627
    this->elementVector.change_dim(this->nBasFcts);
628
629
630
631
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
632
    if (rhs.boundaryManager) {
633
634
635
      if (this->boundaryManager) 
	delete this->boundaryManager; 

636
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
637
    } else {
638
      this->boundaryManager = NULL;
639
640
    }

641
642
643
644
645
646
    return *this;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal)
  {
647
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, T scal)");
648

649
650
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
651

Thomas Witkowski's avatar
Thomas Witkowski committed
652
653
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
654
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
655
      (*vecIterator) *= scal; 
656

657
658
659
660
661
662
663
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
664
665
666
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
667
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
668
669
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
670
671
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
672
673
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
674
675
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
678
679
      *xIterator += *yIterator; 

680
681
682
683
684
685
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y)
  {
686
    FUNCNAME("DOFVector<T>::operator-=(DOFVector<T>& x, const DOFVector<T>& y)");
687

688
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
689
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
690
691
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
692
693
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
694
695
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
696
697
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
698
699
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
700
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
701

702
703
704
705
706
707
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
708
709
710
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
711
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
712
713
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
714
715
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
716
717
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
718
719
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
720
721
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
722
      *xIterator *= *yIterator; 
723
724
725
726
727
728
729
730

    return x;
  }

  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::operator*(DOFVector<T>& x, DOFVector<T>& y)");
731
    const DOFAdmin *admin = NULL;
732

733
    TEST_EXIT(x.getFESpace() && y.getFESpace())
734
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
735
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) && (admin == y.getFESpace()->getAdmin()))
736
737
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
738
739
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
740
    T dot = 0;
741
742
743

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
744
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
745
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
746
      dot += (*xIterator) * (*yIterator);
747

Thomas Witkowski's avatar
Thomas Witkowski committed
748
    return dot;
749
750
751
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
752
753
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
754
  {
755
756
757
    // Unfortunately needed
    using namespace mtl;

758
759
    FUNCNAME("DOFVector<T>::mv()");

Thomas Witkowski's avatar
Thomas Witkowski committed
760
761
    TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && 
	      x.getFESpace() && result.getFESpace())
762
763
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
764
765
766
767
768
769
770
771
772

    const DOFAdmin *rowAdmin = a.getRowFESpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFESpace()->getAdmin();

    TEST_EXIT((rowAdmin && colAdmin) && 
	      (((transpose == NoTranspose) && (colAdmin == x.getFESpace()->getAdmin()) && 
		(rowAdmin == result.getFESpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFESpace()->getAdmin()) && 
		(colAdmin == result.getFESpace()->getAdmin()))))
773
774
775
776
777
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());


778
779
780
781
782
    if (transpose == NoTranspose)
      mult(a.getBaseMatrix(), x, result);
    else if (transpose == Transpose)
      mult(trans(const_cast<DOFMatrix::base_matrix_type&>(a.getBaseMatrix())), x, result);
    else 
783
784
785
786
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
787
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
788
789
790
  {
    FUNCNAME("DOFVector<T>::axpy()");

791
    TEST_EXIT(x.getFESpace() && y.getFESpace())
792
793
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

794
795
796
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT((admin) && (admin == y.getFESpace()->getAdmin()))
797
798
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
799
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
800
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
801
802
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
803
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
       admin->getUsedSize());

      // This is the old implementation of the mv-multiplication. It have been changed
      // because of the OpenMP-parallelization:
//     typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
//     typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
//     for(xIterator.reset(), yIterator.reset();
// 	!xIterator.end();
// 	++xIterator, ++yIterator) 
//       {  
// 	*yIterator += alpha * (*xIterator); 
//       };

      int i;
      int maxI = y.getSize();

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
823
824
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
825
	  y[i] += alpha * x[i];
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
  }

  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d)
  {
    static DOFVector<T> result;
    return mult(d, v, result);
  }

  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v)
  {
    static DOFVector<T> result;
    return mult(d, v, result);
  }

  template<typename T>
  const DOFVector<T>& operator+(const DOFVector<T> &v1 , const DOFVector<T> &v2)
  {
    static DOFVector<T> result;
    return add(v1, v2, result);
  }


  template<typename T>
  void xpay(double alpha, 
	    const DOFVector<T>& x, 
	    DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::xpay()");

857
    TEST_EXIT(x.getFESpace() && y.getFESpace())
858
859
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

860
861
862
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT(admin && (admin == y.getFESpace()->getAdmin()))
863
864
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
865
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
866
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
867
868
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
869
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
       admin->getUsedSize());

    // This is the old implementation of the mv-multiplication. It have been changed
    // because of the OpenMP-parallelization:
    //     typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
    //     typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
    //     for(xIterator.reset(), yIterator.reset();
    // 	!xIterator.end();
    // 	++xIterator, ++yIterator) 
    //       {  
    // 	*yIterator = alpha *(*yIterator)+ (*xIterator); 
    //       };

      int i;
      int maxI = y.getSize();

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
889
890
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
891
	  y[i] = alpha * y[i] + x[i];
892
893
894
895
896
897
898
899
900
  }

  template<typename T>
  inline const DOFVector<T>& mult(double scal,
				  const DOFVector<T>& v,
				  DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
901
902
903
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = scal * (*vIterator);
904
905
906
907
908
909
910
911
912
913
914

    return result;
  }

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
915
916
917
918
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = (*vIterator) + scal;

919
920
921
922
923
924
925
926
927
928
929
    return result;
  }

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v1,
				 const DOFVector<T>& v2,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator v1Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v1)), USED_DOFS);
    typename DOFVector<T>::Iterator v2Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v2)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
930
931
932
    for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	 !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator) 
      *rIterator = (*v1Iterator) + (*v2Iterator);
933

Thomas Witkowski's avatar
Thomas Witkowski committed
934
    return result;
935
936
937
  }

  template<typename T>
938
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
939
  {
940
    static T* localVec = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
941
    static int localVecSize = 0;
942
943
944
    const DOFAdmin* admin = feSpace->getAdmin();
    T *result;
  
945
    if (d) {
946
947
      result = d;
    } else {
948
949
950
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
951
952
953
954
955
956
957
958
      if (localVec && nBasFcts > localVecSize)  {
	delete [] localVec;
	localVec = new T[nBasFcts]; 
      }
      if (!localVec)
	localVec = new T[nBasFcts]; 

      localVecSize = nBasFcts;
959
960
      result = localVec;
    }
961

962
963
964
965
#ifdef _OPENMP
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(el, admin, &(localIndices[0]));
#else
966
967
    const DegreeOfFreedom *localIndices = 
      feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);
968
#endif
969

970
    for (int i = 0; i < nBasFcts; i++)
971
      result[i] = (*this)[localIndices[i]];
972
973

    return result;
974
975
976
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
977
978
  const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo, 
					 const Quadrature *quad,
979
					 const FastQuadrature *quadFast,
Thomas Witkowski's avatar
Thomas Witkowski committed
980
					 T *vecAtQPs) const
981
982
983
  {
    FUNCNAME("DOFVector<T>::getVecAtQPs()");
  
984
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
985

Thomas Witkowski's avatar
Thomas Witkowski committed
986
    if (quad && quadFast)
987
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
988
	("quad != quadFast->quadrature\n");
989

990
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
991
992
      ("invalid basis functions");

993
994
    Element *el = elInfo->getElement();
    const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
995
    const BasisFunction *basFcts = feSpace->getBasisFcts();
996
    int nPoints = quadrature->getNumPoints();
997
    int nBasFcts  = basFcts->getNumber();
998
999
1000
1001
1002
1003
    static T *localvec = NULL;
    T *result;

    if (vecAtQPs) {
      result = vecAtQPs;
    } else {
1004
1005
1006
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
1007
1008
      if (localvec) 
	delete [] localvec;
1009
1010
      localvec = new T[nPoints];
      for (int i = 0; i < nPoints; i++)
1011
	localvec[i] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1012