DOFVector.hh 40.3 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
10
11
12
13
14
15
16
17

#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 "ElementVector.h"
#include "Assembler.h"
18
#include "OpenMP.h"
19
#include "Operator.h"
20
#include "Parameters.h"
21
22
23
24

namespace AMDiS {

  template<typename T>
25
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
26
27
28
29
    : feSpace(f),
      name(n),
      elementVector(NULL),
      boundaryManager(NULL)
30
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
31
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
32
    dim = feSpace->getMesh()->getDim();
33

34
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
35
    localVectors.resize(omp_get_overall_max_threads());
36
    grdPhis.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
37
    grdTmp.resize(omp_get_overall_max_threads());
38
    D2Phis.resize(omp_get_overall_max_threads());
39

40
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
41
42
      localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
      localVectors[i] = GET_MEMORY(T, this->nBasFcts);
43
      grdPhis[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
44
      grdTmp[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
45
      D2Phis[i] = NEW DimMat<double>(dim, NO_INIT);
46
47
    }
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
48
  
49
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
50
51
52
53
  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
54
      FREE_MEMORY(localVectors[i], T, this->nBasFcts);
55
      DELETE grdPhis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
56
      DELETE grdTmp[i];
57
      DELETE D2Phis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
60
61
    }
  }
  
  template<typename T>
62
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
63
    : DOFVectorBase<T>(f, n),
64
      coarsenOperation(NO_OPERATION)
65
66
67
68
69
  {
    init(f, n);
  } 

  template<typename T>
70
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
71
72
73
  {
    name = n;
    feSpace = f;
74
    if (feSpace && feSpace->getAdmin()) {
75
      (feSpace->getAdmin())->addDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
76
77
    }
      
Thomas Witkowski's avatar
Thomas Witkowski committed
78
    this->boundaryManager = NEW BoundaryManager(f);
79
80
81
82
83
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
84
    if (feSpace && feSpace->getAdmin()) {
85
86
87
88
89
90
      (feSpace->getAdmin())->removeDOFIndexed(this);
    }

    if (this->boundaryManager) {
      DELETE this->boundaryManager;
    }
91
92

    vec.clear();
93
94
95
96
97
98
99
100
101
102
103
  }

  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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

    int n_row = elVec.getSize();

    for (DegreeOfFreedom i = 0; i < n_row; i++) {
      BoundaryCondition *condition = 
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;

      if(!(condition && condition->isDirichlet())) {
	DegreeOfFreedom irow = elVec.dofIndices[i];
	(*this)[irow] = (add ? (*this)[irow] : 0.0);
	(*this)[irow] += factor * elVec[i];
      }
    }
  }

126
127
128
129
130
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

131
132
133
    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())
134
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
135
136
137
       feSpace->getAdmin()->getUsedSize());
    
    double nrm = 0.0;
138
139
140
141
142
143
144
145
146
147
148
149
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += (*vecIterator) * (*vecIterator);

    return(sqrt(nrm));
  }

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

150
151
152
    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())
153
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
154
155
156
       feSpace->getAdmin()->getUsedSize());
    
    double nrm = 0.0;
157
158
159
160
161
162
163
164
165
166
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += (*vecIterator) * (*vecIterator);

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
167
    FUNCNAME("DOFVector<T>::asum()");
168

169
170
171
172
173
    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());
174

175
    double nrm = 0.0;
176
177
178
179
180
181
182
183
184
185
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += abs(*vecIterator);

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
186
    FUNCNAME("DOFVector<T>::sum()");
187

188
189
190
191
192
    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());
193

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

    return(nrm);
  }

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

207
208
209
    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())
210
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
211
212
       feSpace->getAdmin()->getUsedSize());
    
213
214
215
216
217
218
219
220
221
222
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
    for(vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
      *vecIterator = alpha ; 
    };
  }


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
223
224
225
226
227
    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()))
228
229
      ("no admin or different admins: %8X, %8X\n",
       feSpace->getAdmin(), x.feSpace->getAdmin());
230
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
231
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
232
233
       feSpace->getAdmin()->getUsedSize());
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >= feSpace->getAdmin()->getUsedSize())
234
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
235
236
       feSpace->getAdmin()->getUsedSize());
    
237
238
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
239
240
241
242
243
244
    for (vecIterator.reset(), xIterator.reset();
	 !vecIterator.end();
	 ++vecIterator, ++xIterator) {
      
      *vecIterator = *xIterator; 
    };
245
246
247
248
249
250
  }


  template<typename T>
  T DOFVector<T>::min() const
  {
251
252
253
254
255
    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())
