DOFVector.hh 31.3 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
113
114
115
  {
    init(f, n);
  } 

  template<typename T>
116
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
117
118
119
  {
    name = n;
    feSpace = f;
Thomas Witkowski's avatar
Thomas Witkowski committed
120
    if (feSpace && feSpace->getAdmin())
121
      (feSpace->getAdmin())->addDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
122

123
    this->boundaryManager = new BoundaryManager(f);
124
125
126
127
128
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
129
    if (feSpace && feSpace->getAdmin())
130
131
      (feSpace->getAdmin())->removeDOFIndexed(this);

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
198
    checkFeSpace(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

Thomas Witkowski's avatar
Thomas Witkowski committed
214
    checkFeSpace(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

Thomas Witkowski's avatar
Thomas Witkowski committed
230
    checkFeSpace(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

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

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
260

    checkFeSpace(feSpace, vec);

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


  template<typename T>
  T DOFVector<T>::min() const
  {
277
278
    FUNCNAME("DOFVector<T>::min()");
    
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    checkFeSpace(feSpace, vec);
280
281

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

286
287
288
289
290
291
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
292
293
    FUNCNAME("DOFVector<T>::max()");
    
Thomas Witkowski's avatar
Thomas Witkowski committed
294
    checkFeSpace(feSpace, vec);
295

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

301
302
303
    return m;
  }

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

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

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

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

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

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

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

  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;

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

    return val;
  }

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

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

    interFct = fct;
393

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

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

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

    traverseVector = this;

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

  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
434
435
436
437
438
439
440
    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);
441
    int number = basFct->getNumber();
442
    for (int i = 0; i < number; i++)
443
444
445
446
447
448
449
450
      (*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
451
    FUNCNAME("DOFVector::Int()");
452
453
454
  
    Mesh* mesh = feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
455
456
457
458
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475

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

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   Int_fct);

    return norm;  
  }

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

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

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

    return 0;
  }

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

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),*q,INIT_PHI);
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
502
    Mesh* mesh = feSpace->getMesh();
503
504
505
506
507
508
509
510
511
512
513
514
515

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   L1Norm_fct);

    return norm;  
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
524
    norm += det * normT;
525
526
527
528
529
530
531
532

    return 0;
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
535
536
537
538
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
539
540
541
542
543
544
545

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, 
						      INIT_PHI);

    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
546
    Mesh* mesh = feSpace->getMesh();
Thomas Witkowski's avatar
Thomas Witkowski committed
547
548
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, 
		   L2Norm_fct);
549
550
551
552
553
554
555

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L2Norm_fct(ElInfo *elinfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
556
557
558
559
    double det = elinfo->getDet();
    const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int numPoints = quad_fast->getNumPoints();
    double normT = 0.0;
560

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

Thomas Witkowski's avatar
Thomas Witkowski committed
564
    norm += det * normT;
565
566
567
568
569
570
571
572

    return 0;
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
579
580
581
582
583
584
585
586
    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; 
587
588
589
590
591
592
593
594

    return 0;
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
597
598
599
600
601
602
603
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree() - 2;
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, INIT_GRD_PHI);
604
605

    norm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
606
607
    traverseVector = const_cast<DOFVector<T>*>(this); 
    Mesh *mesh = feSpace->getMesh();
608
609
610
611
612
613
614
615
616
617

    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
		   Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA,
		   H1Norm_fct);

    return norm;
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
618
					std::vector<DegreeOfFreedom> &newDOF)
619
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
620
621
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
622
	vec[newDOF[i]] = vec[i];
623
624
625
626
627
628
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
631
	 op != this->operators.end(); ++op)
632
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
633

634
635
636
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
637
638
639
640
641
  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
642
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
645
      (*it)->finishAssembling();
  }

