DOFVector.hh 32.1 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 "OpenMP.h"
35
#include "Operator.h"
36
#include "Parameters.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
37
#include "Traverse.h"
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
79
// 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



80
81
82
namespace AMDiS {

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

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
    createTempData();
  }
  

  template<typename T>
  DOFVectorBase<T>::~DOFVectorBase()
  {
    for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
      delete [] localIndices[i];
      delete grdPhis[i];
      delete grdTmp[i];
      delete D2Phis[i];
    }
  }


  template<typename T>
  void DOFVectorBase<T>::createTempData()
  {
111
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
112
    localVectors.resize(omp_get_overall_max_threads());
113
    grdPhis.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
114
    grdTmp.resize(omp_get_overall_max_threads());
115
    D2Phis.resize(omp_get_overall_max_threads());
116

117
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
118
      localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
119
      localVectors[i].change_dim(this->nBasFcts);
120
121
122
      grdPhis[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      grdTmp[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      D2Phis[i] = new DimMat<double>(dim, NO_INIT);
123
    }    
124
  }
125

Thomas Witkowski's avatar
Thomas Witkowski committed
126
127
  
  template<typename T>
128
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
129
    : DOFVectorBase<T>(f, n),
130
      coarsenOperation(NO_OPERATION)
131
  {
132
    vec.resize(0);
133
134
135
    init(f, n);
  } 

136

137
  template<typename T>
138
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
139
  {
140
141
142
143
144
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
    this->boundaryManager = new BoundaryManager(f);
145
146
147
148
149
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
150
151
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
152

Thomas Witkowski's avatar
Thomas Witkowski committed
153
    if (this->boundaryManager)
154
      delete this->boundaryManager;
155
156

    vec.clear();
157
158
159
160
161
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
162
163
164
165
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
166
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
169
170
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

171
172
173
    std::vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
174

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

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

182
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
183
	  (*this)[irow] += factor * elVec[i];
184
185
	else  
	  (*this)[irow] = factor * elVec[i];	
Thomas Witkowski's avatar
Thomas Witkowski committed
186
187
188
189
      }
    }
  }

190
191
192
193
194
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

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

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

    return(sqrt(nrm));
  }

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

211
    checkFeSpace(this->feSpace, vec);
212
213
    
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
216
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
217
218
219
220
221
222
223
224
      nrm += (*vecIterator) * (*vecIterator);

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
225
    FUNCNAME("DOFVector<T>::asum()");
226

227
    checkFeSpace(this->feSpace, vec);
228

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

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
241
    FUNCNAME("DOFVector<T>::sum()");
242

243
    checkFeSpace(this->feSpace, vec);
244

245
    double nrm = 0.0;    
Thomas Witkowski's avatar
Thomas Witkowski committed
246
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
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
260

261
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
262
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
263
264
265
266
267
268
269
      *vecIterator = alpha ; 
  }


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

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

274
275
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >= 
		  this->feSpace->getAdmin()->getUsedSize())
276
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
277
       this->feSpace->getAdmin()->getUsedSize());
278
    
279
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
282
283
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), 
		       USED_DOFS);
    for (vecIterator.reset(), xIterator.reset(); !vecIterator.end();
	 ++vecIterator, ++xIterator)      
284
      *vecIterator = *xIterator; 
285
286
287
288
289
290
  }


  template<typename T>
  T DOFVector<T>::min() const
  {
291
292
    FUNCNAME("DOFVector<T>::min()");
    
293
    checkFeSpace(this->feSpace, vec);
294
295

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

300
301
302
    return m;
  }

303

304
305
306
  template<typename T>
  T DOFVector<T>::max() const 
  {
307
308
    FUNCNAME("DOFVector<T>::max()");
    
309
    checkFeSpace(this->feSpace, vec);
310

311
    T m;    
312
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
313
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
314
      m = std::max(m, *vecIterator);
315

316
    return m;
317

318
319
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
320
321
322
323
324
325
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345

  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;
  }


346
347
348
  template<typename T>
  void DOFVector<T>::print() const
  {
349
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
350

351
    const DOFAdmin *admin = NULL;
352
    const char *format;
353

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

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

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

    return result;
  }
406
407
408
409
410

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

415
    for (int i = 0; i < nBasisFcts; i++)
416
417
418
419
420
421
422
423
      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
424
    FUNCNAME("DOFVector::interpol()");
425

Thomas Witkowski's avatar
Thomas Witkowski committed
426
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
427
428

    interFct = fct;
429

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

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

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

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
456
    TraverseStack stack;
457
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
460
461
462
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
463
464
465
  }

  template<typename T>
466
  void DOFVector<T>::interpolFct(ElInfo* elinfo)
467
  {
468
469
    const BasisFunction *basFct = traverseVector->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFeSpace()->getAdmin();
470
471
472
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
473
474
475

    DegreeOfFreedom *myLocalIndices = this->localIndices[omp_get_thread_num()];
    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), admin, myLocalIndices);
476
477
    int nBasFcts = basFct->getNumber();
    for (int i = 0; i < nBasFcts; i++)
478
      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
479
480
481
482
483
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
488
    if (!q) {
489
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
490
491
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
492

493
    FastQuadrature *quadFast = 
494
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
495

496
497
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
498
    mtl::dense_vector<T> uh_vec(nPoints);
499
500
501
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
502
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
503
504
505
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
506
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
507
508
509
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;
510

511
512
      elInfo = stack.traverseNext(elInfo);
    }
513

514
    return result;  
515
516
  }

517

518
519
520
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
521
    FUNCNAME("DOFVector::L1Norm()");
522
  
523
524
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
525
    if (!q) {
526
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
527
528
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
529

530
531
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
532

533
534
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
535
    mtl::dense_vector<T> uh_vec(nPoints);
536
537
538
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
539
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
540
541
542
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
543
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
544
545
546
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
547

548
549
      elInfo = stack.traverseNext(elInfo);
    }
550

551
    return result;  
552
553
  }

554

555
556
557
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
558
    FUNCNAME("DOFVector::L2NormSquare()");
559

560
561
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
562
    if (!q) {
563
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
564
565
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
566

567
568
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
569

570
571
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
572
    mtl::dense_vector<T> uh_vec(nPoints);
573
574
575
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
576
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
577
578
579
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
580
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
581
582
583
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;
584

585
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
586
    }
587

588
    return result;  
589
590
  }

591

592
  template<typename T>  
593
594
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
595
    FUNCNAME("DOFVector::H1NormSquare()");
596

597
598
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
599
    if (!q) {
600
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
601
602
603
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

604
605
606
607
608
609
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
610
    std::vector<WorldVector<T> > grduh_vec(nPoints);
611
612
613
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
614
615
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
616
617
618
619
620
621
622
623
624
625
626
    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;
      }
627

628
      result += det * normT;
629

630
631
      elInfo = stack.traverseNext(elInfo);
    }
632

633
    return result;  
634
635
636
637
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
638
					std::vector<DegreeOfFreedom> &newDOF)
639
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
640
641
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
642
	vec[newDOF[i]] = vec[i];
643
644
645
646
647
648
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
649
650
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
651
	 op != this->operators.end(); ++op)
652
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
653

654
655
656
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
659
660
661
  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
662
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
663
664
665
      (*it)->finishAssembling();
  }

666
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
667
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
668
  {
669
    this->feSpace = rhs.feSpace;
670
    this->dim = rhs.dim;
671
    this->nBasFcts = rhs.nBasFcts;
672
    vec = rhs.vec;
673
    this->elementVector.change_dim(this->nBasFcts);
674
675
676
677
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
678
    
679
    if (rhs.boundaryManager) {
680
681
682
      if (this->boundaryManager) 
	delete this->boundaryManager; 

683
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
684
    } else {
685
      this->boundaryManager = NULL;
686
687
    }

688
689
690
691
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

692
693
694
695
696
697
    return *this;
  }

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

