DOFVector.hh 31.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include <list>
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include <algorithm>
15
#include <math.h>
16

17
18
19
20
21
#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>

22
23
24
25
26
27
28
29
30
31
32
33
#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"
34
#include "OpenMP.h"
35
#include "Operator.h"
36
#include "Parameters.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
37
#include "Traverse.h"
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
68
69
70
71
72
73
74
75
76
77
78
79
// 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



80
81
82
namespace AMDiS {

  template<typename T>
83
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
84
85
    : feSpace(f),
      name(n),
86
      elementVector(f->getBasisFcts()->getNumber()),
Thomas Witkowski's avatar
Thomas Witkowski committed
87
      boundaryManager(NULL)
88
  {    
Thomas Witkowski's avatar
Thomas Witkowski committed
89
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
90
    dim = feSpace->getMesh()->getDim();
91

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    createTempData();
  }
  

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


  template<typename T>
  void DOFVectorBase<T>::createTempData()
  {
109
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
110
    localVectors.resize(omp_get_overall_max_threads());
111
    D2Phis.resize(omp_get_overall_max_threads());
112

113
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
114
      localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
115
      localVectors[i].change_dim(this->nBasFcts);
116
      D2Phis[i] = new DimMat<double>(dim, NO_INIT);
117
    }    
118
  }
119

Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
  
  template<typename T>
122
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
123
    : DOFVectorBase<T>(f, n),
124
      coarsenOperation(NO_OPERATION)
125
  {
126
    vec.resize(0);
127
128
129
    init(f, n);
  } 

130

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

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
144
145
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
146

Thomas Witkowski's avatar
Thomas Witkowski committed
147
    if (this->boundaryManager)
148
      delete this->boundaryManager;
149
150

    vec.clear();
151
152
153
154
155
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
156
157
158
159
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
160
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
161
162
163
164
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

165
166
167
    std::vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
168

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

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

176
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
177
	  (*this)[irow] += factor * elVec[i];
178
179
	else  
	  (*this)[irow] = factor * elVec[i];	
Thomas Witkowski's avatar
Thomas Witkowski committed
180
181
182
183
      }
    }
  }

184
185
186
187
188
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

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

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

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
219
    FUNCNAME("DOFVector<T>::asum()");
220

221
    checkFeSpace(this->feSpace, vec);
222

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

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
235
    FUNCNAME("DOFVector<T>::sum()");
236

237
    checkFeSpace(this->feSpace, vec);
238

239
    double nrm = 0.0;    
Thomas Witkowski's avatar
Thomas Witkowski committed
240
241
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
242
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
243
244
245
246
247
248
249
250
      nrm += *vecIterator;

    return(nrm);
  }

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

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

255
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
256
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
257
258
259
260
261
262
263
      *vecIterator = alpha ; 
  }


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

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

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


  template<typename T>
  T DOFVector<T>::min() const
  {
285
286
    FUNCNAME("DOFVector<T>::min()");
    
287
    checkFeSpace(this->feSpace, vec);
288
289

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

294
295
296
    return m;
  }

297

298
299
300
  template<typename T>
  T DOFVector<T>::max() const 
  {
301
302
    FUNCNAME("DOFVector<T>::max()");
    
303
    checkFeSpace(this->feSpace, vec);
304

305
    T m;    
306
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
307
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
308
      m = std::max(m, *vecIterator);
309

310
    return m;
311

312
313
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
316
317
318
319
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339

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


340
341
342
  template<typename T>
  void DOFVector<T>::print() const
  {
343
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
344

345
    const DOFAdmin *admin = NULL;
346
    const char *format;
347

348
349
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
350

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

391
392
393
394
395
396
397
398
399
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
400
401
402
403
404

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

409
    for (int i = 0; i < nBasisFcts; i++)
410
411
412
413
414
415
416
417
      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
418
    FUNCNAME("DOFVector::interpol()");
419

Thomas Witkowski's avatar
Thomas Witkowski committed
420
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
421
422

    interFct = fct;
423

424
    if (!this->getFeSpace()) {
425
426
427
428
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
429

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

442
443
444
445
446
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
447
448
449

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
450
    TraverseStack stack;
451
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
454
455
456
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
457
458
459
  }

  template<typename T>
460
  void DOFVector<T>::interpolFct(ElInfo* elinfo)
461
  {
462
463
    const BasisFunction *basFct = traverseVector->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFeSpace()->getAdmin();
464
465
466
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
467
468
469

    DegreeOfFreedom *myLocalIndices = this->localIndices[omp_get_thread_num()];
    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), admin, myLocalIndices);
470
471
    int nBasFcts = basFct->getNumber();
    for (int i = 0; i < nBasFcts; i++)
472
      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
473
474
475
476
477
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
482
    if (!q) {
483
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
484
485
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
486

487
    FastQuadrature *quadFast = 
488
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
489

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

505
506
      elInfo = stack.traverseNext(elInfo);
    }
507

508
    return result;  
509
510
  }

