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
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
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
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
  {
    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
    }
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
710
      if (this->boundaryManager) 
	delete this->boundaryManager; 

711
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
712
    } else {
713
      this->boundaryManager = NULL;
714
715
    }

716
717
718
719
720
721
    return *this;
  }

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

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

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
728
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
729
      (*vecIterator) *= scal; 
730

731
732
733
734
735
736
737
    return x;
  }


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

756
757
758
759
760
761
    return x;
  }

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

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

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

    return x;
  }

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

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

    dot = 0;
821
822
823

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

    return(dot);
  }

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

844
845
    FUNCNAME("DOFVector<T>::mv()");

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

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


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

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

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

881
882
883
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT((admin) && (admin == y.getFESpace()->getAdmin()))
884
885
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
886
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
887
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
888
889
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
890
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
       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];
	}
      }
915
916
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
  }

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

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

949
950
951
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

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

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

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

1034
1035
1036
  }

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

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

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

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

    return result;
1059
1060
1061
  }

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

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

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

1079
1080
1081
1082
    Element *el = elInfo->getElement();

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

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

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

1089
    static T *localvec = NULL;
1090

1091
1092
1093
1094
1095
    T *result;

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

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

1115
1116
1117
    return const_cast<const T*>(result);
  }

1118
1119
1120
1121
1122
1123
1124



  // 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
1125
  {
1126
1127
    return v.getSize();
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
1128

1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
  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
1148
1149
  }

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

    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
1164
1165
1166
    Mesh* mesh = feSpace->getMesh();   
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, 
		   DoubleWell_fct);
1167
1168
1169
1170
1171
1172
1173

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::DoubleWell_fct(ElInfo *elinfo)
  {   
Thomas Witkowski's avatar
Thomas Witkowski committed
1174
1175
    const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int nPoints = quad_fast->getNumPoints();
1176
    double normT = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1177
1178
    for (int iq = 0; iq < nPoints; iq++) {
      normT += quad_fast->getWeight(iq) * sqr(uh_vec[iq]) * sqr(1.0 - uh_vec[iq]);