DOFVector.hh 32.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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
    createTempData();
  }
  

  template<typename T>
  DOFVectorBase<T>::~DOFVectorBase()
  {
    for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
      delete [] localIndices[i];
      delete [] localVectors[i];
      delete grdPhis[i];
      delete grdTmp[i];
      delete D2Phis[i];
    }
  }


  template<typename T>
  void DOFVectorBase<T>::createTempData()
  {
100
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
101
    localVectors.resize(omp_get_overall_max_threads());
102
    grdPhis.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
103
    grdTmp.resize(omp_get_overall_max_threads());
104
    D2Phis.resize(omp_get_overall_max_threads());
105

106
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
107
108
      localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
      localVectors[i] = new T[this->nBasFcts];
109
110
111
      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);
112
    }    
113
  }
114

Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
  
  template<typename T>
117
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
118
    : DOFVectorBase<T>(f, n),
119
      coarsenOperation(NO_OPERATION)
120
  {
121
    vec.resize(0);
122
123
124
    init(f, n);
  } 

125

126
  template<typename T>
127
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
128
  {
129
130
131
132
133
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
    this->boundaryManager = new BoundaryManager(f);
134
135
136
137
138
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
139
140
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
141

Thomas Witkowski's avatar
Thomas Witkowski committed
142
    if (this->boundaryManager)
143
      delete this->boundaryManager;
144
145

    vec.clear();
146
147
148
149
150
  }

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

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

160
161
162
    std::vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
163

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

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

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

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

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

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

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
214
    FUNCNAME("DOFVector<T>::asum()");
215

216
    checkFeSpace(this->feSpace, vec);
217

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

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
230
    FUNCNAME("DOFVector<T>::sum()");
231

232
    checkFeSpace(this->feSpace, vec);
233

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

    return(nrm);
  }

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

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

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


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

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

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


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

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

289
290
291
    return m;
  }

292

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

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

305
    return m;
306

307
308
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
309
310
311
312
313
314
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334

  template<typename T>
  T DOFVector<T>::average() const
  {
    FUNCNAME("DOFVector<T>::average()");
    
    checkFeSpace(this->feSpace, vec);

    int count = 0;
    T m = 0;
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
      m += *vecIterator;
      count++;
    }

    return m / count;
  }


335
336
337
  template<typename T>
  void DOFVector<T>::print() const
  {
338
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
339

340
    const DOFAdmin *admin = NULL;
341
    const char *format;
342

343
344
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
345

346
    MSG("Vec `%s':\n", this->name.c_str());
347
    int j = 0;
348
349
350
351
352
353
354
355
356
    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
357
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
358
	if ((j % 3) == 0) {
359
360
	  if (j) 
	    Msg::print("\n");
361
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
362
	} else {
363
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
364
	}
365
	j++;
366
      }
367
      Msg::print("\n");
368
    } else {
369
370
	MSG("no DOFAdmin, print whole vector.\n");
    
371
372
373
374
375
376
	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 {
377
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
378
	  }
379
380
381
382
383
384
385
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }

386
387
388
389
390
391
392
393
394
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
395
396
397
398
399

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

404
    for (int i = 0; i < nBasisFcts; i++)
405
406
407
408
409
410
411
412
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
413
    FUNCNAME("DOFVector::interpol()");
414

Thomas Witkowski's avatar
Thomas Witkowski committed
415
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
416
417

    interFct = fct;
418

419
    if (!this->getFeSpace()) {
420
421
422
423
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
424

425
    if (!(this->getFeSpace()->getAdmin())) {
426
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
427
	  this->getFeSpace()->getName().c_str());
428
429
      return;
    }
430
  
431
    if (!(this->getFeSpace()->getBasisFcts())) {
432
433
434
435
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
436

437
438
439
440
441
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
442
443
444

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
445
    TraverseStack stack;
446
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
449
450
451
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
452
453
454
  }

  template<typename T>
455
  void DOFVector<T>::interpolFct(ElInfo* elinfo)
456
  {
457
458
    const BasisFunction *basFct = traverseVector->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFeSpace()->getAdmin();
459
460
461
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
462
463
464

    DegreeOfFreedom *myLocalIndices = this->localIndices[omp_get_thread_num()];
    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), admin, myLocalIndices);
465
466
    int nBasFcts = basFct->getNumber();
    for (int i = 0; i < nBasFcts; i++)
467
      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
468
469
470
471
472
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
477
    if (!q) {
478
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
479
480
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
481

482
    FastQuadrature *quadFast = 
483
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
484

485
486
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
Thomas Witkowski's avatar
Thomas Witkowski committed
487
    std::vector<T> uh_vec(nPoints);
488
489
490
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
491
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
492
493
494
495
496
497
498
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getVecAtQPs(elInfo, NULL, quadFast, &(uh_vec[0]));
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;
499

500
501
      elInfo = stack.traverseNext(elInfo);
    }
502

503
    return result;  
504
505
506
507
508
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
513
    if (!q) {
514
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
515
516
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
517

518
519
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
520

521
522
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
523
    std::vector<T> uh_vec(nPoints);
524
525
526
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
527
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
528
529
530
531
532
533
534
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getVecAtQPs(elInfo, NULL, quadFast, &(uh_vec[0]));
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
535

536
537
      elInfo = stack.traverseNext(elInfo);
    }
538

539
    return result;  
540
541
542
543
544
  }

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

547
548
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
549
    if (!q) {
550
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
553

554
555
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
556

557
558
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
559
    std::vector<T> uh_vec(nPoints);
560
561
562
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
563
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
564
565
566
567
568
569
570
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getVecAtQPs(elInfo, NULL, quadFast, &(uh_vec[0]));
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;
571

572
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
573
    }
