DOFVector.hh 33.9 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] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
      localVectors[i] = GET_MEMORY(T, this->nBasFcts);
89
      grdPhis[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
90
      grdTmp[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
91
      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
99
  DOFVectorBase<T>::~DOFVectorBase()
  {
    for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
      FREE_MEMORY(localIndices[i], DegreeOfFreedom, this->nBasFcts);
Thomas Witkowski's avatar
Thomas Witkowski committed
100
      FREE_MEMORY(localVectors[i], T, this->nBasFcts);
101
      DELETE grdPhis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
102
      DELETE grdTmp[i];
103
      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
  {
    init(f, n);
113

114
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
115
116
    applicationOrdering = NULL;
#endif
117
118
119
  } 

  template<typename T>
120
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
121
122
123
  {
    name = n;
    feSpace = f;
124
    if(feSpace && feSpace->getAdmin()) {
125
      (feSpace->getAdmin())->addDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
126
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
127
    this->boundaryManager = NEW BoundaryManager(f);
128
129
130
131
132
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
133
    if (feSpace && feSpace->getAdmin())
134
135
      (feSpace->getAdmin())->removeDOFIndexed(this);

Thomas Witkowski's avatar
Thomas Witkowski committed
136
    if (this->boundaryManager)
137
      DELETE this->boundaryManager;
138
139

    vec.clear();
140
141
142
143
144
145
146
147
148
149
150
  }

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

  template<typename T>
  FastQuadrature *DOFVector<T>::quad_fast = NULL;

  template<typename T>
  double DOFVector<T>::norm = 0.0;

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

160
161
162
163
    Vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
						feSpace->getAdmin(),
						&indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
164

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

Thomas Witkowski's avatar
Thomas Witkowski committed
169
      if (!(condition && condition->isDirichlet())) {
170
	DegreeOfFreedom irow = indices[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
173
174
175
176
	(*this)[irow] = (add ? (*this)[irow] : 0.0);
	(*this)[irow] += factor * elVec[i];
      }
    }
  }

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

182
183
184
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
185
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
186
187
188
       feSpace->getAdmin()->getUsedSize());
    
    double nrm = 0.0;
189
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
190
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
191
192
193
194
195
196
197
198
199
200
      nrm += (*vecIterator) * (*vecIterator);

    return(sqrt(nrm));
  }

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

201
202
203
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
204
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
205
206
207
       feSpace->getAdmin()->getUsedSize());
    
    double nrm = 0.0;
208
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
209
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
210
211
212
213
214
215
216
217
      nrm += (*vecIterator) * (*vecIterator);

    return nrm;
  }

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

220
221
222
223
224
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       feSpace->getAdmin()->getUsedSize());
225

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

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
237
    FUNCNAME("DOFVector<T>::sum()");
238

239
240
241
242
243
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       feSpace->getAdmin()->getUsedSize());
244

245
    double nrm = 0.0;    
246
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
247
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
248
249
250
251
252
253
254
255
      nrm += *vecIterator;

    return(nrm);
  }

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

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


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
274
275
276
277
278
    FUNCNAME("DOFVector<T>::copy()");
    
    TEST_EXIT_DBG(feSpace && x.feSpace)
      ("feSpace is NULL: %8X, %8X\n", feSpace, x.feSpace);
    TEST_EXIT_DBG(feSpace->getAdmin() && (feSpace->getAdmin() == x.feSpace->getAdmin()))
279
280
      ("no admin or different admins: %8X, %8X\n",
       feSpace->getAdmin(), x.feSpace->getAdmin());
281
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
282
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
283
284
       feSpace->getAdmin()->getUsedSize());
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >= feSpace->getAdmin()->getUsedSize())
285
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
286
287
       feSpace->getAdmin()->getUsedSize());
    
288
289
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
290
291
292
293
294
    for (vecIterator.reset(), xIterator.reset();
	 !vecIterator.end();
	 ++vecIterator, ++xIterator) {
      
      *vecIterator = *xIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
295
    }
296
297
298
299
300
301
  }


  template<typename T>
  T DOFVector<T>::min() const
  {
302
303
304
305
306
    FUNCNAME("DOFVector<T>::min()");
    
    TEST_EXIT_DBG(feSpace && feSpace->getAdmin())
      ("pointer is NULL: %8X, %8X\n", this, feSpace->getAdmin());
    TEST_EXIT_DBG((static_cast<int>(vec.size())) >= feSpace->getAdmin()->getUsedSize())
307
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
308
309
310
       feSpace->getAdmin()->getUsedSize());

    T m;    
311
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
312
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
313
      m = std::min(m, *vecIterator);
314
315
    }

