DOFVector.hh 31.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
      nrm += (*vecIterator) * (*vecIterator);

174
175
176
177
178
179
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return sqrt(nrm);
180
181
  }

182

183
184
185
186
187
  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

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

196
197
198
199
200
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

201
202
203
    return nrm;
  }

204

205
206
207
  template<typename T>
  T DOFVector<T>::asum() const
  {
208
    FUNCNAME("DOFVector<T>::asum()");
209

210
    checkFeSpace(this->feSpace, vec);
211

212
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
215
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
216
217
      nrm += abs(*vecIterator);

218
219
220
221
222
223
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
224
225
  }

226

227
228
229
  template<typename T>
  T DOFVector<T>::sum() const
  {
230
    FUNCNAME("DOFVector<T>::sum()");
231

232
    checkFeSpace(this->feSpace, vec);
233

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

240
241
242
243
244
245
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
246
247
  }

248

249
250
251
  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
252
    FUNCNAME("DOFVector<T>::set()");
253

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

256
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
257
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
258
259
260
261
262
263
264
      *vecIterator = alpha ; 
  }


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

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

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


  template<typename T>
  T DOFVector<T>::min() const
  {
286
287
    FUNCNAME("DOFVector<T>::min()");
    
288
    checkFeSpace(this->feSpace, vec);
289
290

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

295
296
297
298
299
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMin = m;
    MPI::COMM_WORLD.Allreduce(&localMin, &m, 1, MPI_DOUBLE, MPI_MIN);
#endif

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
317
318
319
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMax = m;
    MPI::COMM_WORLD.Allreduce(&localMax, &m, 1, MPI_DOUBLE, MPI_MAX);
#endif
320

321
    return m;
322
323
  }

324

Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
327
328
329
330
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346

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

347
348
349
350
351
352
353
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localSum = m;
    int localCount = count;
    MPI::COMM_WORLD.Allreduce(&localSum, &m, 1, MPI_DOUBLE, MPI_SUM);
    MPI::COMM_WORLD.Allreduce(&localCount, &count, 1, MPI_INT, MPI_SUM);
#endif

354
355
356
357
    return m / count;
  }


358
359
360
  template<typename T>
  void DOFVector<T>::print() const
  {
361
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
362

363
    const DOFAdmin *admin = NULL;
364
    const char *format;
365

366
367
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
368

369
    MSG("Vec `%s':\n", this->name.c_str());
370
    int j = 0;
371
372
373
374
375
376
377
378
379
    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
380
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
381
	if ((j % 3) == 0) {
382
383
	  if (j) 
	    Msg::print("\n");
384
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
385
	} else {
386
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
387
	}
388
	j++;
389
      }
390
      Msg::print("\n");
391
    } else {
392
393
	MSG("no DOFAdmin, print whole vector.\n");
    
394
395
396
397
398
399
	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 {
400
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
401
	  }
402
403
404
405
406
407
408
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }

409

410
411
412
413
414
415
416
417
418
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
419

420

421
422
423
424
  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda, 
			     DegreeOfFreedom* dof_indices)
  {
425
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFeSpace()->getBasisFcts());
426
    int nBasisFcts = phi->getNumber();
427
428
    T val = 0.0;

429
    for (int i = 0; i < nBasisFcts; i++)
430
431
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

432
433
434
435
436
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localVal = val;
    MPI::COMM_WORLD.Allreduce(&localVal, &val, 1, MPI_DOUBLE, MPI_SUM);
#endif

437
438
439
    return val;
  }

440

441
442
443
  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
444
    FUNCNAME("DOFVector::interpol()");
445

Thomas Witkowski's avatar
Thomas Witkowski committed
446
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
447
448

    interFct = fct;
449

450
    if (!this->getFeSpace()) {
451
452
453
454
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
455

456
    if (!(this->getFeSpace()->getAdmin())) {
457
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
458
	  this->getFeSpace()->getName().c_str());
459
460
      return;
    }
461
  
462
    if (!(this->getFeSpace()->getBasisFcts())) {
463
464
465
466
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
467

468
469
470
471
472
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
473
474
475

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
476
    TraverseStack stack;
477
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
478
479
480
481
482
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
483
484
  }

485

486
  template<typename T>
487
  void DOFVector<T>::interpolFct(ElInfo* elinfo)
488
  {
489
490
    const BasisFunction *basFct = traverseVector->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFeSpace()->getAdmin();
491
492
493
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
494

495
    int nBasFcts = basFct->getNumber();
496
497
498
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), 
			    admin, myLocalIndices);
499
    for (int i = 0; i < nBasFcts; i++)
500
      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
501
502
  }

503

504
505
506
  template<typename T>
  double DOFVector<T>::Int(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
507
    FUNCNAME("DOFVector::Int()");
508
  
509
    Mesh* mesh = this->feSpace->getMesh();
510

Thomas Witkowski's avatar
Thomas Witkowski committed
511
    if (!q) {
512
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
513
514
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
515

516
    FastQuadrature *quadFast = 
517
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
518

519
520
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
521
    mtl::dense_vector<T> uh_vec(nPoints);
522
523
524
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
525
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
526
527
528
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
529
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
530
531
532
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;
533

534
535
      elInfo = stack.traverseNext(elInfo);
    }
536

537
538
539
540
541
542
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
543
544
  }

545

546
547
548
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
549
    FUNCNAME("DOFVector::L1Norm()");
550
  
551
552
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
553
    if (!q) {
554
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
555
556
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
557

558
559
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
560

561
562
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
563
    mtl::dense_vector<T> uh_vec(nPoints);
564
565
566
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
567
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
568
569
570
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
571
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
572
573
574
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
575

576
577
      elInfo = stack.traverseNext(elInfo);
    }
