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
  {
77
    /*
Thomas Witkowski's avatar
Thomas Witkowski committed
78
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
79
    dim = feSpace->getMesh()->getDim();
80

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

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

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

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

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

    vec.clear();
138
139
140
141
142
  }

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

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

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

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

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

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

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

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

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

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
208
    FUNCNAME("DOFVector<T>::asum()");
209

210
    checkFeSpace(this->feSpace, vec);
211

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

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
224
    FUNCNAME("DOFVector<T>::sum()");
225

226
    checkFeSpace(this->feSpace, vec);
227

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

    return(nrm);
  }

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

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

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


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

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

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


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

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

283
284
285
286
287
288
    return m;
  }

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

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

298
299
300
    return m;
  }

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

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

312
    const DOFAdmin *admin = NULL;
313
    const char *format;
314

315
316
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
317

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

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

    return result;
  }
367
368
369
370
371
372

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

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

    return val;
  }

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

387
388
389
    TEST_EXIT_DBG(fct)("no function to interpolate\n");

    interFct = fct;
390

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

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

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

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
419
420
421
422
423
    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);
    }
424
425
426
  }

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

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

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

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

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

457
458
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
Thomas Witkowski's avatar
Thomas Witkowski committed
459
    std::vector<T> uh_vec(nPoints);
460
461
462
463
464
465
466
467
468
469
470
    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;
471

472
473
      elInfo = stack.traverseNext(elInfo);
    }
474

475
    return result;  
476
477
478
479
480
  }

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

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

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

493
494
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
495
    std::vector<T> uh_vec(nPoints);
496
497
498
499
500
501
502
503
504
505
506
    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;
507

508
509
      elInfo = stack.traverseNext(elInfo);
    }
510

511
    return result;  
512
513
514
515
516
  }

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

519
520
    Mesh* mesh = this->feSpace->getMesh();

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

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

529
530
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
531
    std::vector<T> uh_vec(nPoints);
532
533
534
535
536
537
538
539
540
541
542
    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;
543

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

547
    return result;  
548
549
  }

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

555
556
    Mesh* mesh = this->feSpace->getMesh();

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

562
563
564
565
566
567
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

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

586
      result += det * normT;
587

588
589
      elInfo = stack.traverseNext(elInfo);
    }
590

591
    return result;  
592
593
594
595
  }

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

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

612
613
614
    return fillFlag;
  }

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

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

640
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
641
    } else {
642
      this->boundaryManager = NULL;
643
644
    }

645
646
647
648
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

649
650
651
652
653
654
    return *this;
  }

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

657
658
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
659

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

665
666
667
668
669
670
671
    return x;
  }


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

688
689
690
691
692
693
    return x;
  }

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

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

710
711
712
713
714
715
    return x;
  }

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

    return x;
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
748
    T dot = 0;
749
750
751

    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
752
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
753
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
754
      dot += (*xIterator) * (*yIterator);
755

Thomas Witkowski's avatar
Thomas Witkowski committed
756
    return dot;
757
758
759
  }

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

766
767
    FUNCNAME("DOFVector<T>::mv()");

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

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


786
787
788
789
790
    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 
791
792
793
794
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
795
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
796
797
798
  {
    FUNCNAME("DOFVector<T>::axpy()");

799
    TEST_EXIT(x.getFESpace() && y.getFESpace())
800
801
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

802
803
804
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

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

865
    TEST_EXIT(x.getFESpace() && y.getFESpace())
866
867
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

868
869
870
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

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

    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
923
924
925
926
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = (*vIterator) + scal;

927
928
929
930
931
932
933
934
935
936
937
    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
938
939
940
    for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	 !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator) 
      *rIterator = (*v1Iterator) + (*v2Iterator);
941

Thomas Witkowski's avatar
Thomas Witkowski committed
942
    return result;
943
944
945
  }

  template<typename T>
946
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
947
  {