700
701
    TEST_EXIT_DBG(x.getFeSpace() && x.getFeSpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFeSpace(), x.getFeSpace()->getAdmin());
702

Thomas Witkowski's avatar
Thomas Witkowski committed
703
704
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
705
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
706
      (*vecIterator) *= scal; 
707

708
709
710
711
712
713
714
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
715
716
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
717
718
719
720
    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()))
721
      ("no admin or different admins: %8X, %8X\n",
722
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
723
724
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
725
726
    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
727
728
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
729
730
      *xIterator += *yIterator; 

731
732
733
734
735
736
    return x;
  }

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

739
740
741
742
    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()))
743
      ("no admin or different admins: %8X, %8X\n",
744
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
745
746
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
747
748
    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
749
750
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
751
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
752

753
754
755
756
757
758
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
759
760
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
761
762
763
764
    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()))
765
      ("no admin or different admins: %8X, %8X\n",
766
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
767
768
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
769
770
    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
771
772
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
773
      *xIterator *= *yIterator; 
774
775
776
777
778
779
780
781

    return x;
  }

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

784
785
786
    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()))
787
      ("no admin or different admins: %8X, %8X\n",
788
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
789
790
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
791
    T dot = 0;
792
793
794

    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
795
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
796
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
797
      dot += (*xIterator) * (*yIterator);
798

Thomas Witkowski's avatar
Thomas Witkowski committed
799
    return dot;
800
801
802
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
803
804
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
805
  {
806
807
808
    // Unfortunately needed
    using namespace mtl;

809
810
    FUNCNAME("DOFVector<T>::mv()");

811
812
813
814
    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());
815

816
817
    const DOFAdmin *rowAdmin = a.getRowFeSpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFeSpace()->getAdmin();
818
819

    TEST_EXIT((rowAdmin && colAdmin) && 
820
821
822
823
	      (((transpose == NoTranspose) && (colAdmin == x.getFeSpace()->getAdmin()) && 
		(rowAdmin == result.getFeSpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFeSpace()->getAdmin()) && 
		(colAdmin == result.getFeSpace()->getAdmin()))))
824
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
825
826
       a.getRowFeSpace()->getAdmin(), a.getColFeSpace()->getAdmin(), 
       x.getFeSpace()->getAdmin(), result.getFeSpace()->getAdmin());
827
828


829
830
831
832
833
    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 
834
835
836
837
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
838
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
839
840
841
  {
    FUNCNAME("DOFVector<T>::axpy()");

842
843
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
844

845
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
846

847
    TEST_EXIT((admin) && (admin == y.getFeSpace()->getAdmin()))
848
      ("no admin or different admins: %8X, %8X\n",
849
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
850
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
851
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
852
853
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
854
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
       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
Thomas Witkowski's avatar
Thomas Witkowski committed
874
      for (i = 0; i < maxI; i++)
875
	if (!admin->isDofFree(i))
876
	  y[i] += alpha * x[i];
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
  }

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

908
909
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
910

911
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
912

913
    TEST_EXIT(admin && (admin == y.getFeSpace()->getAdmin()))
914
      ("no admin or different admins: %8X, %8X\n",
915
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
916
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
917
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
918
919
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
920
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
       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
Thomas Witkowski's avatar
Thomas Witkowski committed
940
      for (i = 0; i < maxI; i++)
941
	if (!admin->isDofFree(i))
942
	  y[i] = alpha * y[i] + x[i];
943
944
945
946
947
948
949
950
951
  }

  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);
Thomas Witkowski's avatar
Thomas Witkowski committed
952
953
954
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = scal * (*vIterator);
955
956
957
958
959
960
961
962
963
964
965

    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);
Thomas Witkowski's avatar
Thomas Witkowski committed
966
967
968
969
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = (*vIterator) + scal;

970
971
972
    return result;
  }

973

974
  template<typename T>
975
  inline const DOFVector<T>& add(const DOFVector<T>& v1, 
976
977
978
979
980
981
				 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);