578

579
580
581
582
583
584
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
585
586
  }

587

588
589
590
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
591
    FUNCNAME("DOFVector::L2NormSquare()");
592

593
594
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
595
    if (!q) {
596
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
597
598
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
599

600
601
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
602

603
604
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
605
    mtl::dense_vector<T> uh_vec(nPoints);
606
607
608
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
609
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
610
611
612
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
613
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
614
615
616
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;
617

618
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
619
    }
620

621
622
623
624
625
626
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
627
628
  }

629

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

635
636
    Mesh* mesh = this->feSpace->getMesh();

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

642
643
644
645
646
647
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
648
    std::vector<WorldVector<T> > grduh_vec(nPoints);
649
650
651
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
652
653
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
654
655
656
657
658
659
660
661
662
663
664
    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;
      }
665

666
      result += det * normT;
667

668
669
      elInfo = stack.traverseNext(elInfo);
    }
670

671
672
673
674
675
676
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
677
678
  }

679

680
681
  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
682
					std::vector<DegreeOfFreedom> &newDOF)
683
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
684
685
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
686
	vec[newDOF[i]] = vec[i];
687
688
  }

689

690
691
692
693
  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
694
695
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
696
	 op != this->operators.end(); ++op)
697
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
698

699
700
701
    return fillFlag;
  }

702

Thomas Witkowski's avatar
Thomas Witkowski committed
703
704
705
706
707
  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
708
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
709
710
711
      (*it)->finishAssembling();
  }

712

713
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
714
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
715
  {
716
    this->feSpace = rhs.feSpace;
717
    this->dim = rhs.dim;
718
    this->nBasFcts = rhs.nBasFcts;
719
    vec = rhs.vec;
720
    this->elementVector.change_dim(this->nBasFcts);
721
722
723
724
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
725
    
726
    if (rhs.boundaryManager) {
727
728
729
      if (this->boundaryManager) 
	delete this->boundaryManager; 

730
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
731
    } else {
732
      this->boundaryManager = NULL;
733
734
    }

735
736
737
738
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    this->rankDofs = rhs.rankDofs;
#endif

739
740
741
    return *this;
  }

742

743
744
745
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal)
  {
746
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, T scal)");
747

748
749
    TEST_EXIT_DBG(x.getFeSpace() && x.getFeSpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFeSpace(), x.getFeSpace()->getAdmin());
750

Thomas Witkowski's avatar
Thomas Witkowski committed
751
752
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
753
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
754
      (*vecIterator) *= scal; 
755

756
757
758
759
760
761
762
    return x;
  }


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

779
780
781
    return x;
  }

782

783
784
785
  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y)
  {
786
    FUNCNAME("DOFVector<T>::operator-=(DOFVector<T>& x, const DOFVector<T>& y)");
787

788
789
790
791
    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()))
792
      ("no admin or different admins: %8X, %8X\n",
793
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
794
795
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
796
797
    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
798
799
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
800
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
801

802
803
804
    return x;
  }

805

806
807
808
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
809
810
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
811
812
813
814
    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()))
815
      ("no admin or different admins: %8X, %8X\n",
816
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
817
818
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
819
820
    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
821
822
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
823
      *xIterator *= *yIterator; 
824
825
826
827

    return x;
  }

828

829
830
831
832
  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::operator*(DOFVector<T>& x, DOFVector<T>& y)");
833
    const DOFAdmin *admin = NULL;
834

835
836
837
    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()))
838
      ("no admin or different admins: %8X, %8X\n",
839
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
840
841
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
842
    T dot = 0;
843
844
845

    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
846
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
847
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
848
      dot += (*xIterator) * (*yIterator);
849

Thomas Witkowski's avatar
Thomas Witkowski committed
850
    return dot;
851
852
  }

853

854
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
855
856
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
857
  {
858
859
860
    // Unfortunately needed
    using namespace mtl;

861
862
    FUNCNAME("DOFVector<T>::mv()");

863
864
865
866
    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());
867

868
869
    const DOFAdmin *rowAdmin = a.getRowFeSpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFeSpace()->getAdmin();
870
871

    TEST_EXIT((rowAdmin && colAdmin) && 
872
873
874
875
	      (((transpose == NoTranspose) && (colAdmin == x.getFeSpace()->getAdmin()) && 
		(rowAdmin == result.getFeSpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFeSpace()->getAdmin()) && 
		(colAdmin == result.getFeSpace()->getAdmin()))))
876
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
877
878
       a.getRowFeSpace()->getAdmin(), a.getColFeSpace()->getAdmin(), 
       x.getFeSpace()->getAdmin(), result.getFeSpace()->getAdmin());
879
880


881
882
883
884
885
    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 
886
887
888
      ERROR_EXIT("transpose = %d\n", transpose);
  }

889

890
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
891
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
892
893
894
  {
    FUNCNAME("DOFVector<T>::axpy()");

895
896
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
897

898
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
899

900
    TEST_EXIT((admin) && (admin == y.getFeSpace()->getAdmin()))
901
      ("no admin or different admins: %8X, %8X\n",
902
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
903
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
904
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
905
906
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
907
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
908
909
       admin->getUsedSize());

910
911
912
    for (int i = 0; i < y.getSize(); i++)
      if (!admin->isDofFree(i))
	y[i] += alpha * x[i];
913
914
  }

915

916
917
918
919
920
921
922
  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d)
  {
    static DOFVector<T> result;
    return mult(d, v, result);
  }

923

924
925
926
927
928
929
930
  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v)
  {
    static DOFVector<T> result;
    return mult(d, v, result);
  }

Thomas Witkowski's avatar