DOFVector.hh 39.8 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
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
727
	vec[newDOF[i]] = vec[i];
728
729
730
731
732
733
  }

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

740
741
742
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
745
746
747
748
749
750
751
752
  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();
    }
  }

753
  template<typename T>
754
755
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs )
  {
756
    feSpace = rhs.feSpace;
757
    DOFVectorBase<T>::feSpace = rhs.feSpace;
758
    this->nBasFcts = rhs.nBasFcts;
759
    vec = rhs.vec;
760
    this->elementVector = new ElementVector(this->nBasFcts);
761
762
763
764
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
765
    if (rhs.boundaryManager) {
766
767
      if (this->boundaryManager) 
	delete this->boundaryManager; 
768
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
769

770
      //    boundaryManager->setDOFVector(this);
771
    } else {
772
      this->boundaryManager = NULL;
773
774
    }

775
776
777
778
779
780
    return *this;
  }

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

783
784
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
785
786

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

791
792
793
794
795
796
797
    return x;
  }


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

816
817
818
819
820
821
    return x;
  }

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

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

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
845
846
847
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
848
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
849
850
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
851
852
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
853
854
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
855
856
    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);
857
858
859
860
861
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {
      *xIterator *= *yIterator; 
    }
862
863
864
865
866
867
868
869
870

    return x;
  }

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

871
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
872
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
873
874
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() && 
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
875
876
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
877
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
878

879
    T dot = 0;    
880
881
    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
882
    for (xIterator.reset(), yIterator.reset(); !xIterator.end(); 
883
884
	 ++xIterator, ++yIterator) {
      dot += (*xIterator) * (*yIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
885
    }
886
887
888
889
890
891
892
893
894
895
896
897
898

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

899
    TEST_EXIT_DBG(a.getRowFESpace() && a.getColFESpace() && x.getFESpace() && result.getFESpace())
900
901
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
902
903
904
905
906
907
908
909
    
    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()))))
910
911
912
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());
913
    
914
    if (transpose == NoTranspose) {
915
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
916
	("x.size = %d too small: admin->sizeUsed = %d\n",
917
918
	 x.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( result.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
919
	("size = %d too small: admin->sizeUsed = %d\n",
920
921
	 result.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
922
	("a.size = %d too small: admin->sizeUsed = %d\n",
923
924
	 a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      
925

Thomas Witkowski's avatar
Thomas Witkowski committed
926
927
      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
928
      for (vecIterator.reset(), rowIterator.reset(); !rowIterator.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
929
930
931
932
933
      	  ++rowIterator, ++vecIterator) { 
      	
	double sum = 0;
      	if (!add) 
	  *vecIterator = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
934
935
936
      	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
	     colIterator != rowIterator->end();
	     colIterator++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
937
938
939
940
941
942
943
944
945
946
      	  
	  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;
947
948
949
      }

    } else if (transpose == Transpose) {
950
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
951
	("x.size = %d too small: admin->sizeUsed = %d\n",
952
953
	 x.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( result.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
954
	("size = %d too small: admin->sizeUsed = %d\n",
955
956
	 result.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
957
	("a.size = %d too small: admin->sizeUsed = %d\n",
958
959
	 a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      
Thomas Witkowski's avatar
Thomas Witkowski committed
960
      if (!add) {
961
	typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
962
	for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
963
	  *vecIterator = 0; 
Thomas Witkowski's avatar
Thomas Witkowski committed
964
	}
965
966
967
968
      }

      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
969
970
      for (xIterator.reset(), rowIterator.reset(); !rowIterator.end();
	   ++rowIterator, ++xIterator) { 
971
	T ax = (*xIterator);
972
	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
	    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()");

994
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
995
996
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

997
998
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() && 
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
999
1000
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
1001
    TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= x.getFESpace()->getAdmin()->getUsedSize())
1002
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
1003
1004
       x.getFESpace()->getAdmin()->getUsedSize());
    TEST_EXIT_DBG(static_cast<int>(y.getSize()) >= x.getFESpace()->getAdmin()->getUsedSize())
1005
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
1006
1007
       x.getFESpace()-&g