511

512
513
514
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
515
    FUNCNAME("DOFVector::L1Norm()");
516
  
517
518
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
519
    if (!q) {
520
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
523

524
525
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
526

527
528
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
529
    mtl::dense_vector<T> uh_vec(nPoints);
530
531
532
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
533
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
534
535
536
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
537
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
538
539
540
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
541

542
543
      elInfo = stack.traverseNext(elInfo);
    }
544

545
    return result;  
546
547
  }

548

549
550
551
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
552
    FUNCNAME("DOFVector::L2NormSquare()");
553

554
555
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
556
    if (!q) {
557
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
558
559
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
560

561
562
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
563

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

579
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
580
    }
581

582
    return result;  
583
584
  }

585

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

591
592
    Mesh* mesh = this->feSpace->getMesh();

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

598
599
600
601
602
603
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

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

622
      result += det * normT;
623

624
625
      elInfo = stack.traverseNext(elInfo);
    }
626

627
    return result;  
628
629
630
631
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
632
					std::vector<DegreeOfFreedom> &newDOF)
633
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
634
635
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
636
	vec[newDOF[i]] = vec[i];
637
638
639
640
641
642
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
645
	 op != this->operators.end(); ++op)
646
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
647

648
649
650
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
651
652
653
654
655
  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
656
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
659
      (*it)->finishAssembling();
  }

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

677
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
678
    } else {
679
      this->boundaryManager = NULL;
680
681
    }

682
683
684
685
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

686
687
688
689
690
691
    return *this;
  }

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

694
695
    TEST_EXIT_DBG(x.getFeSpace() && x.getFeSpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFeSpace(), x.getFeSpace()->getAdmin());
696

Thomas Witkowski's avatar
Thomas Witkowski committed
697
698
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
699
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
700
      (*vecIterator) *= scal; 
701

702
703
704
705
706
707
708
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
709
710
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
711
712
713
714
    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()))
715
      ("no admin or different admins: %8X, %8X\n",
716
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
717
718
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
719
720
    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
721
722
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
723
724
      *xIterator += *yIterator; 

725
726
727
728
729
730
    return x;
  }

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

733
734
735
736
    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()))
737
      ("no admin or different admins: %8X, %8X\n",
738
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
739
740
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
741
742
    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
743
744
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
745
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
746

747
748
749
750
751
752
    return x;
  }

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

    return x;
  }

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

778
779
780
    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()))
781
      ("no admin or different admins: %8X, %8X\n",
782
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
783
784
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
785
    T dot = 0;
786
787
788

    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
789
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
790
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
791
      dot += (*xIterator) * (*yIterator);
792

Thomas Witkowski's avatar
Thomas Witkowski committed
793
    return dot;
794
795
796
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
797
798
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
799
  {
800
801
802
    // Unfortunately needed
    using namespace mtl;

803
804
    FUNCNAME("DOFVector<T>::mv()");

805
806
807
808
    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());
809

810
811
    const DOFAdmin *rowAdmin = a.getRowFeSpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFeSpace()->getAdmin();
812
813

    TEST_EXIT((rowAdmin && colAdmin) && 
814
815
816
817
	      (((transpose == NoTranspose) && (colAdmin == x.getFeSpace()->getAdmin()) && 
		(rowAdmin == result.getFeSpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFeSpace()->getAdmin()) && 
		(colAdmin == result.getFeSpace()->getAdmin()))))
818
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
819
820
       a.getRowFeSpace()->getAdmin(), a.getColFeSpace()->getAdmin(), 
       x.getFeSpace()->getAdmin(), result.getFeSpace()->getAdmin());
821
822


823
824
825
826
827
    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 
828
829
830
831
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
832
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
833
834
835
  {
    FUNCNAME("DOFVector<T>::axpy()");

836
837
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
838

839
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
840

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

  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>&