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

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

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

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

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

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

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

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


} // namespace mtl



68
69
70
namespace AMDiS {

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

80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
    createTempData();
  }
  

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


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

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

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

124

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

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

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

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

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

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

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

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

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

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

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

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

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

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

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

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

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

    return(nrm);
  }

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

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

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

    return(nrm);
  }

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

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

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


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

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

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


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

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

288
289
290
    return m;
  }

291

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

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

304
    return m;
305

306
307
  }

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

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

  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;
  }


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

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

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

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

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

    return result;
  }
394
395
396
397
398

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

403
    for (int i = 0; i < nBasisFcts; i++)
404
405
406
407
408
409
410
411
      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
412
    FUNCNAME("DOFVector::interpol()");
413

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

    interFct = fct;
417

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

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

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

    traverseVector = this;

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

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

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

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

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

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

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

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

502
    return result;  
503
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
    mtl::dense_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
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
531
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
532
533
534
      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
545
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
546
    FUNCNAME("DOFVector::L2NormSquare()");
547

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

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

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

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

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

576
    return result;  
577
578
  }

579

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

585
586
    Mesh* mesh = this->feSpace->getMesh();

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

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

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
598
    std::vector<WorldVector<T> > grduh_vec(nPoints);
599
600
601
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
602
603
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
604
605
606
607
608
609
610
611
612
613
614
    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;
      }
615

616
      result += det * normT;
617

618
619
      elInfo = stack.traverseNext(elInfo);
    }
620

621
    return result;  
622
623
624
625
  }

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

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

642
643
644
    return fillFlag;
  }

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

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

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

676
677
678
679
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

680
681
682
683
684
685
    return *this;
  }

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

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

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

696
697
698
699
700
701
702
    return x;
  }


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

719
720
721
722
723
724
    return x;
  }

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

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

741
742
743
744
745
746
    return x;
  }

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

    return x;
  }

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

772
773
774
    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()))
775
      ("no admin or different admins: %8X, %8X\n",
776
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
777
778
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
779
    T dot = 0;
780
781
782

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

Thomas Witkowski's avatar
Thomas Witkowski committed
787
    return dot;
788
789
790
  }

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

797
798
    FUNCNAME("DOFVector<T>::mv()");

799
800
801
802
    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());
803

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

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


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

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

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

833
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
834

835
    TEST_EXIT((admin) && (admin == y.getFeSpace()->getAdmin()))
836
      ("no admin or different admins: %8X, %8X\n",
837
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
838
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
839
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
840
841
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
842
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
       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
862
      for (i = 0; i < maxI; i++)
863
	if (!admin->isDofFree(i))
864
	  y[i] += alpha * x[i];
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