DOFVector.hh 33.9 KB
Newer Older
1
#include <list>
Thomas Witkowski's avatar
Thomas Witkowski committed
2
#include <algorithm>
3
#include <math.h>
4
5
6
7
8
9
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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/utility/category.hpp>
#include <boost/numeric/linear_algebra/identity.hpp>


// Defining the interface for MTL4
namespace mtl {

  // Let MTL4 know that DOFVector it is a column vector
  namespace traits {
    template <typename T>
    struct category< AMDiS::DOFVector<T> > 
    {
      typedef tag::dense_col_vector type;
    };
  }

  namespace ashape {
    template <typename T>
    struct ashape< AMDiS::DOFVector<T> > 
    {
      typedef cvec<typename ashape<T>::type> type;
    };
  }

  // Modelling Collection and MutableCollection
  template <typename T>
  struct Collection< AMDiS::DOFVector<T> > 
  {
    typedef T               value_type;
    typedef const T&        const_reference;
    typedef std::size_t     size_type;
  };

  template <typename T>
  struct MutableCollection< AMDiS::DOFVector<T> > 
    : public Collection< AMDiS::DOFVector<T> > 
  {
    typedef T&              reference;
  };


} // namespace mtl



70
71
72
namespace AMDiS {

  template<typename T>
73
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
74
75
    : feSpace(f),
      name(n),
76
      elementVector(NULL),
Thomas Witkowski's avatar
Thomas Witkowski committed
77
      boundaryManager(NULL)
78
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
79
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
80
    dim = feSpace->getMesh()->getDim();
81

82
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
83
    localVectors.resize(omp_get_overall_max_threads());
84
    grdPhis.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
85
    grdTmp.resize(omp_get_overall_max_threads());
86
    D2Phis.resize(omp_get_overall_max_threads());
87

88
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
89
90
      localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
      localVectors[i] = GET_MEMORY(T, this->nBasFcts);
91
      grdPhis[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
92
      grdTmp[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
93
      D2Phis[i] = NEW DimMat<double>(dim, NO_INIT);
94
    }
95
96

    elementVector = NEW ElementVector(this->nBasFcts);
97
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
98
  
99
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
100
101
102
103
  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
104
      FREE_MEMORY(localVectors[i], T, this->nBasFcts);
105
      DELETE grdPhis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
106
      DELETE grdTmp[i];
107
      DELETE D2Phis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
108
109
110
111
    }
  }
  
  template<typename T>
112
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
113
    : DOFVectorBase<T>(f, n),
114
      //elementVector(NULL), 
115
      coarsenOperation(NO_OPERATION)
116
117
  {
    init(f, n);
118
119
120
121

#ifdef HAVE_PARALLEL_AMDIS
    applicationOrdering = NULL;
#endif
122
123
124
  } 

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
141
    if (this->boundaryManager)
142
      DELETE this->boundaryManager;
143
144

    vec.clear();
145
146
147
148
149
150
151
152
153
154
155
  }

  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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
  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
170
      if (!(condition && condition->isDirichlet())) {
Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
173
174
175
176
177
	DegreeOfFreedom irow = elVec.dofIndices[i];
	(*this)[irow] = (add ? (*this)[irow] : 0.0);
	(*this)[irow] += factor * elVec[i];
      }
    }
  }

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

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

    return(sqrt(nrm));
  }

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

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

    return nrm;
  }

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

221
222
223
224
225
    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());
226

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

    return(nrm);
  }

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

240
241
242
243
244
    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());
245

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

    return(nrm);
  }

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

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


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


  template<typename T>
  T DOFVector<T>::min() const
  {
303
304
305
306
307
    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())
308
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
309
310
311
       feSpace->getAdmin()->getUsedSize());

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

317
318
319
320
321
322
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
323
324
325
326
327
    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())
328
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
329
       feSpace->getAdmin()->getUsedSize());
330

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

337
338
339
    return m;
  }

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

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

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

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

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

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

  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;

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

    return val;
  }

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

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

    interFct = fct;
428

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

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

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

    traverseVector = this;

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

  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
469
470
471
472
473
474
475
    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);
476
477

    int number = basFct->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
478
    for (int i = 0; i < number; i++)
479
480
481
482
483
484
485
486
      (*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
487
    FUNCNAME("DOFVector::Int()");
488
489
490
  
    Mesh* mesh = feSpace->getMesh();

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

    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();
517
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
518
519

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

    return 0;
  }

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

    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
541
    Mesh* mesh = feSpace->getMesh();
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559

    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();
560
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
561
562

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

    return 0;
  }


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

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

    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
588
589
    Mesh* mesh = feSpace->getMesh();
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, L2Norm_fct);
590
591
592
593
594
595
596
597
598
599
600
601

    return norm;  
  }

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

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

    int numPoints = quad_fast->getNumPoints();