256
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
257
258
259
       feSpace->getAdmin()->getUsedSize());

    T m;    
260
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
261
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
262
      m = std::min(m, *vecIterator);
263
264
    }

265
266
267
268
269
270
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
271
272
273
274
275
    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())
276
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
277
       feSpace->getAdmin()->getUsedSize());
278

279
    T m;    
280
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
281
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
282
      m = std::max(m, *vecIterator);
283
284
    }

285
286
287
    return m;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
290
291
292
293
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

294
295
296
297
298
  template<typename T>
  void gemv(MatrixTranspose transpose, T alpha,
	    const DOFMatrix& a, const DOFVector<T>& x,
	    T beta, DOFVector<T>& y)
  {
299
300
301
302
    FUNCNAME("DOFVector<T>::gemv()");

    int j, jcol, ysize;
    T sum, ax;
303
    const DOFMatrix::MatrixRow *row;
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
    
    TEST_EXIT_DBG(a.getRowFESpace() && 
		  a.getColFESpace() && 
		  x.getFESpace() && 
		  y.getFESpace())
      ("feSpace is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), y.getFESpace());
    TEST_EXIT_DBG(a.getRowFESpace()->getAdmin() && 
		  a.getColFESpace()->getAdmin() && 
		  (((transpose == NoTranspose) && 
		    (a.getColFESpace()->getAdmin() == x.getFESpace()->getAdmin()) && 
		    (a.getRowFESpace()->getAdmin() == y.getFESpace()->getAdmin()))||
		   ((transpose == Transpose) && 
		    (a.getRowFESpace()->getAdmin() == x.getFESpace()->getAdmin()) && 
		    (a.getColFESpace()->getAdmin() == y.getFESpace()->getAdmin()))))
319
320
321
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
322
    
323
    if (transpose == NoTranspose) {
324
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
325
	("x.size = %d too small: admin->sizeUsed = %d\n",
326
327
	 x.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>(y.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
328
	("y.size = %d too small: admin->sizeUsed = %d\n",
329
330
	 y.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
331
	("a.size = %d too small: admin->sizeUsed = %d\n",
332
	a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
333
334
    }
    else  if (transpose == Transpose) {
335
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
336
	("x.size = %d too small: admin->sizeUsed = %d\n",
337
338
	 x.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>(y.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
339
	("y.size = %d too small: admin->sizeUsed = %d\n",
340
341
	 y.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
342
	("a.size = %d too small: admin->sizeUsed = %d\n",
343
	a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
    }

    ysize = y.getSize();

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), FREE_DOFS);
    for(vecIterator.reset();
	!vecIterator.end(); ++vecIterator) { 
      *vecIterator = 0;
    };

    if (transpose == NoTranspose) {
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS); 
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
      for(vecIterator.reset(), rowIterator.reset();
	  !rowIterator.end();
	  ++rowIterator, ++vecIterator) { 
	sum = 0;
	row = &(a[rowIterator.getDOFIndex()]);
362
	for(std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
	    colIterator != rowIterator->end();
	    colIterator++) {
	  jcol = colIterator->col;
	  if (jcol >= 0) { // entry used? 
	    sum += (static_cast<T>(colIterator->entry)) * x[jcol];
	  } else {
	    if (jcol == DOFMatrix::NO_MORE_ENTRIES)
	      break;
	  }
	}
	*vecIterator *= beta;
	*vecIterator += alpha * sum;
      };
    } else if (transpose == Transpose) {
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
      for(vecIterator.reset();
	  !vecIterator.end(); 
	  ++vecIterator) {
	*vecIterator  *= beta ; 
      };
    
      typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
      for(xIterator.reset(), rowIterator.reset();
	  !rowIterator.end();
	  ++rowIterator, ++xIterator) { 
	ax = alpha * (*xIterator);
	row = &(a[rowIterator.getDOFIndex()]);
391
	for(std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
	    colIterator != rowIterator->end();
	    colIterator++) {
	  jcol = colIterator->col;
	  if (jcol >= 0) // entry used?
	    y[jcol] += ax * (static_cast<T>(colIterator->entry));
	  else 
	    if (jcol == DOFMatrix::NO_MORE_ENTRIES) break;
	}
      }
    }
    else {
      ERROR_EXIT("transpose=%d\n", transpose);
    }
  }

  template<typename T>
  void DOFVector<T>::print() const
  {
410
    FUNCNAME("DOFVector<T>::print()");
411
    const DOFAdmin *admin = NULL;
412
    const char *format;
413

414
415
    if (feSpace) 
      admin = feSpace->getAdmin();
416
417

    MSG("Vec `%s':\n", name.c_str());
418
    int j = 0;
419
420
421
422
423
424
425
426
427
428
429
430
    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);
      for(vecIterator.reset();
	  !vecIterator.end(); ++vecIterator) {
	if ((j % 3) == 0) {
431
432
	  if (j) 
	    Msg::print("\n");
433
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
434
	} else 
435
436
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
	j++;
437
      }
438
      Msg::print("\n");
439
    } else {
440
441
	MSG("no DOFAdmin, print whole vector.\n");
    
442
443
444
445
446
447
	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 {
448
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
449
	  }
450
451
452
453
454
455
456
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }

457
458
459
460
461
462
463
464
465
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
466
467
468
469
470
471
472
473
474

  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;

475
    for (int i = 0; i < numberOfBasFcts; i++)
476
477
478
479
480
481
482
483
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

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

486
487
488
    TEST_EXIT_DBG(fct)("no function to interpolate\n");

    interFct = fct;
489

490
491
492
493
494
    if (!this->getFESpace()) {
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
495

496
497
498
499
500
    if (!(this->getFESpace()->getAdmin())) {
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
	  this->getFESpace()->getName().c_str());
      return;
    }
501
  
502
503
504
505
506
    if (!(this->getFESpace()->getBasisFcts())) {
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
507

508
509
510
511
512
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548

    traverseVector = this;
    this->getFESpace()->getMesh()->traverse(-1, 
					    Mesh::CALL_LEAF_EL |
					    Mesh::FILL_COORDS, 
					    interpolFct);

    return;
  }

  template<typename T>
  int DOFVector<T>::interpolFct(ElInfo* elinfo)
  {
    int i;
    const BasisFunction *basFct = traverseVector->getFESpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFESpace()->getAdmin();
    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);

    int number = basFct->getNumber();
    for (i = 0; i < number; i++)
      (*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
549
    FUNCNAME("DOFVector::Int()");
550
551
552
  
    Mesh* mesh = feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
553
554
555
556
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578

    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();
579
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
580
581

    int numPoints = quad_fast->getNumPoints();
582
583
584
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*(uh_vec[iq]);
    }
585
586
587
588
589
590
591
592
    norm += det*normT;

    return 0;
  }

  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
593
    FUNCNAME("DOFVector::L1Norm()");
594
  
Thomas Witkowski's avatar
Thomas Witkowski committed
595
596
597
598
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
599
600
601
602

    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
603
    Mesh* mesh = feSpace->getMesh();
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621

    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();
622
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
623
624

    int numPoints = quad_fast->getNumPoints();
625
626
627
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*abs(uh_vec[iq]);
    }
628
629
630
631
632
633
634
635
636
    norm += det*normT;

    return 0;
  }


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

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

    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
650
651
    Mesh* mesh = feSpace->getMesh();
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, L2Norm_fct);
652
653
654
655
656
657
658
659
660
661
662
663

    return norm;  
  }

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

    det = elinfo->getDet();