316
317
318
319
320
321
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
322
323
324
325
326
    FUNCNAME("DOFVector<T>::max()");
    
    TEST_EXIT_DBG(feSpace && feSpace->getAdmin())
      ("pointer is NULL: %8X, %8X\n", this, feSpace->getAdmin());
    TEST_EXIT_DBG((static_cast<int>(vec.size())) >= feSpace->getAdmin()->getUsedSize())
327
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
328
       feSpace->getAdmin()->getUsedSize());
329

330
    T m;    
331
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
332
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
333
      m = std::max(m, *vecIterator);
334
335
    }

336
337
338
    return m;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
341
342
343
344
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

345
346
347
  template<typename T>
  void DOFVector<T>::print() const
  {
348
    FUNCNAME("DOFVector<T>::print()");
349
    const DOFAdmin *admin = NULL;
350
    const char *format;
351

352
353
    if (feSpace) 
      admin = feSpace->getAdmin();
354
355

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

395
396
397
398
399
400
401
402
403
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
404
405
406
407
408
409
410
411
412

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

413
    for (int i = 0; i < numberOfBasFcts; i++)
414
415
416
417
418
419
420
421
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

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

424
425
426
    TEST_EXIT_DBG(fct)("no function to interpolate\n");

    interFct = fct;
427

428
429
430
431
432
    if (!this->getFESpace()) {
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
433

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

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

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
454
455
456
457
458
459
460
    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);
    }
461
462
463
464
465
466
467
  }

  template<typename T>
  int DOFVector<T>::interpolFct(ElInfo* elinfo)
  {
    const BasisFunction *basFct = traverseVector->getFESpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFESpace()->getAdmin();
Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
470
471
472
473
474
    const DegreeOfFreedom *dof =  basFct->getLocalIndices(const_cast<Element *>(elinfo->getElement()),
							  admin, NULL);
    const T *inter_val = const_cast<BasisFunction*>(basFct)->interpol(elinfo, 
								      0, 
								      NULL,
								      traverseVector->interFct, 
								      NULL);
475
476

    int number = basFct->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
477
    for (int i = 0; i < number; i++)
478
479
480
481
482
483
484
485
      (*traverseVector)[dof[i]] = inter_val[i];

    return 0;
  }

  template<typename T>
  double DOFVector<T>::Int(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
486
    FUNCNAME("DOFVector::Int()");
487
488
489
  
    Mesh* mesh = feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
490
491
492
493
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *q, INIT_PHI);
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   Int_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::Int_fct(ElInfo *elinfo)
  {
    double det, normT;
    const T *uh_vec;
    int iq;

    det = elinfo->getDet();
516
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
517
518

    int numPoints = quad_fast->getNumPoints();
519
520
521
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*(uh_vec[iq]);
    }
522
523
524
525
526
527
528
529
    norm += det*normT;

    return 0;
  }

  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
530
    FUNCNAME("DOFVector::L1Norm()");
531
  
Thomas Witkowski's avatar
Thomas Witkowski committed
532
533
534
535
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
536
537
538
539

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),*q,INIT_PHI);
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
540
    Mesh* mesh = feSpace->getMesh();
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   L1Norm_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L1Norm_fct(ElInfo *elinfo)
  {
    double det, normT;
    const T *uh_loc, *uh_vec;
    int iq;

    det = elinfo->getDet();
559
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
560
561

    int numPoints = quad_fast->getNumPoints();
562
563
564
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*abs(uh_vec[iq]);
    }
565
566
567
568
569
570
571
572
573
    norm += det*normT;

    return 0;
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
576
577
578
579
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
580
581
582
583
584
585
586

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, 
						      INIT_PHI);

    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
587
588
    Mesh* mesh = feSpace->getMesh();
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, L2Norm_fct);
589
590
591
592
593
594
595
596
597
598
599
600

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L2Norm_fct(ElInfo *elinfo)
  {
    double det, normT;
    const T *uh_vec;
    int iq;

    det = elinfo->getDet();
601
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
602
603

    int numPoints = quad_fast->getNumPoints();
604
605
606
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq]);
    }
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
    norm += det*normT;

    return 0;
  }

  template<typename T>
  int DOFVector<T>::H1Norm_fct(ElInfo *elinfo)
  {
    double norm2, normT;
    const WorldVector<T> *grduh_vec;
    int iq, j;

    double det = elinfo->getDet();

    grduh_vec = traverseVector->getGrdAtQPs(elinfo, NULL, quad_fast, NULL);

    int dimOfWorld = Global::getGeo(WORLD);

    int numPoints = quad_fast->getNumPoints();
    for (normT = iq = 0; iq < numPoints; iq++)
      {
	for (norm2 = j = 0; j < dimOfWorld; j++)
	  norm2 += sqr(grduh_vec[iq][j]);

	normT += quad_fast->getWeight(iq)*norm2;
      }
    norm += det*normT;

    return 0;
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
646
647
648
649
650
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree() - 2;
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, INIT_GRD_PHI);
651
652

    norm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
    traverseVector = const_cast<DOFVector<T>*>(this); 
    Mesh *mesh = feSpace->getMesh();
