DOFVector.hh 31.7 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
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
123
    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
153
    Vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
						feSpace->getAdmin(),
						&indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
154

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

Thomas Witkowski's avatar
Thomas Witkowski committed
159
      if (!(condition && condition->isDirichlet())) {
160
	DegreeOfFreedom irow = indices[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
161
162
163
164
165
166

	if (add) {
	  (*this)[irow] += factor * elVec[i];
	} else {  
	  (*this)[irow] = factor * elVec[i];
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
169
170
      }
    }
  }

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

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

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

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
206
    FUNCNAME("DOFVector<T>::asum()");
207

208
    checkFeSpace(this->feSpace, vec);
209

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

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
222
    FUNCNAME("DOFVector<T>::sum()");
223

224
    checkFeSpace(this->feSpace, vec);
225

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

    return(nrm);
  }

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

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

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


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

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

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


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

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

281
282
283
284
285
286
    return m;
  }

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

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

296
297
298
    return m;
  }

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

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

310
    const DOFAdmin *admin = NULL;
311
    const char *format;
312

313
314
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
315

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

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

    return result;
  }
365
366
367
368
369
370

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

374
    for (int i = 0; i < nBasisFcts; i++)
375
376
377
378
379
380
381
382
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

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

385
386
387
    TEST_EXIT_DBG(fct)("no function to interpolate\n");

    interFct = fct;
388

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

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

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

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
419
420
421
    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);
    }
422
423
424
  }

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

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

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

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

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

455
456
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
Thomas Witkowski's avatar
Thomas Witkowski committed
457
    std::vector<T> uh_vec(nPoints);
458
459
460
461
462
463
464
465
466
467
468
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
    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
495
496
497
498
499
500
501
502
503
504
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
    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;
505

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

509
    return result;  
510
511
512
513
514
  }

  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
515
    FUNCNAME("DOFVector::L2NormSquare()");
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
    std::vector<T> uh_vec(nPoints);
530
531
532
533
534
535
536
537
538
539
540
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
    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;
541

542
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
543
    }
544

545
    return result;  
546
547
  }

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

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

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

560
561
562
563
564
565
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

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

584
      result += det * normT;
585

586
587
      elInfo = stack.traverseNext(elInfo);
    }
588

589
    return result;  
590
591
592
593
  }

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

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

610
611
612
    return fillFlag;
  }

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

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

637
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
638
    } else {
639
      this->boundaryManager = NULL;
640
641
    }

642
643
644
645
646
647
    return *this;
  }

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

650
651
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
652

Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
655
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
656
      (*vecIterator) *= scal; 
657

658
659
660
661
662
663
664
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
665
666
667
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
668
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
669
670
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
671
672
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
673
674
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
675
676
    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
677
678
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
679
680
      *xIterator += *yIterator; 

681
682
683
684
685
686
    return x;
  }

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

689
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
690
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
691
692
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
693
694
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
695
696
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
697
698
    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
699
700
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
701
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
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
711
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
712
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
713
714
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
715
716
      ("no admin or different admins: %8X, %8X\n",
       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
      *xIterator *= *yIterator; 
724
725
726
727
728
729
730
731

    return x;
  }

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

734
    TEST_EXIT(x.getFESpace() && y.getFESpace())
735
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
736
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) && (admin == y.getFESpace()->getAdmin()))
737
738
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
739
740
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
741
    T dot = 0;
742
743
744

    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
745
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
746
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
747
      dot += (*xIterator) * (*yIterator);
748

Thomas Witkowski's avatar
Thomas Witkowski committed
749
    return dot;
750
751
752
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
753
754
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
755
  {
756
757
758
    // Unfortunately needed
    using namespace mtl;

759
760
    FUNCNAME("DOFVector<T>::mv()");

Thomas Witkowski's avatar
Thomas Witkowski committed
761
762
    TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && 
	      x.getFESpace() && result.getFESpace())
763
764
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
765
766
767
768
769
770
771
772
773

    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()))))
774
775
776
777
778
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());


779
780
781
782
783
    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 
784
785
786
787
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
788
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
789
790
791
  {
    FUNCNAME("DOFVector<T>::axpy()");

792
    TEST_EXIT(x.getFESpace() && y.getFESpace())
793
794
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

795
796
797
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

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

858
    TEST_EXIT(x.getFESpace() && y.getFESpace())
859
860
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

861
862
863
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

  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
902
903
904
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = scal * (*vIterator);
905
906
907
908
909
910
911
912
913
914
915

    return result;
  }

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
916
917
918
919
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = (*vIterator) + scal;

920
921
922
923
924
925
926
927
928
929
930
    return result;
  }

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v1,
				 const DOFVector<T>& v2,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator v1Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v1)), USED_DOFS);
    typename DOFVector<T>::Iterator v2Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v2)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
931
932
933
    for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	 !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator) 
      *rIterator = (*v1Iterator) + (*v2Iterator);
934

Thomas Witkowski's avatar
Thomas Witkowski committed
935
    return result;
936
937
938
  }

  template<typename T>
939
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
940
  {
941
    static T* localVec = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
942
    static int localVecSize = 0;
943
944
945
    const DOFAdmin* admin = feSpace->getAdmin();
    T *result;
  
946
    if (d) {
947
948
      result = d;
    } else {
949
950
951
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
952
953
954
955
956
957
958
959
      if (localVec && nBasFcts > localVecSize)  {
	delete [] localVec;
	localVec = new T[nBasFcts]; 
      }
      if (!localVec)
	localVec = new T[nBasFcts]; 

      localVecSize = nBasFcts;
960
961
      result = localVec;
    }
962

963
964
965
966
#ifdef _OPENMP
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(el, admin, &(localIndices[0]));
#else
967
968
    const DegreeOfFreedom *localIndices = 
      feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);
969
#endif
970

971
    for (int i = 0; i < nBasFcts; i++)
972
      result[i] = (*this)[localIndices[i]];
973
974

    return result;
975
976
977
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
978
979
  const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo, 
					 const Quadrature *quad,
980
					 const FastQuadrature *quadFast,
Thomas Witkowski's avatar
Thomas Witkowski committed
981
					 T *vecAtQPs) const
982
983
984
  {
    FUNCNAME("DOFVector<T>::getVecAtQPs()");
  
985
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
986

Thomas Witkowski's avatar
Thomas Witkowski committed
987
    if (quad && quadFast)
988
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
989
	("quad != quadFast->quadrature\n");
990

991
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
992
993
      ("invalid basis functions");

994
995
    Element *el = elInfo->getElement();
    const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
996
    const BasisFunction *basFcts = feSpace->getBasisFcts();
997
    int nPoints = quadrature->getNumPoints();
998
    int nBasFcts  = basFcts->getNumber();
Peter Gottschling's avatar