DOFVector.hh 39.6 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"
Thomas Witkowski's avatar
Thomas Witkowski committed
21
#include "Traverse.h"
22
23
24
25

namespace AMDiS {

  template<typename T>
26
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
27
28
29
    : feSpace(f),
      name(n),
      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
48

    elementVector = NEW ElementVector(this->nBasFcts);
49
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
50
  
51
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
54
55
  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
56
      FREE_MEMORY(localVectors[i], T, this->nBasFcts);
57
      DELETE grdPhis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
58
      DELETE grdTmp[i];
59
      DELETE D2Phis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
60
61
62
63
    }
  }
  
  template<typename T>
64
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
65
    : DOFVectorBase<T>(f, n),
66
      coarsenOperation(NO_OPERATION)
67
68
69
70
71
  {
    init(f, n);
  } 

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

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

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

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

  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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
  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;

Thomas Witkowski's avatar
Thomas Witkowski committed
119
      if (!(condition && condition->isDirichlet())) {
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
122
123
124
125
126
	DegreeOfFreedom irow = elVec.dofIndices[i];
	(*this)[irow] = (add ? (*this)[irow] : 0.0);
	(*this)[irow] += factor * elVec[i];
      }
    }
  }

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

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

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

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

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

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

    return(nrm);
  }

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

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

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

    return(nrm);
  }

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

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


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


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

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

266
267
268
269
270
271
    return m;
  }

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

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

286
287
288
    return m;
  }

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

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

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

    ysize = y.getSize();

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), FREE_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
350
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) { 
351
      *vecIterator = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
352
    }
353
354
355
356

    if (transpose == NoTranspose) {
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS); 
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
357
358
      for (vecIterator.reset(), rowIterator.reset(); !rowIterator.end();
	   ++rowIterator, ++vecIterator) { 
359
360
	sum = 0;
	row = &(a[rowIterator.getDOFIndex()]);
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
363
	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
	     colIterator != rowIterator->end();
	     colIterator++) {
364
365
366
367
368
369
370
371
372
373
	  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;
Thomas Witkowski's avatar
Thomas Witkowski committed
374
      }
375
376
    } else if (transpose == Transpose) {
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
377
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
378
	*vecIterator  *= beta ; 
Thomas Witkowski's avatar
Thomas Witkowski committed
379
      }
380
381
382
    
      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);
Thomas Witkowski's avatar
Thomas Witkowski committed
383
384
      for (xIterator.reset(), rowIterator.reset(); !rowIterator.end();
	   ++rowIterator, ++xIterator) { 
385
386
	ax = alpha * (*xIterator);
	row = &(a[rowIterator.getDOFIndex()]);
Thomas Witkowski's avatar
Thomas Witkowski committed
387
	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
	    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
  {
406
    FUNCNAME("DOFVector<T>::print()");
407
    const DOFAdmin *admin = NULL;
408
    const char *format;
409

410
411
    if (feSpace) 
      admin = feSpace->getAdmin();
412
413

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

453
454
455
456
457
458
459
460
461
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
462
463
464
465
466
467
468
469
470

  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;

471
    for (int i = 0; i < numberOfBasFcts; i++)
472
473
474
475
476
477
478
479
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

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

482
483
484
    TEST_EXIT_DBG(fct)("no function to interpolate\n");

    interFct = fct;
485

486
487
488
489
490
    if (!this->getFESpace()) {
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
491

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

504
505
506
507
508
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
509
510
511

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
514
515
516
517
518
    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);
    }
519
520
521
522
523
524
525
  }

  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
526
527
528
529
530
531
532
    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);
533
534

    int number = basFct->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
535
    for (int i = 0; i < number; i++)
536
537
538
539
540
541
542
543
      (*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
544
    FUNCNAME("DOFVector::Int()");
545
546
547
  
    Mesh* mesh = feSpace->getMesh();

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

    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();
574
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
575
576

    int numPoints = quad_fast->getNumPoints();
577
578
579
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*(uh_vec[iq]);
    }
580
581
582
583
584
585
586
587
    norm += det*normT;

    return 0;
  }

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

    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
598
    Mesh* mesh = feSpace->getMesh();
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616

    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();
617
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
618
619

    int numPoints = quad_fast->getNumPoints();
620
621
622
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*abs(uh_vec[iq]);
    }
623
624
625
626
627
628
629
630
631
    norm += det*normT;

    return 0;
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
634
635
636
637
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
638
639
640
641
642
643
644

    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
645
646
    Mesh* mesh = feSpace->getMesh();
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, L2Norm_fct);
647
648
649
650
651
652
653
654
655
656
657
658

    return norm;  
  }

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

    det = elinfo->getDet();
659
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
660
661

    int numPoints = quad_fast->getNumPoints();
662
663
664
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq]);
    }
665
666
667
668
669
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
    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
700
    FUNCNAME("DOFVector::H1NormSquare()");
701

Thomas Witkowski's avatar
Thomas Witkowski committed
702
703
704
705
706
707
708
    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);