655
656
657
658
659
660
661
662
663
664

    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
		   Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA,
		   H1Norm_fct);

    return norm;
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
665
					std::vector<DegreeOfFreedom> &newDOF)
666
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
667
668
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
669
	vec[newDOF[i]] = vec[i];
670
671
672
673
674
675
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
678
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
	 op != this->operators.end(); ++op) {
679
680
      fillFlag |= (*op)->getFillFlag();
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
681

682
683
684
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
685
686
687
688
689
690
691
692
693
694
  template<typename T>
  void DOFVectorBase<T>::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = this->operators.begin();
	 it != this->operators.end(); ++it) {     
      (*it)->finishAssembling();
    }
  }

695
  template<typename T>
696
697
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs )
  {
698
    feSpace = rhs.feSpace;
699
    DOFVectorBase<T>::feSpace = rhs.feSpace;
700
    this->nBasFcts = rhs.nBasFcts;
701
    vec = rhs.vec;
702
    this->elementVector.change_dim(this->nBasFcts);
703
704
705
706
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
707
    if (rhs.boundaryManager) {
708
709
      if (this->boundaryManager) 
	delete this->boundaryManager; 
710
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
711

712
      //    boundaryManager->setDOFVector(this);
713
    } else {
714
      this->boundaryManager = NULL;
715
716
    }

717
718
719
720
721
722
    return *this;
  }

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

725
726
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
727
728

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
729
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
730
      (*vecIterator) *= scal; 
Thomas Witkowski's avatar
Thomas Witkowski committed
731
    }
732

733
734
735
736
737
738
739
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
740
741
742
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
743
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
744
745
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
746
747
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
748
749
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
750
751
    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);
752
753
754
755
756
757
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {     
      *xIterator += *yIterator; 
    }

758
759
760
761
762
763
    return x;
  }

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

766
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
767
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
768
769
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
770
771
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
772
773
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
774
775
    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);
776
777
778
779
780
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {
      *xIterator -= *yIterator; 
    }
781
782
783
784
785
786
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
787
788
789
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
790
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
791
792
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
793
794
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
795
796
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
797
798
    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);
799
800
801
802
803
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {
      *xIterator *= *yIterator; 
    }
804
805
806
807
808
809
810
811

    return x;
  }

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

815
    TEST_EXIT(x.getFESpace() && y.getFESpace())
816
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
817
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) && (admin == y.getFESpace()->getAdmin()))
818
819
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
820
821
822
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

    dot = 0;
823
824
825

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
826
827
828
829
830
831
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) 
      {
	dot += (*xIterator) * (*yIterator);
      };
832
833
834
835
836
837
838
839
840
841
842

    return(dot);
  }

  template<typename T>
  void mv(MatrixTranspose transpose,
	  const DOFMatrix &a, 
	  const DOFVector<T>&x,
	  DOFVector<T> &result,
	  bool add)
  {
843
844
845
    // Unfortunately needed
    using namespace mtl;

846
847
    FUNCNAME("DOFVector<T>::mv()");

848
    TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && x.getFESpace() && result.getFESpace())
849
850
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
851
852
853
854
855
856
857
858
859

    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()))))
860
861
862
863
864
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());


865
866
867
868
869
    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 
870
871
872
873
874
875
876
877
878
879
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
  void axpy(double alpha, 
	    const DOFVector<T>& x, 
	    DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::axpy()");

880
    TEST_EXIT(x.getFESpace() && y.getFESpace())
881
882
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

883
884
885
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT((admin) && (admin == y.getFESpace()->getAdmin()))
886
887
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
888
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
889
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
890
891
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
892
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
       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
      for (i = 0; i < maxI; i++) {
	if (!admin->isDOFFree(i)) {
	  y[i] += alpha * x[i];
	}
      }
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
  }

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

948
    TEST_EXIT(x.getFESpace() && y.getFESpace())
949
950
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

951
952
953
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT(admin && (admin == y.getFESpace()->getAdmin()))
954
955
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
956
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
957
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
958
959
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
960
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
       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
      for (i = 0; i < maxI; i++) {
	if (!admin->isDOFFree(i)) {
	  y[i] = alpha * y[i] + x[i];
	}
      }
985
986
987
988
989
990
991
992
993
  }

  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);