646
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
647
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
648
  {
649
    feSpace = rhs.feSpace;
650
    DOFVectorBase<T>::feSpace = rhs.feSpace;
651
    this->nBasFcts = rhs.nBasFcts;
652
    vec = rhs.vec;
653
    this->elementVector.change_dim(this->nBasFcts);
654
655
656
657
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
658
    if (rhs.boundaryManager) {
659
660
661
      if (this->boundaryManager) 
	delete this->boundaryManager; 

662
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
663
    } else {
664
      this->boundaryManager = NULL;
665
666
    }

667
668
669
670
671
672
    return *this;
  }

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

675
676
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
677

Thomas Witkowski's avatar
Thomas Witkowski committed
678
679
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
680
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
681
      (*vecIterator) *= scal; 
682

683
684
685
686
687
688
689
    return x;
  }


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

706
707
708
709
710
711
    return x;
  }

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

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

728
729
730
731
732
733
    return x;
  }

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

    return x;
  }

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

759
    TEST_EXIT(x.getFESpace() && y.getFESpace())
760
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
761
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) && (admin == y.getFESpace()->getAdmin()))
762
763
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
764
765
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
766
    T dot = 0;
767
768
769

    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
770
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
771
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
772
      dot += (*xIterator) * (*yIterator);
773

Thomas Witkowski's avatar
Thomas Witkowski committed
774
    return dot;
775
776
777
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
778
779
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
780
  {
781
782
783
    // Unfortunately needed
    using namespace mtl;

784
785
    FUNCNAME("DOFVector<T>::mv()");

Thomas Witkowski's avatar
Thomas Witkowski committed
786
787
    TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && 
	      x.getFESpace() && result.getFESpace())
788
789
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
790
791
792
793
794
795
796
797
798

    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()))))
799
800
801
802
803
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());


804
805
806
807
808
    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 
809
810
811
812
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
813
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
814
815
816
  {
    FUNCNAME("DOFVector<T>::axpy()");

817
    TEST_EXIT(x.getFESpace() && y.getFESpace())
818
819
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

820
821
822
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

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

883
    TEST_EXIT(x.getFESpace() && y.getFESpace())
884
885
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

886
887
888
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

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

    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
941
942
943
944
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = (*vIterator) + scal;

945
946
947
948
949
950
951
952
953
954
955
    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
956
957
958
    for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	 !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator) 
      *rIterator = (*v1Iterator) + (*v2Iterator);
959

Thomas Witkowski's avatar
Thomas Witkowski committed
960
    return result;
961
962
963
  }

  template<typename T>
964
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
965
  {
966
    static T* localVec = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
967
    static int localVecSize = 0;
968
    const DOFAdmin* admin = feSpace->getAdmin();
969

970
971
    T *result;
  
972
    if (d) {
973
974
      result = d;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
975
976
977
978
979
980
981
982
      if (localVec && nBasFcts > localVecSize)  {
	delete [] localVec;
	localVec = new T[nBasFcts]; 
      }
      if (!localVec)
	localVec = new T[nBasFcts]; 

      localVecSize = nBasFcts;
983
984
      result = localVec;
    }
985

986
987
    const DegreeOfFreedom *localIndices = 
      feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);
988

989
    for (int i = 0; i < nBasFcts; i++)
990
      result[i] = (*this)[localIndices[i]];
991
992

    return result;
993
994
995
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
996
997
  const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo, 
					 const Quadrature *quad,
998
					 const FastQuadrature *quadFast,
Thomas Witkowski's avatar
Thomas Witkowski committed
999
					 T *vecAtQPs) const
1000
1001
1002
  {
    FUNCNAME("DOFVector<T>::getVecAtQPs()");
  
1003
    TEST_EXIT(quad || quadFast)("neither quad nor quadFast defined\n");
1004

Thomas Witkowski's avatar
Thomas Witkowski committed
1005
    if (quad && quadFast)
1006
1007
      TEST_EXIT(quad == quadFast->getQuadrature())
	("quad != quadFast->quadrature\n");
1008

1009
    TEST_EXIT(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
1010
1011
      ("invalid basis functions");

1012
1013
    Element *el = elInfo->getElement();
    const Quadrature *quadrature = quadFast ? quadFast->getQuadrature