DOFVector.hh 32.2 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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
    createTempData();
  }
  

  template<typename T>
  DOFVectorBase<T>::~DOFVectorBase()
  {
    for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
      delete [] localIndices[i];
      delete [] localVectors[i];
      delete grdPhis[i];
      delete grdTmp[i];
      delete D2Phis[i];
    }
  }


  template<typename T>
  void DOFVectorBase<T>::createTempData()
  {
100
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
101
    localVectors.resize(omp_get_overall_max_threads());
102
    grdPhis.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
103
    grdTmp.resize(omp_get_overall_max_threads());
104
    D2Phis.resize(omp_get_overall_max_threads());
105

106
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
107
108
      localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
      localVectors[i] = new T[this->nBasFcts];
109
110
111
      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);
112
    }    
113
  }
114

Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
  
  template<typename T>
117
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
118
    : DOFVectorBase<T>(f, n),
119
      coarsenOperation(NO_OPERATION)
120
  {
121
    vec.resize(0);
122
123
124
    init(f, n);
  } 

125

126
  template<typename T>
127
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
128
  {
129
130
131
132
133
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
    this->boundaryManager = new BoundaryManager(f);
134
135
136
137
138
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
139
140
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
141

Thomas Witkowski's avatar
Thomas Witkowski committed
142
    if (this->boundaryManager)
143
      delete this->boundaryManager;
144
145

    vec.clear();
146
147
148
149
150
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153
154
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
155
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
156
157
158
159
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

160
161
162
    std::vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
163

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

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

171
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
172
	  (*this)[irow] += factor * elVec[i];
173
174
	else  
	  (*this)[irow] = factor * elVec[i];	
Thomas Witkowski's avatar
Thomas Witkowski committed
175
176
177
178
      }
    }
  }

179
180
181
182
183
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

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

186
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
189
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
190
191
192
193
194
195
196
197
198
199
      nrm += (*vecIterator) * (*vecIterator);

    return(sqrt(nrm));
  }

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

200
    checkFeSpace(this->feSpace, vec);
201
202
    
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
203
204
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
205
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
206
207
208
209
210
211
212
213
      nrm += (*vecIterator) * (*vecIterator);

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
214
    FUNCNAME("DOFVector<T>::asum()");
215

216
    checkFeSpace(this->feSpace, vec);
217

218
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
219
220
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
221
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
222
223
224
225
226
227
228
229
      nrm += abs(*vecIterator);

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
230
    FUNCNAME("DOFVector<T>::sum()");
231

232
    checkFeSpace(this->feSpace, vec);
233

234
    double nrm = 0.0;    
Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
237
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
238
239
240
241
242
243
244
245
      nrm += *vecIterator;

    return(nrm);
  }

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

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

250
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
251
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
252
253
254
255
256
257
258
      *vecIterator = alpha ; 
  }


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

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

263
264
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >= 
		  this->feSpace->getAdmin()->getUsedSize())
265
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
266
       this->feSpace->getAdmin()->getUsedSize());
267
    
268
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270
271
272
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), 
		       USED_DOFS);
    for (vecIterator.reset(), xIterator.reset(); !vecIterator.end();
	 ++vecIterator, ++xIterator)      
273
      *vecIterator = *xIterator; 
274
275
276
277
278
279
  }


  template<typename T>
  T DOFVector<T>::min() const
  {
280
281
    FUNCNAME("DOFVector<T>::min()");
    
282
    checkFeSpace(this->feSpace, vec);
283
284

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

289
290
291
292
293
294
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
295
296
    FUNCNAME("DOFVector<T>::max()");
    
297
    checkFeSpace(this->feSpace, vec);
298

299
    T m;    
300
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
301
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
302
      m = std::max(m, *vecIterator);
303

304
305
306
    return m;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
307
308
309
310
311
312
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

313
314
315
  template<typename T>
  void DOFVector<T>::print() const
  {
316
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
317

318
    const DOFAdmin *admin = NULL;
319
    const char *format;
320

321
322
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
323

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

364
365
366
367
368
369
370
371
372
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
373
374
375
376
377

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

382
    for (int i = 0; i < nBasisFcts; i++)
383
384
385
386
387
388
389
390
      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
391
    FUNCNAME("DOFVector::interpol()");
392

Thomas Witkowski's avatar
Thomas Witkowski committed
393
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
394
395

    interFct = fct;
396

397
    if (!this->getFeSpace()) {
398
399
400
401
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
402

403
    if (!(this->getFeSpace()->getAdmin())) {
404
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
405
	  this->getFeSpace()->getName().c_str());
406
407
      return;
    }
408
  
409
    if (!(this->getFeSpace()->getBasisFcts())) {
410
411
412
413
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
414

415
416
417
418
419
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
420
421
422

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
423
    TraverseStack stack;
424
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
427
428
429
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
430
431
432
  }

  template<typename T>
433
  void DOFVector<T>::interpolFct(ElInfo* elinfo)
434
  {
435
436
    const BasisFunction *basFct = traverseVector->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFeSpace()->getAdmin();
437
438
439
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
440
441
442

    DegreeOfFreedom *myLocalIndices = this->localIndices[omp_get_thread_num()];
    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), admin, myLocalIndices);
443
444
    int nBasFcts = basFct->getNumber();
    for (int i = 0; i < nBasFcts; i++)
445
      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
446
447
448
449
450
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
455
    if (!q) {
456
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
457
458
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
459

460
    FastQuadrature *quadFast = 
461
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
462

463
464
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
Thomas Witkowski's avatar
Thomas Witkowski committed
465
    std::vector<T> uh_vec(nPoints);
466
467
468
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
469
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
470
471
472
473
474
475
476
    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;
477

478
479
      elInfo = stack.traverseNext(elInfo);
    }
480

481
    return result;  
482
483
484
485
486
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
491
    if (!q) {
492
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
493
494
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
495

496
497
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
498

499
500
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
501
    std::vector<T> uh_vec(nPoints);
502
503
504
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
505
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
506
507
508
509
510
511
512
    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;
513

514
515
      elInfo = stack.traverseNext(elInfo);
    }
516

517
    return result;  
518
519
520
521
522
  }

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

