DOFVector.hh 31.4 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();
Thomas Witkowski's avatar
Thomas Witkowski committed
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);
442
    int number = basFct->getNumber();
443
    for (int i = 0; i < number; i++)
444
445
446
447
448
449
450
451
      (*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
452
    FUNCNAME("DOFVector::Int()");
453
  
454
    Mesh* mesh = this->feSpace->getMesh();
455

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

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

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

    return norm;  
  }

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

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

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

    return 0;
  }

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

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

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

    return norm;  
  }

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

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

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

    return 0;
  }


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

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

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

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

    return norm;  
  }

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

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

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

    return 0;
  }

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

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

    return 0;
  }


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

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

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

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

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

    return norm;
  }

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

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

628
629
630
    return fillFlag;
  }

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

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

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

660
661
662
663
664
665
    return *this;
  }

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

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

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

676
677
678
679
680
681
682
    return x;
  }


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

699
700
701
702
703
704
    return x;
  }

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

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

721
722
723
724
725
726
    return x;
  }

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

    return x;
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
759
    T dot = 0;
760
761
762

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

Thomas Witkowski's avatar
Thomas Witkowski committed
767
    return dot;
768
769
770
  }

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

777
778
    FUNCNAME("DOFVector<T>::mv()");

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

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


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

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

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

813
814
815
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

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

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

879
880
881
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
953
    return result;
954
955
956
  }

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

963
964
    T *result;
  
965
    if (d) {
966
967
      result = d;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
968
969
970
971
972
973
974
975
      if (localVec && nBasFcts > localVecSize)  {
	delete [] localVec;
	localVec = new T[nBasFcts]; 
      }
      if (!localVec)
	localVec = new T[nBasFcts]; 

      localVecSize = nBasFcts;
976
977
      result = localVec;
    }
978

979
980
    const DegreeOfFreedom *localIndices = 
      feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);
981

982
    for (int i = 0; i < nBasFcts; i++)
983
      result[i] = (*this)[localIndices[i]];
984
985

    return result;
986
987
988
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
989
990
  const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo, 
					 const Quadrature *quad,
991
					 const FastQuadrature *quadFast,
Thomas Witkowski's avatar
Thomas Witkowski committed
992
					 T *vecAtQPs) const
993
994
995
  {
    FUNCNAME("DOFVector<T>::getVecAtQPs()");
  
996
    TEST_EXIT(quad || quadFast)("neither quad nor quadFast defined\n");
997

Thomas Witkowski's avatar
Thomas Witkowski committed
998
    if (quad && quadFast)
999
1000
      TEST_EXIT(quad == quadFast->getQuadrature())
	("quad != quadFast->quadrature\n");
1001

1002
    TEST_EXIT(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
1003
1004
      ("invalid basis functions");

1005
1006
    Element *el = elInfo->getElement();
    const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
1007
    const BasisFunction *basFcts = feSpace->getBasisFcts();
1008
1009
    int numPoints = quadrature->getNumPoints();
    int nBasFcts  = basFcts->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
1010
    int j;
1011
1012
1013
1014
1015
1016
    static T *localvec = NULL;
    T *result;

    if (vecAtQPs) {
      result = vecAtQPs;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1017
1018
      if (localvec) 
	delete [] localvec;
1019
      localvec = new T[numPoints];
Thomas Witkowski's avatar
Thomas Witkowski committed
1020
      for (int i = 0; i < numPoints; i++)
1021
	localvec[i] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed