DOFVector.hh 29.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include <list>
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include <algorithm>
15
#include <math.h>
16

17
18
19
20
21
#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>

22
23
24
25
26
27
28
29
30
31
32
33
#include "FixVec.h"
#include "Boundary.h"
#include "DOFAdmin.h"
#include "ElInfo.h"
#include "Error.h"
#include "FiniteElemSpace.h"
#include "Global.h"
#include "Mesh.h"
#include "Quadrature.h"
#include "AbstractFunction.h"
#include "BoundaryManager.h"
#include "Assembler.h"
34
#include "Operator.h"
35
#include "Initfile.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
36
#include "Traverse.h"
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
70
71
72
73
74
75
76
77
78
// 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



79
80
81
namespace AMDiS {

  template<typename T>
82
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
83
84
    : feSpace(f),
      name(n),
85
      elementVector(f->getBasisFcts()->getNumber()),
Thomas Witkowski's avatar
Thomas Witkowski committed
86
      boundaryManager(NULL)
87
  {    
Thomas Witkowski's avatar
Thomas Witkowski committed
88
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
89
    dim = feSpace->getMesh()->getDim();
90
91
92
93
94
  }
  

  template<typename T>
  DOFVectorBase<T>::~DOFVectorBase()
95
  {}
96

97
 
Thomas Witkowski's avatar
Thomas Witkowski committed
98
  template<typename T>
99
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
100
    : DOFVectorBase<T>(f, n),
101
      coarsenOperation(COARSE_INTERPOL)
102
  {
103
    vec.resize(0);
104
105
106
    init(f, n);
  } 

107

108
  template<typename T>
109
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
110
  {
111
112
113
114
115
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
    this->boundaryManager = new BoundaryManager(f);
116
117
118
119
120
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
121
122
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
123

Thomas Witkowski's avatar
Thomas Witkowski committed
124
    if (this->boundaryManager)
125
      delete this->boundaryManager;
126
127

    vec.clear();
128
129
130
131
132
  }

  template<typename T>
  DOFVector<T> * DOFVector<T>::traverseVector = NULL;

Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
135
136
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
137
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
138
139
140
141
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

142
143
144
    std::vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
145

146
    for (DegreeOfFreedom i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
147
148
149
      BoundaryCondition *condition = 
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;

Thomas Witkowski's avatar
Thomas Witkowski committed
150
      if (!(condition && condition->isDirichlet())) {
151
	DegreeOfFreedom irow = indices[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
152

153
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
154
	  (*this)[irow] += factor * elVec[i];
155
156
	else  
	  (*this)[irow] = factor * elVec[i];	
Thomas Witkowski's avatar
Thomas Witkowski committed
157
158
159
160
      }
    }
  }

161
162
163
164
165
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

166
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
167

168
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
169
170
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
171
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
172
173
174
175
176
177
178
179
180
181
      nrm += (*vecIterator) * (*vecIterator);

    return(sqrt(nrm));
  }

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

182
    checkFeSpace(this->feSpace, vec);
183
184
    
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
185
186
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
187
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
188
189
190
191
192
193
194
195
      nrm += (*vecIterator) * (*vecIterator);

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
196
    FUNCNAME("DOFVector<T>::asum()");
197

198
    checkFeSpace(this->feSpace, vec);
199

200
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
201
202
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
203
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
204
205
206
207
208
209
210
211
      nrm += abs(*vecIterator);

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
212
    FUNCNAME("DOFVector<T>::sum()");
213

214
    checkFeSpace(this->feSpace, vec);
215

216
    double nrm = 0.0;    
Thomas Witkowski's avatar
Thomas Witkowski committed
217
218
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
219
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
220
221
222
223
224
225
226
227
      nrm += *vecIterator;

    return(nrm);
  }

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

230
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
231

232
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
233
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
234
235
236
237
238
239
240
      *vecIterator = alpha ; 
  }


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
241
    FUNCNAME("DOFVector<T>::copy()");
Thomas Witkowski's avatar
Thomas Witkowski committed
242

243
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
244

245
246
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >= 
		  this->feSpace->getAdmin()->getUsedSize())
247
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
248
       this->feSpace->getAdmin()->getUsedSize());
249
    
250
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
251
252
253
254
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), 
		       USED_DOFS);
    for (vecIterator.reset(), xIterator.reset(); !vecIterator.end();
	 ++vecIterator, ++xIterator)      
255
      *vecIterator = *xIterator; 
256
257
258
259
260
261
  }


  template<typename T>
  T DOFVector<T>::min() const
  {
262
263
    FUNCNAME("DOFVector<T>::min()");
    
264
    checkFeSpace(this->feSpace, vec);
265
266

    T m;    
267
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
268
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
269
      m = std::min(m, *vecIterator);
270

271
272
273
    return m;
  }

274