994
995
996
997
998
999
    for(vIterator.reset(), rIterator.reset();
	!vIterator.end();
	++vIterator, ++rIterator) 
      {  
	*rIterator = scal * (*vIterator);
      };
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010

    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);
1011
1012
1013
1014
1015
1016
    for(vIterator.reset(), rIterator.reset();
	!vIterator.end();
	++vIterator, ++rIterator) 
      {  
	*rIterator = (*vIterator) + scal;
      };
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
    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);
1028
1029
1030
1031
1032
1033
    for(v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	!v1Iterator.end();
	++v1Iterator, ++v2Iterator, ++rIterator) 
      {  
	*rIterator = (*v1Iterator) + (*v2Iterator);
      };
1034
    return result;
1035

1036
1037
1038
  }

  template<typename T>
1039
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
1040
  {
1041
1042
    static T* localVec = NULL;
    const DOFAdmin* admin = feSpace->getAdmin();
1043

1044
1045
    T *result;
  
1046
    if (d) {
1047
1048
1049
1050
1051
1052
      result = d;
    } else {
      if(localVec) delete [] localVec;
      localVec = new T[nBasFcts]; 
      result = localVec;
    }
1053

1054
1055
    const DegreeOfFreedom *localIndices = 
      feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);
1056

1057
    for (int i = 0; i < nBasFcts; i++)
1058
      result[i] = (*this)[localIndices[i]];
1059
1060

    return result;
1061
1062
1063
  }

  template<typename T>
1064
1065
  const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo         *elInfo, 
					 const Quadrature     *quad,
1066
					 const FastQuadrature *quadFast,
1067
					 T                    *vecAtQPs) const
1068
1069
1070
  {
    FUNCNAME("DOFVector<T>::getVecAtQPs()");
  
1071
    TEST_EXIT(quad || quadFast)("neither quad nor quadFast defined\n");
1072

1073
1074
1075
    if(quad && quadFast) {
      TEST_EXIT(quad == quadFast->getQuadrature())
	("quad != quadFast->quadrature\n");
1076
1077
    }

1078
    TEST_EXIT(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
1079
1080
      ("invalid basis functions");

1081
1082
1083
1084
    Element *el = elInfo->getElement();

    const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;

1085
    const BasisFunction *basFcts = feSpace->getBasisFcts();
1086
1087
1088
1089
1090

    int numPoints = quadrature->getNumPoints();
    int nBasFcts  = basFcts->getNumber();
    int i, j;

1091
    static T *localvec = NULL;
1092

1093
1094
1095
1096
1097
    T *result;

    if (vecAtQPs) {
      result = vecAtQPs;
    } else {
1098
1099
1100
      if(localvec) delete [] localvec;
      localvec = new T[numPoints];
      for(i = 0; i < numPoints; i++) {
1101
1102
1103
1104
	localvec[i] = 0.0;
      }
      result = localvec;
    }
1105
1106
1107
1108
1109
1110
  
    const T *localVec = getLocalVector(el, NULL);

    for(i = 0; i < numPoints; i++) {
      for(result[i] = j = 0; j < nBasFcts; j++) {
	result[i] +=  localVec[j] * 
1111
1112
1113
1114
1115
	  (quadFast ? 
	   (quadFast->getPhi(i, j)) : 
	   ((*(basFcts->getPhi(j)))(quad->getLambda(i))));
      }
    }
1116

1117
1118
1119
    return const_cast<const T*>(result);
  }

1120
1121
1122
1123
1124
1125
1126



  // Some free functions used in MTL4

  template <typename T>
  inline std::size_t size(const AMDiS::DOFVector<T>& v)
Thomas Witkowski's avatar
Thomas Witkowski committed
1127
  {
1128
1129
    return v.getSize();
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
1130

1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
  template <typename T>
  inline std::size_t num_rows(const AMDiS::DOFVector<T>& v)
  {
    return v.getSize();
  }

  template <typename T>
  inline std::size_t num_cols(const AMDiS::DOFVector<T>& v)
  {
    return 1;
  }
 
  template <typename T>
  inline void set_to_zero(AMDiS::DOFVector<T>& v)
  {
    using math::zero;
    T  ref, my_zero(zero(ref));

    std::fill(v.begin(), v.end(), my_zero);
Thomas Witkowski's avatar
Thomas Witkowski committed
1150
1151
  }

1152
1153
1154
1155
1156
1157
  template<typename T>
  double DOFVector<T>::DoubleWell(Quadrature* q) const
  {
    FUNCNAME("DOFVector::DoubleWell()");
  
    if (!q) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1158
1159
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
1160
1161
1162
1163
1164
1165
    }

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, INIT_PHI);
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
1166
1167
1168
    Mesh* mesh = feSpace->getMesh();   
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, 
		   DoubleWell_fct);