664
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
665
666

    int numPoints = quad_fast->getNumPoints();
667
668
669
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq]);
    }
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
    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
705
    FUNCNAME("DOFVector::H1NormSquare()");
706

Thomas Witkowski's avatar
Thomas Witkowski committed
707
708
709
710
711
712
713
    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);
714
715

    norm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
    traverseVector = const_cast<DOFVector<T>*>(this); 
    Mesh *mesh = feSpace->getMesh();
718
719
720
721
722
723
724
725
726
727

    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, 
728
					std::vector<DegreeOfFreedom> &newDOF)
729
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
730
731
732
    for (int i = first; i <= last; i++) {
      if (newDOF[i] >= 0) {
	vec[newDOF[i]] = vec[i];
733
734
735
736
737
738
739
740
741
      }
    }
  }

  template<typename T>
  ElementVector* DOFVectorBase<T>::assemble(T factor, ElInfo *elInfo,
					    const BoundaryType *bound, 
					    Operator *op)
  {
742
    FUNCNAME("DOFVector::assemble()");
743

744
745
    if (!(op || this->operators.size())) 
      return NULL;
746
747
748
749

    Operator *operat = op ? op : this->operators[0];

    this->elementVector = 
750
      operat->getAssembler(omp_get_thread_num())->initElementVector(this->elementVector, elInfo);
751

752
    if (op) {
753
754
      op->getElementVector(elInfo, this->elementVector);
    } else {
755
756
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
757
758
759
760
761
762
763
      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	   it != this->operators.end(); 
	   ++it, ++factorIt) {
	(*it)->getElementVector(elInfo, 
				this->elementVector, 
				*factorIt ? **factorIt : 1.0);
      }
764
765
766
    }

    addElementVector(factor,
767
768
    		     *this->elementVector, 
    		     bound);