525
526
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
527
    if (!q) {
528
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
529
530
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
531

532
533
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
534

535
536
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
537
    std::vector<T> uh_vec(nPoints);
538
539
540
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
541
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
542
543
544
545
546
547
548
    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;
549

550
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
551
    }
552

553
    return result;  
554
555
  }

556
  template<typename T>  
557
558
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
559
    FUNCNAME("DOFVector::H1NormSquare()");
560

561
562
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
563
    if (!q) {
564
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
565
566
567
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

568
569
570
571
572
573
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
574
    std::vector<WorldVector<T> > grduh_vec(nPoints);
575
576
577
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
578
579
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
580
581
582
583
584
585
586
587
588
589
590
    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;
      }
591

592
      result += det * normT;
593

594
595
      elInfo = stack.traverseNext(elInfo);
    }
596

597
    return result;  
598
599
600
601
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
602
					std::vector<DegreeOfFreedom> &newDOF)
603
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
604
605
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
606
	vec[newDOF[i]] = vec[i];
607
608
609
610
611
612
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
613
614
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
615
	 op != this->operators.end(); ++op)
616
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
617

618
619
620
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
621
622
623
624
625
  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
626
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
627
628
629
      (*it)->finishAssembling();
  }

630
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
631
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
632
  {
633
    this->feSpace = rhs.feSpace;
634
    this->dim = rhs.dim;
635
    this->nBasFcts = rhs.nBasFcts;
636
    vec = rhs.vec;
637
    this->elementVector.change_dim(this->nBasFcts);
638
639
640
641
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
642
    
643
    if (rhs.boundaryManager) {
644
645
646
      if (this->boundaryManager) 
	delete this->boundaryManager; 

647
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
648
    } else {
649
      this->boundaryManager = NULL;
650
651
    }

652
653
654
655
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

656
657
658
659
660
661
    return *this;
  }

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

664
665
    TEST_EXIT_DBG(x.getFeSpace() && x.getFeSpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFeSpace(), x.getFeSpace()->getAdmin());
666

Thomas Witkowski's avatar
Thomas Witkowski committed
667
668
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
669
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
670
      (*vecIterator) *= scal; 
671

672
673
674
675
676
677
678
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
679
680
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
681
682
683
684
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
685
      ("no admin or different admins: %8X, %8X\n",
686
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
687
688
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
689
690
    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
691
692
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
693
694
      *xIterator += *yIterator; 

695
696
697
698
699
700
    return x;
  }

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

703
704
705
706
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
707
      ("no admin or different admins: %8X, %8X\n",
708
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
709
710
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
711
712
    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
713
714
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
715
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
716

717
718
719
720
721
722
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
723
724
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
725
726
727
728
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
729
      ("no admin or different admins: %8X, %8X\n",
730
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
731
732
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
733
734
    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
735
736
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
737
      *xIterator *= *yIterator; 
738
739
740
741
742
743
744
745

    return x;
  }

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

748
749
750
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT((admin = x.getFeSpace()->getAdmin()) && (admin == y.getFeSpace()->getAdmin()))
751
      ("no admin or different admins: %8X, %8X\n",
752
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
753
754
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
755
    T dot = 0;
756
757
758

    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
759
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
760
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
761
      dot += (*xIterator) * (*yIterator);
762

Thomas Witkowski's avatar
Thomas Witkowski committed
763
    return dot;
764
765
766
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
767
768
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
769
  {
770
771
772
    // Unfortunately needed
    using namespace mtl;

773
774
    FUNCNAME("DOFVector<T>::mv()");

775
776
777
778
    TEST_EXIT(a.getRowFeSpace() && a.getColFeSpace() && 
	      x.getFeSpace() && result.getFeSpace())
      ("getFeSpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFeSpace(), a.getColFeSpace(), x.getFeSpace(), result.getFeSpace());
779

780
781
    const DOFAdmin *rowAdmin = a.getRowFeSpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFeSpace()->getAdmin();
782
783

    TEST_EXIT((rowAdmin && colAdmin) && 
784
785
786
787
	      (((transpose == NoTranspose) && (colAdmin == x.getFeSpace()->getAdmin()) && 
		(rowAdmin == result.getFeSpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFeSpace()->getAdmin()) && 
		(colAdmin == result.getFeSpace()->getAdmin()))))
788
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
789
790
       a.getRowFeSpace()->getAdmin(), a.getColFeSpace()->getAdmin(), 
       x.getFeSpace()->getAdmin(), result.getFeSpace()->getAdmin());
791
792


793
794
795
796
797
    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 
798
799
800
801
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
802
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
803
804
805
  {
    FUNCNAME("DOFVector<T>::axpy()");

806
807
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
808

809
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
810

811
    TEST_EXIT((admin) && (admin == y.getFeSpace()->getAdmin()))
812
      ("no admin or different admins: %8X, %8X\n",
813
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
814
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
815
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
816
817
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
818
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
       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
838
839
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
840
	  y[i] += alpha * x[i];
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
  }

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

872
873
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
874

875
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
876

877
    TEST_EXIT(admin && (admin == y.getFeSpace()->getAdmin()))
878
      ("no admin or different admins: %8X, %8X\n",
879
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
880
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
881
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
882
883
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
884
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903