275
276
277
  template<typename T>
  T DOFVector<T>::max() const 
  {
278
279
    FUNCNAME("DOFVector<T>::max()");
    
280
    checkFeSpace(this->feSpace, vec);
281

282
    T m;    
283
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
284
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
285
      m = std::max(m, *vecIterator);
286

287
    return m;
288

289
290
  }

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

297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316

  template<typename T>
  T DOFVector<T>::average() const
  {
    FUNCNAME("DOFVector<T>::average()");
    
    checkFeSpace(this->feSpace, vec);

    int count = 0;
    T m = 0;
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
      m += *vecIterator;
      count++;
    }

    return m / count;
  }


317
318
319
  template<typename T>
  void DOFVector<T>::print() const
  {
320
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
321

322
    const DOFAdmin *admin = NULL;
323
    const char *format;
324

325
326
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
327

328
    MSG("Vec `%s':\n", this->name.c_str());
329
    int j = 0;
330
331
332
333
334
335
336
337
338
    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
339
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
340
	if ((j % 3) == 0) {
341
342
	  if (j) 
	    Msg::print("\n");
343
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
344
	} else {
345
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
346
	}
347
	j++;
348
      }
349
      Msg::print("\n");
350
    } else {
351
352
	MSG("no DOFAdmin, print whole vector.\n");
    
353
354
355
356
357
358
	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 {
359
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
360
	  }
361
362
363
364
365
366
367
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }

368
369
370
371
372
373
374
375
376
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
377
378
379
380
381

  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda, 
			     DegreeOfFreedom* dof_indices)
  {
382
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFeSpace()->getBasisFcts());
383
    int nBasisFcts = phi->getNumber();
384
385
    T val = 0.0;

386
    for (int i = 0; i < nBasisFcts; i++)
387
388
389
390
391
392
393
394
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
395
    FUNCNAME("DOFVector::interpol()");
396

Thomas Witkowski's avatar
Thomas Witkowski committed
397
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
398
399

    interFct = fct;
400

401
    if (!this->getFeSpace()) {
402
403
404
405
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
406

407
    if (!(this->getFeSpace()->getAdmin())) {
408
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
409
	  this->getFeSpace()->getName().c_str());
410
411
      return;
    }
412
  
413
    if (!(this->getFeSpace()->getBasisFcts())) {
414
415
416
417
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
418

419
420
421
422
423
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
424
425
426

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
427
    TraverseStack stack;
428
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
429
430
431
432
433
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
434
435
436
  }

  template<typename T>
437
  void DOFVector<T>::interpolFct(ElInfo* elinfo)
438
  {
439
440
    const BasisFunction *basFct = traverseVector->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFeSpace()->getAdmin();
441
442
443
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
444

445
    int nBasFcts = basFct->getNumber();
446
447
448
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), 
			    admin, myLocalIndices);
449
    for (int i = 0; i < nBasFcts; i++)
450
      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
451
452
453
454
455
  }

  template<typename T>
  double DOFVector<T>::Int(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
456
    FUNCNAME("DOFVector::Int()");
457
  
458
    Mesh* mesh = this->feSpace->getMesh();
459

Thomas Witkowski's avatar
Thomas Witkowski committed
460
    if (!q) {
461
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
462
463
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
464

465
    FastQuadrature *quadFast = 
466
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
467

468
469
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
470
    mtl::dense_vector<T> uh_vec(nPoints);
471
472
473
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
474
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
475
476
477
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
478
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
479
480
481
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;
482

483
484
      elInfo = stack.traverseNext(elInfo);
    }
485

486
    return result;  
487
488
  }

489

490
491
492
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
493
    FUNCNAME("DOFVector::L1Norm()");
494
  
495
496
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
497
    if (!q) {
498
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
499
500
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
501

502
503
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
504

505
506
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
507
    mtl::dense_vector<T> uh_vec(nPoints);
508
509
510
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
511
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
512
513
514
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
515
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
516
517
518
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
519

520
521
      elInfo = stack.traverseNext(elInfo);
    }
522

523
    return result;  
524
525
  }

526

527
528
529
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
530
    FUNCNAME("DOFVector::L2NormSquare()");
531

532
533
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
534
    if (!q) {
535
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
536
537
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
538

539
540
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
541

542
543
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
544
    mtl::dense_vector<T> uh_vec(nPoints);
545
546
547
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
548
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
549
550
551
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
552
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
553
554
555
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;
556

557
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
558
    }
559

560
    return result;  
561
562
  }

563

564
  template<typename T>  
565
566
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
567
    FUNCNAME("DOFVector::H1NormSquare()");
568

569
570
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
571
    if (!q) {
572
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
573
574
575
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

576
577
578
579
580
581
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
582
    std::vector<WorldVector<T> > grduh_vec(nPoints);
583
584
585
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
586
587
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
588
589
590
591
592
593
594
595
596
597
598
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getGrdAtQPs(elInfo, NULL, quadFast, &(grduh_vec[0]));

      for (int iq = 0; iq < nPoints; iq++) {
	double norm2 = 0.0;
	for (int j = 0; j < dimOfWorld; j++)
	  norm2 += sqr(grduh_vec[iq][j]);
	normT += quadFast->getWeight(iq) * norm2;
      }
599

600
      result += det * normT;
601

602
603
      elInfo = stack.traverseNext(elInfo);
    }