769
770
771
772
773
774
775
776

    return this->elementVector;
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
777
    std::vector<Operator*>::iterator op;
778
779
780
781
782
783
784
    for(op = this->operators.begin(); op != this->operators.end(); ++op) {
      fillFlag |= (*op)->getFillFlag();
    }
    return fillFlag;
  }

  template<typename T>
785
786
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs )
  {
787
788
789
790
791
792
793
    feSpace = rhs.feSpace;
    vec = rhs.vec;
    this->elementVector = NULL;
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
794
    if (rhs.boundaryManager) {
795
796
      if (this->boundaryManager) 
	delete this->boundaryManager; 
797
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
798

799
      //    boundaryManager->setDOFVector(this);
800
    } else {
801
      this->boundaryManager=NULL;
802
803
    }

804
805
806
807
808
809
    return *this;
  }

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

812
813
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
814
815

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
816
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
817
818
      (*vecIterator) *= scal; 
    };
819

820
821
822
823
824
825
826
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
827
828
829
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
830
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
831
832
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
833
834
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
835
836
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
837
838
    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);
839
840
841
842
843
844
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {     
      *xIterator += *yIterator; 
    }

845
846
847
848
849
850
    return x;
  }

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

853
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
854
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
855
856
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
857
858
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
859
860
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
861
862
    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);
863
864
865
866
867
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {
      *xIterator -= *yIterator; 
    }
868
869
870
871
872
873
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
874
875
876
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
877
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
878
879
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
880
881
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
882
883
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
884
885
    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);
886
887
888
889
890
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {
      *xIterator *= *yIterator; 
    }
891
892
893
894
895
896
897
898
899

    return x;
  }

  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::operator*(DOFVector<T>& x, DOFVector<T>& y)");

900
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
901
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
902
903
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() && 
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
904
905
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
906
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
907

908
    T dot = 0;    
909
910
911
912
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
913
914
915
	 ++xIterator, ++yIterator) {
      dot += (*xIterator) * (*yIterator);
    };
916
917
918
919
920
921
922
923
924
925
926
927
928

    return(dot);
  }

  template<typename T>
  void mv(MatrixTranspose transpose,
	  const DOFMatrix &a, 
	  const DOFVector<T>&x,
	  DOFVector<T> &result,
	  bool add)
  {
    FUNCNAME("DOFVector<T>::mv()");

929
    TEST_EXIT_DBG(a.getRowFESpace() && a.getColFESpace() && x.getFESpace() && result.getFESpace())
930
931
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
932
933
934
935
936
937
938
939
    
    TEST_EXIT_DBG((a.getRowFESpace()->getAdmin() && a.getColFESpace()->getAdmin()) && 
	      (((transpose == NoTranspose) && 
		(a.getColFESpace()->getAdmin() == x.getFESpace()->getAdmin()) && 
		(a.getRowFESpace()->getAdmin() == result.getFESpace()->getAdmin())) ||
	       ((transpose == Transpose) && 
		(a.getRowFESpace()->getAdmin() == x.getFESpace()->getAdmin()) && 
		(a.getColFESpace()->getAdmin() == result.getFESpace()->getAdmin()))))
940
941
942
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());
943
    
