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

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

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

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

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

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

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

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


} // namespace mtl



68
69
70
namespace AMDiS {

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

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

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

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

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

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

    vec.clear();
136
137
138
139
140
141
142
143
144
145
146
  }

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

  template<typename T>
  FastQuadrature *DOFVector<T>::quad_fast = NULL;

  template<typename T>
  double DOFVector<T>::norm = 0.0;

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

156
157
158
159
    Vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
						feSpace->getAdmin(),
						&indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
160

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

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

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

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

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

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

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
212
    FUNCNAME("DOFVector<T>::asum()");
213

214
    checkFeSpace(this->feSpace, vec);
215

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

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
228
    FUNCNAME("DOFVector<T>::sum()");
229

230
    checkFeSpace(this->feSpace, vec);
231

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

    return(nrm);
  }

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

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

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


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

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

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


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

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

287
288
289
290
291
292
    return m;
  }

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

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

302
303
304
    return m;
  }

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

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

316
    const DOFAdmin *admin = NULL;
317
    const char *format;
318

319
320
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
321

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

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

    return result;
  }
371
372
373
374
375
376
377
378
379

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

380
    for (int i = 0; i < numberOfBasFcts; i++)
381
382
383
384
385
386
387
388
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

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

391
392
393
    TEST_EXIT_DBG(fct)("no function to interpolate\n");

    interFct = fct;
394

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

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

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

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
423
424
425
426
427
    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);
    }
428
429
430
  }

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

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

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

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

458
459
    quad_fast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
460
461
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
462
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, Int_fct);
463
464
465
466
467
468
469

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::Int_fct(ElInfo *elinfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
472
473
    double det = elinfo->getDet();
    const T* uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int numPoints = quad_fast->getNumPoints();
    double normT = 0.0;
474

Thomas Witkowski's avatar
Thomas Witkowski committed
475
476
    for (int iq = 0; iq < numPoints; iq++)
      normT += quad_fast->getWeight(iq) * (uh_vec[iq]);
477

Thomas Witkowski's avatar
Thomas Witkowski committed
478
    norm += det * normT;
479
480
481
482
483
484
485

    return 0;
  }

  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
486
    FUNCNAME("DOFVector::L1Norm()");
487
  
Thomas Witkowski's avatar
Thomas Witkowski committed
488
    if (!q) {
489
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
490
491
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
492

493
494
    quad_fast = FastQuadrature::provideFastQuadrature(this-> feSpace->getBasisFcts(),
						      *q, INIT_PHI);
495
496
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
497
    Mesh* mesh = this->feSpace->getMesh();
498

499
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, 
500
501
502
503
504
505
506
507
		   L1Norm_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L1Norm_fct(ElInfo *elinfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
508
509
510
511
    double det = elinfo->getDet();
    const T* uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int numPoints = quad_fast->getNumPoints();
    double normT = 0.0;
512

Thomas Witkowski's avatar
Thomas Witkowski committed
513
514
    for (int iq = 0; iq < numPoints; iq++)
      normT += quad_fast->getWeight(iq) * abs(uh_vec[iq]);
515

Thomas Witkowski's avatar
Thomas Witkowski committed
516
    norm += det * normT;
517
518
519
520
521
522
523
524

    return 0;
  }


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

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
    quad_fast = FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), 
						      *q, INIT_PHI);
534
535
536

    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
537
    Mesh* mesh = this->feSpace->getMesh();
Thomas Witkowski's avatar
Thomas Witkowski committed
538
539
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, 
		   L2Norm_fct);