604

605
    return result;  
606
607
608
609
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
610
					std::vector<DegreeOfFreedom> &newDOF)
611
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
612
613
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
614
	vec[newDOF[i]] = vec[i];
615
616
617
618
619
620
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
621
622
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
623
	 op != this->operators.end(); ++op)
624
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
625

626
627
628
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
631
632
633
  template<typename T>
  void DOFVectorBase<T>::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = this->operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
634
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
637
      (*it)->finishAssembling();
  }

638
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
639
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
640
  {
641
    this->feSpace = rhs.feSpace;
642
    this->dim = rhs.dim;
643
    this->nBasFcts = rhs.nBasFcts;
644
    vec = rhs.vec;
645
    this->elementVector.change_dim(this->nBasFcts);
646
647
648
649
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
650
    
651
    if (rhs.boundaryManager) {
652
653
654
      if (this->boundaryManager) 
	delete this->boundaryManager; 

655
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
656
    } else {
657
      this->boundaryManager = NULL;
658
659
    }

660
661
662
663
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

664
665
666
667
668
669
    return *this;
  }

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

672
673
    TEST_EXIT_DBG(x.getFeSpace() && x.getFeSpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFeSpace(), x.getFeSpace()->getAdmin());
674

Thomas Witkowski's avatar
Thomas Witkowski committed
675
676
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
677
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
678
      (*vecIterator) *= scal; 
679

680
681
682
683
684
685
686
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
687
688
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
689
690
691
692
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
693
      ("no admin or different admins: %8X, %8X\n",
694
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
695
696
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
697
698
    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);
Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
701
702
      *xIterator += *yIterator; 

703
704
705
706
707
708
    return x;
  }

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

711
712
713
714
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
715
      ("no admin or different admins: %8X, %8X\n",
716
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
717
718
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
719
720
    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);
Thomas Witkowski's avatar
Thomas Witkowski committed
721
722
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
723
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
724

725
726
727
728
729
730
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
731
732
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
733
734
735
736
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
737
      ("no admin or different admins: %8X, %8X\n",
738
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
739
740
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
741
742
    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);
Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
745
      *xIterator *= *yIterator; 
746
747
748
749
750
751
752
753

    return x;
  }

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

756
757
758
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT((admin = x.getFeSpace()->getAdmin()) && (admin == y.getFeSpace()->getAdmin()))
759
      ("no admin or different admins: %8X, %8X\n",
760
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
761
762
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
763
    T dot = 0;
764
765
766

    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
767
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
768
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
769
      dot += (*xIterator) * (*yIterator);
770

Thomas Witkowski's avatar
Thomas Witkowski committed
771
    return dot;
772
773
774
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
777
  {
778
779
780
    // Unfortunately needed
    using namespace mtl;

781
782
    FUNCNAME("DOFVector<T>::mv()");

783
784
785
786
    TEST_EXIT(a.getRowFeSpace() && a.getColFeSpace() && 
	      x.getFeSpace() && result.getFeSpace())
      ("getFeSpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFeSpace(), a.getColFeSpace(), x.getFeSpace(), result.getFeSpace());
787

788
789
    const DOFAdmin *rowAdmin = a.getRowFeSpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFeSpace()->getAdmin();
790
791

    TEST_EXIT((rowAdmin && colAdmin) && 
792
793
794
795
	      (((transpose == NoTranspose) && (colAdmin == x.getFeSpace()->getAdmin()) && 
		(rowAdmin == result.getFeSpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFeSpace()->getAdmin()) && 
		(colAdmin == result.getFeSpace()->getAdmin()))))
796
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
797
798
       a.getRowFeSpace()->getAdmin(), a.getColFeSpace()->getAdmin(), 
       x.getFeSpace()->getAdmin(), result.getFeSpace()->getAdmin());
799
800


801
802
803
804
805
    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 
806
807
808
809
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
810
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
811
812
813
  {
    FUNCNAME("DOFVector<T>::axpy()");

814
815
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
816

817
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
818

819
    TEST_EXIT((admin) && (admin == y.getFeSpace()->getAdmin()))
820
      ("no admin or different admins: %8X, %8X\n",
821
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
822
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
823
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
824
825
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
826
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
827
828
       admin->getUsedSize());

829
830
831
    for (int i = 0; i < y.getSize(); i++)
      if (!admin->isDofFree(i))
	y[i] += alpha * x[i];
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
  }

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

863
864
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
865

866
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
867

868
    TEST_EXIT(admin && (admin == y.getFeSpace()->getAdmin()))
869
      ("no admin or different admins: %8X, %8X\n",
870
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
871
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
872
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
873
874
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())