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

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

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

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

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

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

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

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


} // namespace mtl



68
69
70
namespace AMDiS {

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

80
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
81
    localVectors.resize(omp_get_overall_max_threads());
82
    grdPhis.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
83
    grdTmp.resize(omp_get_overall_max_threads());
84
    D2Phis.resize(omp_get_overall_max_threads());
85

86
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
87
88
      localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
      localVectors[i] = new T[this->nBasFcts];
89
90
91
      grdPhis[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      grdTmp[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      D2Phis[i] = new DimMat<double>(dim, NO_INIT);
92
    }    
93
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
94
  
95
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
96
97
98
  DOFVectorBase<T>::~DOFVectorBase()
  {
    for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
99
100
      delete [] localIndices[i];
      delete [] localVectors[i];
101
102
103
      delete grdPhis[i];
      delete grdTmp[i];
      delete D2Phis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
104
105
106
107
    }
  }
  
  template<typename T>
108
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
109
    : DOFVectorBase<T>(f, n),
110
      coarsenOperation(NO_OPERATION)
111
  {
112
    vec.resize(0);
113
114
115
116
    init(f, n);
  } 

  template<typename T>
117
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
118
  {
119
120
121
122
123
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
    this->boundaryManager = new BoundaryManager(f);
124
125
126
127
128
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
129
130
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
131

Thomas Witkowski's avatar
Thomas Witkowski committed
132
    if (this->boundaryManager)
133
      delete this->boundaryManager;
134
135

    vec.clear();
136
137
138
139
140
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
143
144
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
145
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
146
147
148
149
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

150
151
152
    std::vector<DegreeOfFreedom> &indices = elInfo->getLocalIndices(feSpace);
    TEST_EXIT_DBG(static_cast<int>(indices.size()) == nBasFcts)
      ("Local index vector is too small!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
153

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

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

161
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
162
	  (*this)[irow] += factor * elVec[i];
163
164
	else  
	  (*this)[irow] = factor * elVec[i];	
Thomas Witkowski's avatar
Thomas Witkowski committed
165
166
167
168
      }
    }
  }

169
170
171
172
173
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

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

176
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
177
178
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
179
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
180
181
182
183
184
185
186
187
188
189
      nrm += (*vecIterator) * (*vecIterator);

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
204
    FUNCNAME("DOFVector<T>::asum()");
205

206
    checkFeSpace(this->feSpace, vec);
207

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

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
220
    FUNCNAME("DOFVector<T>::sum()");
221

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

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

    return(nrm);
  }

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

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

240
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
241
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
242
243
244
245
246
247
248
      *vecIterator = alpha ; 
  }


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

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

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


  template<typename T>
  T DOFVector<T>::min() const
  {
270
271
    FUNCNAME("DOFVector<T>::min()");
    
272
    checkFeSpace(this->feSpace, vec);
273
274

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

279
280
281
282
283
284
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
285
286
    FUNCNAME("DOFVector<T>::max()");
    
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::max(m, *vecIterator);
293

294
295
296
    return m;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
297
298
299
300
301
302
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

303
304
305
  template<typename T>
  void DOFVector<T>::print() const
  {
306
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
307

308
    const DOFAdmin *admin = NULL;
309
    const char *format;
310

311
312
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
313

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

354
355
356
357
358
359
360
361
362
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
363
364
365
366
367
368

  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda, 
			     DegreeOfFreedom* dof_indices)
  {
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFESpace()->getBasisFcts());
369
    int nBasisFcts = phi->getNumber();
370
371
    T val = 0.0;

372
    for (int i = 0; i < nBasisFcts; i++)
373
374
375
376
377
378
379
380
      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
381
    FUNCNAME("DOFVector::interpol()");
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
384
385

    interFct = fct;
386

387
388
389
390
391
    if (!this->getFESpace()) {
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
392

393
394
395
396
397
    if (!(this->getFESpace()->getAdmin())) {
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
	  this->getFESpace()->getName().c_str());
      return;
    }
398
  
399
400
401
402
403
    if (!(this->getFESpace()->getBasisFcts())) {
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
404

405
406
407
408
409
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
410
411
412

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
415
416
417
418
419
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this->getFESpace()->getMesh(), -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
420
421
422
  }

  template<typename T>
423
  void DOFVector<T>::interpolFct(ElInfo* elinfo)
424
425
426
  {
    const BasisFunction *basFct = traverseVector->getFESpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFESpace()->getAdmin();
427
428
429
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
430
431
432

    DegreeOfFreedom *myLocalIndices = this->localIndices[omp_get_thread_num()];
    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), admin, myLocalIndices);
433
434
    int nBasFcts = basFct->getNumber();
    for (int i = 0; i < nBasFcts; i++)
435
      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
436
437
438
439
440
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
445
    if (!q) {
446
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
449

450
    FastQuadrature *quadFast = 
451
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
452

453
454
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
Thomas Witkowski's avatar
Thomas Witkowski committed
455
    std::vector<T> uh_vec(nPoints);
456
    TraverseStack stack;
457
    stack.addFeSpace(this->feSpace);
458
459
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
460
461
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
462
463
464
465
466
467
468
    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;
469

470
471
      elInfo = stack.traverseNext(elInfo);
    }