605
606
607
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq]);
    }
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
    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
643
    FUNCNAME("DOFVector::H1NormSquare()");
644

Thomas Witkowski's avatar
Thomas Witkowski committed
645
646
647
648
649
650
651
    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);
652
653

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

    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, 
666
					std::vector<DegreeOfFreedom> &newDOF)
667
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
668
669
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
670
	vec[newDOF[i]] = vec[i];
671
672
673
674
675
676
  }

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

683
684
685
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
686
687
688
689
690
691
692
693
694
695
  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();
    }
  }

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

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

718
719
720
721
722
723
    return *this;
  }

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

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

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

734
735
736
737
738
739
740
    return x;
  }


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

759
760
761
762
763
764
    return x;
  }

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

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

    return x;
  }

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

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

    dot = 0;
824
825
826

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

    return(dot);
  }

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

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

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

    const DOFAdmin *rowAdmin = a.getRowFESpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFESpace()->getAdmin();

    TEST_EXIT((rowAdmin && colAdmin) && 
	      (((transpose == NoTranspose) && (colAdmin == x.getFESpace()->getAdmin()) && 
		(rowAdmin == result.getFESpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFESpace()->getAdmin()) && 
		(colAdmin == result.getFESpace()->getAdmin()))))
861
862
863
864
865
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());


866
867
868
869
870
    if (transpose == NoTranspose)
      mult(a.getBaseMatrix(), x, result);
    else if (transpose == Transpose)
      mult(trans(const_cast<DOFMatrix::base_matrix_type&>(a.getBaseMatrix())), x, result);
    else 
871
872
873
874
875
876
877
878
879
880
      ERROR_EXIT("transpose = %d\n", transpose);
  }

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

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

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

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

      // This is the old implementation of the mv-multiplication. It have been changed
      // because of the OpenMP-parallelization:
//     typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
//     typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
//     for(xIterator.reset(), yIterator.reset();
// 	!xIterator.end();
// 	++xIterator, ++yIterator) 
//       {  
// 	*yIterator += alpha * (*xIterator); 
//       };

      int i;
      int maxI = y.getSize();

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
      for (i = 0; i < maxI; i++) {
	if (!admin->isDOFFree(i)) {
	  y[i] += alpha * x[i];
	}
      }
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
  }

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

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

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

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

    // This is the old implementation of the mv-multiplication. It have been changed
    // because of the OpenMP-parallelization:
    //     typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
    //     typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
    //     for(xIterator.reset(), yIterator.reset();
    // 	!xIterator.end();
    // 	++xIterator, ++yIterator) 
    //       {  
    // 	*yIterator = alpha *(*yIterator)+ (*xIterator); 
    //       };

      int i;
      int maxI = y.getSize();

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
      for (i = 0; i < maxI; i++) {
	if (!admin->isDOFFree(i)) {
	  y[i] = alpha * y[i] + x[i];
	}
      }
986
987
988
989
990
991
992
993
994
  }

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

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

1037
1038
1039
  }

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

1045
1046
    int i;
    int nBasFcts = feSpace->getBasisFcts()->getNumber();
1047

1048
1049
1050
1051
1052
1053
1054
1055
1056
    T *result;
  
    if(d) {
      result = d;
    } else {
      if(localVec) delete [] localVec;
      localVec = new T[nBasFcts]; 
      result = localVec;
    }
1057

1058
1059
    const DegreeOfFreedom *localIndices = 
      feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);
1060

1061
1062
    for(i = 0; i < nBasFcts; i++) {
      result[i] = (*this)[localIndices[i]];
1063
1064
1065
    }

    return result;
1066
1067
1068
  }

  template<typename T>
1069
1070
  const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo         *elInfo, 
					 const Quadrature     *quad,
1071
					 const FastQuadrature *quadFast,
1072
					 T                    *vecAtQPs) const
1073
1074
1075
  {
    FUNCNAME("DOFVector<T>::getVecAtQPs()");
  
1076
    TEST_EXIT(quad || quadFast)("neither quad nor quadFast defined\n");
1077

1078
1079
1080
    if(quad && quadFast) {
      TEST_EXIT(quad == quadFast->getQuadrature())
	("quad != quadFast->quadrature\n");
1081
1082
    }

1083
    TEST_EXIT(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
1084
1085
      ("invalid basis functions");

1086
1087
1088
1089
    Element *el = elInfo->getElement();

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

1090
    const BasisFunction *basFcts = feSpace->getBasisFcts();
1091
1092
1093
1094
1095

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

1096
    static T *localvec = NULL;
1097

1098
1099
1100
1101
1102
    T *result;

    if (vecAtQPs) {
      result = vecAtQPs;
    } else {
1103
1104
1105
      if(localvec) delete [] localvec;
      localvec = new T[numPoints];
      for(i = 0; i < numPoints; i++) {