944
    if (transpose == NoTranspose) {
945
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
946
	("x.size = %d too small: admin->sizeUsed = %d\n",
947
948
	 x.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( result.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
949
	("size = %d too small: admin->sizeUsed = %d\n",
950
951
	 result.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
952
	("a.size = %d too small: admin->sizeUsed = %d\n",
953
954
	 a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      
955

Thomas Witkowski's avatar
Thomas Witkowski committed
956
957
958
959
960
961
962
963
964
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS); 
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
      for(vecIterator.reset(), rowIterator.reset();
      	  !rowIterator.end();
      	  ++rowIterator, ++vecIterator) { 
      	
	double sum = 0;
      	if (!add) 
	  *vecIterator = 0.0;
965
      	for(std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
966
967
968
969
970
971
972
973
974
975
976
977
      	    colIterator != rowIterator->end();
      	    colIterator++) {
      	  
	  int jcol = colIterator->col;
      	  if (jcol >= 0) { // entry used? 
      	    sum += (static_cast<double>(colIterator->entry)) * x[jcol];
      	  } else {
      	    if (jcol == DOFMatrix::NO_MORE_ENTRIES)
      	      break;
      	  }
      	}
      	*vecIterator += sum;
978
979
980
      }

    } else if (transpose == Transpose) {
981
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
982
	("x.size = %d too small: admin->sizeUsed = %d\n",
983
984
	 x.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( result.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
985
	("size = %d too small: admin->sizeUsed = %d\n",
986
987
	 result.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
988
	("a.size = %d too small: admin->sizeUsed = %d\n",
989
990
	 a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
      if(!add) {
	typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
	for(vecIterator.reset();
	    !vecIterator.end(); ++vecIterator) {
	  *vecIterator = 0; 
	};
      }

      typename DOFVector<T>::Iterator xIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(&x)), USED_DOFS);
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
      for(xIterator.reset(), rowIterator.reset();
	  !rowIterator.end();
	  ++rowIterator, ++xIterator) { 
	T ax = (*xIterator);
1005
	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
	    colIterator != rowIterator->end();
	    colIterator++) {
	  int jcol = colIterator->col;
	  if (jcol >= 0) // entry used?
	    result[jcol] += (static_cast<double>(colIterator->entry)) * ax;
	  else 
	    if (jcol == DOFMatrix::NO_MORE_ENTRIES) break;
	}
      }
    } else {
      ERROR_EXIT("transpose = %d\n", transpose);
    }
  }

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

1027
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
1028
1029
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

1030
1031
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() && 
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
1032
1033
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
1034
    TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= x.getFESpace()->getAdmin()->getUsedSize())
1035
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
1036
1037
       x.getFESpace()->getAdmin()->getUsedSize());
    TEST_EXIT_DBG(static_cast<int>(y.getSize()) >= x.getFESpace()->getAdmin()->getUsedSize())
1038
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
1039
1040
       x.getFESpace()->getAdmin()->getUsedSize());
    
1041
1042
1043
1044
1045
1046
1047
    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); 
    };
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
  }

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

1079
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
1080
1081
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

1082
1083
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() && 
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
1084
1085
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
1086
    TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= x.getFESpace()->getAdmin()->getUsedSize())
1087
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
1088
1089
       x.getFESpace()->getAdmin()->getUsedSize());
    TEST_EXIT_DBG(static_cast<int>(y.getSize()) >= x.getFESpace()->getAdmin()->getUsedSize())
1090
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
1091
1092
       x.getFESpace()->getAdmin()->getUsedSize());
    
1093
1094
1095
1096
1097
1098
1099
    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); 
    };
1100
1101
1102
1103
1104
1105
1106
1107
1108
  }

  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);
1109
1110
1111
1112
1113
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end();
	 ++vIterator, ++rIterator) {
      *rIterator = scal * (*vIterator);
    };
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148

    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);
    for(vIterator.reset(), rIterator.reset();
	!vIterator.end();
	++vIterator, ++rIterator) 
      {  
	*rIterator = (*vIterator) + scal;
      };
    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);
    for(v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	!v1Iterator.end();
	++v1Iterator, ++v2Iterator, ++rIterator) 
      {  
	*rIterator = (*v1Iterator) + (*v2Iterator);
      };
    return result;
1149
1150
1151
1152
1153
1154
1155
1156
  }

  template<typename T>
  inline const DOFVector<T>& mod(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);
1157

1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end();
	 ++vIterator, ++rIterator) {
      *rIterator = fmod((*vIterator), 1.0);
    }

    return result;
  }

  template<typename T>
  inline const DOFVector<T>& Tanh(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);

    double eps;
    GET_PARAMETER(1, "vecphase->epsilon", "%f", &eps);
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end();
	 ++vIterator, ++rIterator) {
      *rIterator = 0.5 * (1.0 - tanh(3.0 * (*vIterator) / eps));
    }

    return result;
1183
1184
  }

1185

1186
1187
1188
1189
  template<typename T>
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
  {
    static T* localVec = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
1190
    T *result;   
1191
    if (d) {
1192
1193
      result = d;
    } else {
1194
1195
      if (localVec) 
	delete [] localVec;
1196
      localVec = new T[nBasFcts]; 
1197
      result = localVec;     
1198
    }
1199
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1200
1201
1202
    DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
    feSpace->getBasisFcts()->getLocalIndices(el, feSpace->getAdmin(), 
					     myLocalIndices);
1203

1204
    for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1205
      result[i] = (*this)[myLocalIndices[i]];
1206
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1207
       
1208
1209
1210
1211
    return result;
  }

  template<typename T>