DOFVector.hh 32 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
123
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
    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
    std::vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
153

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

Thomas Witkowski's avatar
Thomas Witkowski committed
158
      if (!(condition && condition->isDirichlet())) {
159
	DegreeOfFreedom irow = indices[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
160

161
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
162
	  (*this)[irow] += factor * elVec[i];
163
164
	else  
	  (*this)[irow] = factor * elVec[i];	
Thomas Witkowski's avatar
Thomas Witkowski committed
165
166
167
168
      }
    }
  }

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

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

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

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
204
    FUNCNAME("DOFVector<T>::asum()");
205

206
    checkFeSpace(this->feSpace, vec);
207

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

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
220
    FUNCNAME("DOFVector<T>::sum()");
221

222
    checkFeSpace(this->feSpace, vec);
223

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

    return(nrm);
  }

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

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

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


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

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

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


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

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

279
280
281
282
283
284
    return m;
  }

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

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

294
295
296
    return m;
  }

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

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

308
    const DOFAdmin *admin = NULL;
309
    const char *format;
310

311
312
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
313

314
    MSG("Vec `%s':\n", this->name.c_str());
315
    int j = 0;
316
317
318
319
320
321
322
323
324
    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
325
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
326
	if ((j % 3) == 0) {
327
328
	  if (j) 
	    Msg::print("\n");
329
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
330
	} else {
331
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
332
	}
333
	j++;
334
      }
335
      Msg::print("\n");
336
    } else {
337
338
	MSG("no DOFAdmin, print whole vector.\n");
    
339
340
341
342
343
344
	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 {
345
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
346
	  }
347
348
349
350
351
352
353
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }

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

    return result;
  }
363
364
365
366
367
368

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

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

    return val;
  }

  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
381
    FUNCNAME("DOFVector::interpol()");
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
384
385

    interFct = fct;
386

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

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

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

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
415
416
417
418
419
    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);
    }
420
421
422
  }

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

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

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

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

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

453
454
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
Thomas Witkowski's avatar
Thomas Witkowski committed
455
    std::vector<T> uh_vec(nPoints);
456
457
458
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
459
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
460
461
462
463
464
465
466
    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;
467

468
469
      elInfo = stack.traverseNext(elInfo);
    }
470

471
    return result;  
472
473
474
475
476
  }

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

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

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

489
490
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
491
    std::vector<T> uh_vec(nPoints);
492
493
494
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
495
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
496
497
498
499
500
501
502
    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;
503

504
505
      elInfo = stack.traverseNext(elInfo);
    }
506

507
    return result;  
508
509
510
511
512
  }

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

515
516
    Mesh* mesh = this->feSpace->getMesh();

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

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

525
526
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
527
    std::vector<T> uh_vec(nPoints);
528
529
530
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
531
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
532
533
534
535
536
537
538
    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;
539

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

543
    return result;  
544
545
  }

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

551
552
    Mesh* mesh = this->feSpace->getMesh();

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

558
559
560
561
562
563
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
564
    std::vector<WorldVector<T> > grduh_vec(nPoints);
565
566
567
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
568
569
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
570
571
572
573
574
575
576
577
578
579
580
    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;
      }
581

582
      result += det * normT;
583

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

587
    return result;  
588
589
590
591
  }

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

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

608
609
610
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
611
612
613
614
615
  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
616
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
617
618
619
      (*it)->finishAssembling();
  }

620
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
621
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
622
  {
623
    this->feSpace = rhs.feSpace;
624
    this->nBasFcts = rhs.nBasFcts;
625
    vec = rhs.vec;
626
    this->elementVector.change_dim(this->nBasFcts);
627
628
629
630
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
631
    
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
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

645
646
647
648
649
650
    return *this;
  }

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

653
654
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
655

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

661
662
663
664
665
666
667
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
668
669
670
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
671
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
672
673
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
674
675
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
676
677
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
678
679
    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
680
681
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
682
683
      *xIterator += *yIterator; 

684
685
686
687
688
689
    return x;
  }

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

692
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
693
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
694
695
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
696
697
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
698
699
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
700
701
    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
702
703
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
704
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
705

706
707
708
709
710
711
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
712
713
714
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
715
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
716
717
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
718
719
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
720
721
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
722
723
    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
724
725
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
726
      *xIterator *= *yIterator; 
727
728
729
730
731
732
733
734

    return x;
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
744
    T dot = 0;
745
746
747

    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
748
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
749
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
750
      dot += (*xIterator) * (*yIterator);
751

Thomas Witkowski's avatar
Thomas Witkowski committed
752
    return dot;
753
754
755
  }

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

762
763
    FUNCNAME("DOFVector<T>::mv()");

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

    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()))))
777
778
779
780
781
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());


782
783
784
785
786
    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 
787
788
789
790
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
791
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
792
793
794
  {
    FUNCNAME("DOFVector<T>::axpy()");

795
    TEST_EXIT(x.getFESpace() && y.getFESpace())
796
797
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

798
799
800
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT((admin) && (admin == y.getFESpace()->getAdmin()))
801
802
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
803
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
804
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
805
806
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
807
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
       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
827
828
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
829
	  y[i] += alpha * x[i];
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
857
858
859
860
  }

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

861
    TEST_EXIT(x.getFESpace() && y.getFESpace())
862
863
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

864
865
866
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT(admin && (admin == y.getFESpace()->getAdmin()))
867
868
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
869
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
870
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
871
872
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
873
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
       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
893
894
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
895
	  y[i] = alpha * y[i] + x[i];
896
897
898
899
900
901
902
903
904
  }

  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
905
906
907
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = scal * (*vIterator);
908
909
910
911
912
913
914
915
916
917
918

    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
919
920
921
922
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = (*vIterator) + scal;

923
924
925
    return result;
  }

926

927
  template<typename T>