574

575
    return result;  
576
577
  }

578
  template<typename T>  
579
580
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
581
    FUNCNAME("DOFVector::H1NormSquare()");
582

583
584
    Mesh* mesh = this->feSpace->getMesh();

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

590
591
592
593
594
595
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
596
    std::vector<WorldVector<T> > grduh_vec(nPoints);
597
598
599
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
600
601
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
602
603
604
605
606
607
608
609
610
611
612
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getGrdAtQPs(elInfo, NULL, quadFast, &(grduh_vec[0]));

      for (int iq = 0; iq < nPoints; iq++) {
	double norm2 = 0.0;
	for (int j = 0; j < dimOfWorld; j++)
	  norm2 += sqr(grduh_vec[iq][j]);
	normT += quadFast->getWeight(iq) * norm2;
      }
613

614
      result += det * normT;
615

616
617
      elInfo = stack.traverseNext(elInfo);
    }
618

619
    return result;  
620
621
622
623
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
624
					std::vector<DegreeOfFreedom> &newDOF)
625
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
626
627
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
628
	vec[newDOF[i]] = vec[i];
629
630
631
632
633
634
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
637
	 op != this->operators.end(); ++op)
638
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
639

640
641
642
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
645
646
647
  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
648
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
649
650
651
      (*it)->finishAssembling();
  }

652
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
653
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
654
  {
655
    this->feSpace = rhs.feSpace;
656
    this->dim = rhs.dim;
657
    this->nBasFcts = rhs.nBasFcts;
658
    vec = rhs.vec;
659
    this->elementVector.change_dim(this->nBasFcts);
660
661
662
663
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
664
    
665
    if (rhs.boundaryManager) {
666
667
668
      if (this->boundaryManager) 
	delete this->boundaryManager; 

669
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
670
    } else {
671
      this->boundaryManager = NULL;
672
673
    }

674
675
676
677
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

678
679
680
681
682
683
    return *this;
  }

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

686
687
    TEST_EXIT_DBG(x.getFeSpace() && x.getFeSpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFeSpace(), x.getFeSpace()->getAdmin());
688

Thomas Witkowski's avatar
Thomas Witkowski committed
689
690
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
691
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
692
      (*vecIterator) *= scal; 
693

694
695
696
697
698
699
700
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
701
702
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
703
704
705
706
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
707
      ("no admin or different admins: %8X, %8X\n",
708
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
709
710
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
711
712
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
713
714
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
715
716
      *xIterator += *yIterator; 

717
718
719
720
721
722
    return x;
  }

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

725
726
727
728
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
729
      ("no admin or different admins: %8X, %8X\n",
730
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
731
732
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
733
734
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
737
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
738

739
740
741
742
743
744
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
745
746
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
747
748
749
750
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
751
      ("no admin or different admins: %8X, %8X\n",
752
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
753
754
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
755
756
    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
757
758
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
759
      *xIterator *= *yIterator; 
760
761
762
763
764
765
766
767

    return x;
  }

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

770
771
772
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT((admin = x.getFeSpace()->getAdmin()) && (admin == y.getFeSpace()->getAdmin()))
773
      ("no admin or different admins: %8X, %8X\n",
774
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
775
776
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
777
    T dot = 0;
778
779
780

    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
781
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
782
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
783
      dot += (*xIterator) * (*yIterator);
784

Thomas Witkowski's avatar
Thomas Witkowski committed
785
    return dot;
786
787
788
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
789
790
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
791
  {
792
793
794
    // Unfortunately needed
    using namespace mtl;

795
796
    FUNCNAME("DOFVector<T>::mv()");

797
798
799
800
    TEST_EXIT(a.getRowFeSpace() && a.getColFeSpace() && 
	      x.getFeSpace() && result.getFeSpace())
      ("getFeSpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFeSpace(), a.getColFeSpace(), x.getFeSpace(), result.getFeSpace());
801

802
803
    const DOFAdmin *rowAdmin = a.getRowFeSpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFeSpace()->getAdmin();
804
805

    TEST_EXIT((rowAdmin && colAdmin) && 
806
807
808
809
	      (((transpose == NoTranspose) && (colAdmin == x.getFeSpace()->getAdmin()) && 
		(rowAdmin == result.getFeSpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFeSpace()->getAdmin()) && 
		(colAdmin == result.getFeSpace()->getAdmin()))))
810
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
811
812
       a.getRowFeSpace()->getAdmin(), a.getColFeSpace()->getAdmin(), 
       x.getFeSpace()->getAdmin(), result.getFeSpace()->getAdmin());
813
814


815
816
817
818
819
    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 
820
821
822
823
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
824
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
825
826
827
  {
    FUNCNAME("DOFVector<T>::axpy()");

828
829
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
830

831
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
832

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

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

894
895
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
896

897
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
898

899
    TEST_EXIT(admin && (admin == y.getFeSpace()->getAdmin()))
900
      ("no admin or different admins: %8X, %8X\n",
901
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
902
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
903
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
904
905
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
906
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
       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
926
927
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
928
	  y[i] = alpha * y[i] + x[i];
929
930
931
932
933
934
935
936
937
  }

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

    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
952
953
954
955
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = (*vIterator) + scal;

956
957
958
    return result;
  }

959

960
  template<typename T>
961
  inline const DOFVector<T>& add(const DOFVector<T>& v1, 
962
963
964
965
966
967
				 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
968
969
970
    for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	 !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator) 
      *rIterator = (*v1Iterator) + (*v2Iterator);
971

Thomas Witkowski's avatar
Thomas Witkowski committed
972
    return result;
973
974
  }

975

976
  template<typename T>
977
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const