472

473
    return result;  
474
475
476
477
478
  }

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

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

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

491
492
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
493
    std::vector<T> uh_vec(nPoints);
494
    TraverseStack stack;
495
    stack.addFeSpace(this->feSpace);
496
497
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
498
499
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
500
501
502
503
504
505
506
    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;
507

508
509
      elInfo = stack.traverseNext(elInfo);
    }
510

511
    return result;  
512
513
514
515
516
  }

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

519
520
    Mesh* mesh = this->feSpace->getMesh();

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

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

529
530
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
531
    std::vector<T> uh_vec(nPoints);
532
    TraverseStack stack;
533
    stack.addFeSpace(this->feSpace);
534
535
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
536
537
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
538
539
540
541
542
543
544
    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;
545

546
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
547
    }
548

549
    return result;  
550
551
  }

552
  template<typename T>  
553
554
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
555
    FUNCNAME("DOFVector::H1NormSquare()");
556

557
558
    Mesh* mesh = this->feSpace->getMesh();

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

564
565
566
567
568
569
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
570
    std::vector<WorldVector<T> > grduh_vec(nPoints);
571
572
573
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
574
575
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
576
577
578
579
580
581
582
583
584
585
586
    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;
      }
587

588
      result += det * normT;
589

590
591
      elInfo = stack.traverseNext(elInfo);
    }
592

593
    return result;  
594
595
596
597
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
598
					std::vector<DegreeOfFreedom> &newDOF)
599
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
600
601
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
602
	vec[newDOF[i]] = vec[i];
603
604
605
606
607
608
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
609
610
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
611
	 op != this->operators.end(); ++op)
612
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
613

614
615
616
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
617
618
619
620
621
  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
622
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
623
624
625
      (*it)->finishAssembling();
  }

626
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
627
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
628
  {
629
    this->feSpace = rhs.feSpace;
630
    this->nBasFcts = rhs.nBasFcts;
631
    vec = rhs.vec;
632
    this->elementVector.change_dim(this->nBasFcts);
633
634
635
636
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
637
    
638
    if (rhs.boundaryManager) {
639
640
641
      if (this->boundaryManager) 
	delete this->boundaryManager; 

642
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
643
    } else {
644
      this->boundaryManager = NULL;
645
646
    }

647
648
649
650
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

651
652
653
654
655
656
    return *this;
  }

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

659
660
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
661

Thomas Witkowski's avatar
Thomas Witkowski committed
662
663
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
664
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
665
      (*vecIterator) *= scal; 
666

667
668
669
670
671
672
673
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
674
675
676
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
677
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
678
679
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
680
681
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
682
683
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
684
685
    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
686
687
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
688
689
      *xIterator += *yIterator; 

690
691
692
693
694
695
    return x;
  }

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

698
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
699
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
700
701
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
702
703
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
704
705
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
706
707
    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
708
709
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
710
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
711

712
713
714
715
716
717
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
718
719
720
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
721
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
722
723
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
724
725
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
726
727
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
728
729
    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
730
731
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
732
      *xIterator *= *yIterator; 
733
734
735
736
737
738
739
740

    return x;
  }

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

743
    TEST_EXIT(x.getFESpace() && y.getFESpace())
744
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
745
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) && (admin == y.getFESpace()->getAdmin()))
746
747
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
748
749
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
750
    T dot = 0;
751
752
753

    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
754
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
755
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
756
      dot += (*xIterator) * (*yIterator);
757

Thomas Witkowski's avatar
Thomas Witkowski committed
758
    return dot;
759
760
761
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
762
763
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
764
  {
765
766
767
    // Unfortunately needed
    using namespace mtl;

768
769
    FUNCNAME("DOFVector<T>::mv()");

Thomas Witkowski's avatar
Thomas Witkowski committed
770
771
    TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && 
	      x.getFESpace() && result.getFESpace())
772
773
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
774
775
776
777
778
779
780
781
782

    const DOFAdmin *rowAdmin = a.getRowFESpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFESpace()->getAdmin();

    TEST_EXIT((rowAdmin && colAdmin) && 
	      (((transpose == NoTranspose) && (colAdmin == x.getFESpace()->getAdmin()) && 
		(rowAdmin == result.getFESpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFESpace()->getAdmin()) && 
		(colAdmin == result.getFESpace()->getAdmin()))))
783
784
785
786
787
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());


788
789
790
791
792
    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 
793
794
795
796
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
797
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
798
799
800
  {
    FUNCNAME("DOFVector<T>::axpy()");

801
    TEST_EXIT(x.getFESpace() && y.getFESpace())
802
803
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

804
805
806
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT((admin) && (admin == y.getFESpace()->getAdmin()))
807
808
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
809
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
810
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
811
812
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
813
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
       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
833
834
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
835
	  y[i] += alpha * x[i];
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
  }

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

867
    TEST_EXIT(x.getFESpace() && y.getFESpace())
868
869
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

870
871
872
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

  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
911
912
913
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = scal * (*vIterator);
914
915
916
917
918
919
920
921
922
923
924

    return result;