709
710

    norm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
711
712
    traverseVector = const_cast<DOFVector<T>*>(this); 
    Mesh *mesh = feSpace->getMesh();
713
714
715
716
717
718
719
720
721
722

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

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
736
    std::vector<Operator*>::iterator op;
Thomas Witkowski's avatar
Thomas Witkowski committed
737
    for (op = this->operators.begin(); op != this->operators.end(); ++op) {
738
739
740
741
742
743
      fillFlag |= (*op)->getFillFlag();
    }
    return fillFlag;
  }

  template<typename T>
744
745
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs )
  {
746
    feSpace = rhs.feSpace;
747
    DOFVectorBase<T>::feSpace = rhs.feSpace;
748
    this->nBasFcts = rhs.nBasFcts;
749
    vec = rhs.vec;
750
    this->elementVector = new ElementVector(this->nBasFcts);
751
752
753
754
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
755
    if (rhs.boundaryManager) {
756
757
      if (this->boundaryManager) 
	delete this->boundaryManager; 
758
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
759

760
      //    boundaryManager->setDOFVector(this);
761
    } else {
762
      this->boundaryManager = NULL;
763
764
    }

765
766
767
768
769
770
    return *this;
  }

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

773
774
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
775
776

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

781
782
783
784
785
786
787
    return x;
  }


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

806
807
808
809
810
811
    return x;
  }

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

814
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
815
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
816
817
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
818
819
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
820
821
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
822
823
    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);
824
825
826
827
828
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {
      *xIterator -= *yIterator; 
    }
829
830
831
832
833
834
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
835
836
837
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
838
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
839
840
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
841
842
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
843
844
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
845
846
    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);
847
848
849
850
851
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {
      *xIterator *= *yIterator; 
    }
852
853
854
855
856
857
858
859
860

    return x;
  }

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

861
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
862
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
863
864
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() && 
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
865
866
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
867
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
868

869
    T dot = 0;    
870
871
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
872
    for (xIterator.reset(), yIterator.reset(); !xIterator.end(); 
873
874
	 ++xIterator, ++yIterator) {
      dot += (*xIterator) * (*yIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
875
    }
876
877
878
879
880
881
882
883
884
885
886
887
888

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

889
    TEST_EXIT_DBG(a.getRowFESpace() && a.getColFESpace() && x.getFESpace() && result.getFESpace())
890
891
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
892
893
894
895
896
897
898
899
    
    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()))))
900
901
902
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());
903
    
904
    if (transpose == NoTranspose) {
905
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
906
	("x.size = %d too small: admin->sizeUsed = %d\n",
907
908
	 x.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( result.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
909
	("size = %d too small: admin->sizeUsed = %d\n",
910
911
	 result.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
912
	("a.size = %d too small: admin->sizeUsed = %d\n",
913
914
	 a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      
915

Thomas Witkowski's avatar
Thomas Witkowski committed
916
917
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS); 
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
918
      for (vecIterator.reset(), rowIterator.reset(); !rowIterator.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
919
920
921
922
923
      	  ++rowIterator, ++vecIterator) { 
      	
	double sum = 0;
      	if (!add) 
	  *vecIterator = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
924
925
926
      	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
	     colIterator != rowIterator->end();
	     colIterator++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
927
928
929
930
931
932
933
934
935
936
      	  
	  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;
937
938
939
      }

    } else if (transpose == Transpose) {
940
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
941
	("x.size = %d too small: admin->sizeUsed = %d\n",
942
943
	 x.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( result.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
944
	("size = %d too small: admin->sizeUsed = %d\n",
945
946
	 result.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
947
	("a.size = %d too small: admin->sizeUsed = %d\n",
948
949
	 a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      
Thomas Witkowski's avatar
Thomas Witkowski committed
950
      if (!add) {
951
	typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
952
	for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
953
	  *vecIterator = 0; 
Thomas Witkowski's avatar
Thomas Witkowski committed
954
	}
955
956
957
958
      }

      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);
Thomas Witkowski's avatar
Thomas Witkowski committed
959
960
      for (xIterator.reset(), rowIterator.reset(); !rowIterator.end();
	   ++rowIterator, ++xIterator) { 
961
	T ax = (*xIterator);
962
	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
	    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()");

984
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
985
986
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

987
988
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() && 
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
989
990
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
991
    TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= x.getFESpace()->getAdmin()->getUsedSize())
992
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
993
994
       x.getFESpace()->getAdmin()->getUsedSize());
    TEST_EXIT_DBG(static_cast<int>(y.getSize()) >= x.getFESpace()->getAdmin()->getUsedSize())
995
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
996
997
       x.getFESpace()->getAdmin()->getUsedSize());
    
998
999
    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);
Thomas Witkowski's avatar
Thomas Witkowski committed
1000
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
1001
1002
	 ++xIterator, ++yIterator) {
      *yIterator += alpha * (*xIterator); 
Thomas Witkowski's avatar
Thomas Witkowski committed
1003