DOFVector.hh 31.6 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
431
432
433
434
  }

  template<typename T>
  int DOFVector<T>::interpolFct(ElInfo* elinfo)
  {
    const BasisFunction *basFct = traverseVector->getFESpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFESpace()->getAdmin();
435
436
437
438
439
440
441
    const DegreeOfFreedom *dof = 
      basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), admin, NULL);
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
    int nBasFcts = basFct->getNumber();
    for (int i = 0; i < nBasFcts; i++)
442
443
444
445
446
447
448
449
      (*traverseVector)[dof[i]] = inter_val[i];

    return 0;
  }

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

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

459
460
    quad_fast = FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), 
						      *q, INIT_PHI);
461
462
463
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);

464
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, Int_fct);
465
466
467
468
469
470
471

    return norm;  
  }

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

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

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

    return 0;
  }

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

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

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

    return norm;  
  }

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

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

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

    return 0;
  }


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

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

534
535
    quad_fast = FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), 
						      *q, INIT_PHI);
536
537
538

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

    return norm;  
  }

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

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

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

    return 0;
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
572
573
574
575
576
577
578
579
    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; 
580
581
582
583
584
585
586
587

    return 0;
  }


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

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

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

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

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

    return norm;
  }

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

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

626
627
628
    return fillFlag;
  }

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

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

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

658
659
660
661
662
663
    return *this;
  }

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

666
667
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
668

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

674
675
676
677
678
679
680
    return x;
  }


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

697
698
699
700
701
702
    return x;
  }

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

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

719
720
721
722
723
724
    return x;
  }

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

    return x;
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
757
    T dot = 0;
758
759
760

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

Thomas Witkowski's avatar
Thomas Witkowski committed
765
    return dot;
766
767
768
  }

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

775
776
    FUNCNAME("DOFVector<T>::mv()");

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

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


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

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

808
    TEST_EXIT(x.getFESpace() && y.getFESpace())
809
810
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

811
812
813
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT((admin) && (admin == y.getFESpace()->getAdmin()))
814
815
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
816
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
817
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
818
819
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
820
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
       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
840
841
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
842
	  y[i] += alpha * x[i];
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
872
873
  }

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

874
    TEST_EXIT(x.getFESpace() && y.getFESpace())
875
876
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

877
878
879
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
951
    return result;
952
953
954
  }

  template<typename T>
955
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
956
  {
957
    static T* localVec = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
958
    static int localVecSize = 0;
959
960
961
    const DOFAdmin* admin = feSpace->getAdmin();
    T *result;
  
962
    if (d) {
963
964
      result = d;
    } else {
965
966
967
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
968
969
970
971
972
973
974
975
      if (localVec && nBasFcts > localVecSize)  {
	delete []