540
541
542
543
544
545
546

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L2Norm_fct(ElInfo *elinfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
547
548
549
550
    double det = elinfo->getDet();
    const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int numPoints = quad_fast->getNumPoints();
    double normT = 0.0;
551

Thomas Witkowski's avatar
Thomas Witkowski committed
552
553
    for (int iq = 0; iq < numPoints; iq++)
      normT += quad_fast->getWeight(iq) * sqr(uh_vec[iq]);
554

Thomas Witkowski's avatar
Thomas Witkowski committed
555
    norm += det * normT;
556
557
558
559
560
561
562
563

    return 0;
  }

  template<typename T>
  int DOFVector<T>::H1Norm_fct(ElInfo *elinfo)
  {
    double det = elinfo->getDet();
Thomas Witkowski's avatar
Thomas Witkowski committed
564
565
    const WorldVector<T> *grduh_vec = 
      traverseVector->getGrdAtQPs(elinfo, NULL, quad_fast, NULL);
566
567
    int dimOfWorld = Global::getGeo(WORLD);
    int numPoints = quad_fast->getNumPoints();
Thomas Witkowski's avatar
Thomas Witkowski committed
568
    double normT = 0.0;
569

Thomas Witkowski's avatar
Thomas Witkowski committed
570
571
572
573
574
575
576
577
    for (int iq = 0; iq < numPoints; iq++) {
      double norm2 = 0.0;
      for (int j = 0; j < dimOfWorld; j++)
	norm2 += sqr(grduh_vec[iq][j]);
      
      normT += quad_fast->getWeight(iq) * norm2;
    }
    norm += det * normT; 
578
579
580
581
582
583
584
585

    return 0;
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
588
    if (!q) {
589
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
590
591
592
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

593
    quad_fast = FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
594
						      *q, INIT_GRD_PHI);
595
596

    norm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
597
    traverseVector = const_cast<DOFVector<T>*>(this); 
598
    Mesh *mesh = this->feSpace->getMesh();
599

600
601
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | 
		   Mesh::FILL_GRD_LAMBDA, H1Norm_fct);
602
603
604
605
606
607

    return norm;
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
608
					std::vector<DegreeOfFreedom> &newDOF)
609
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
610
611
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
612
	vec[newDOF[i]] = vec[i];
613
614
615
616
617
618
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
619
620
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
621
	 op != this->operators.end(); ++op)
622
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
623

624
625
626
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
627
628
629
630
631
  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
632
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
633
634
635
      (*it)->finishAssembling();
  }

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

651
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
652
    } else {
653
      this->boundaryManager = NULL;
654
655
    }

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
681
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
682
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
683
684
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
685
686
      ("no admin or different admins: %8X, %8X\n",
       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
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
704
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
705
706
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
707
708
      ("no admin or different admins: %8X, %8X\n",
       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
725
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
726
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
727
728
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
729
730
      ("no admin or different admins: %8X, %8X\n",
       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
    TEST_EXIT(x.getFESpace() && y.getFESpace())
749
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
750
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) && (admin == y.getFESpace()->getAdmin()))
751
752
      ("no admin or different admins: %8X, %8X\n",
       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()");

Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
    TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && 
	      x.getFESpace() && result.getFESpace())
777
778
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
779
780
781
782
783
784
785
786
787

    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()))))
788
789
790
791
792
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());


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
    TEST_EXIT(x.getFESpace() && y.getFESpace())
807
808
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

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

    TEST_EXIT((admin) && (admin == y.getFESpace()->getAdmin()))
812
813
      ("no admin or different admins: %8X, %8X\n",
       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
    TEST_EXIT(x.getFESpace() && y.getFESpace())
873
874
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

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

    TEST_EXIT(admin && (admin == y.getFESpace()->getAdmin()))
878
879
      ("no admin or different admins: %8X, %8X\n",
       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
       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
904
905
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
906
	  y[i] = alpha * y[i] + x[i];
907
908
909
910
911
912
913
914
915
  }

  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
916
917
918
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = scal * (*vIterator);
919
920
921
922
923
924
925
926
927
928
929

    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
930
931
932
933
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = (*vIterator) + scal;

934
935
936
937
938
939
940
941
942
943
944
    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
945
946
947
    for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	 !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator) 
      *rIterator = (*v1Iterator) + (*v2Iterator);
948

Thomas Witkowski's avatar
Thomas Witkowski committed
949
    return result;
950
951
952
  }

  template<typename T>
953
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
954
  {
955
    static T* localVec = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
956
    static int localVecSize = 0;
957
958
959
    const DOFAdmin* admin = feSpace->getAdmin();
    T *result;
  
960
    if (d) {