DOFVector.hh 56.6 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
#include "DualTraverse.h"

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
40
#include "parallel/MpiHelper.h"
41
#endif
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
80
81
82
83
// 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



84
85
86
namespace AMDiS {

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

  template<typename T>
  DOFVectorBase<T>::~DOFVectorBase()
100
  {}
101

102
 
Thomas Witkowski's avatar
Thomas Witkowski committed
103
  template<typename T>
104
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
105
    : DOFVectorBase<T>(f, n),
106
      coarsenOperation(COARSE_INTERPOL)
107
  {
108
    vec.resize(0);
109
110
111
    init(f, n);
  } 

112

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

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
126
127
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
128
    
Thomas Witkowski's avatar
Thomas Witkowski committed
129
    if (this->boundaryManager)
130
      delete this->boundaryManager;
131
    
132
    vec.clear();
133
134
135
136
137
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
138
139
140
141
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
142
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
143
144
145
146
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

147
148
149
    std::vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
150

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

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

158
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
159
	  (*this)[irow] += factor * elVec[i];
160
161
	else  
	  (*this)[irow] = factor * elVec[i];	
Thomas Witkowski's avatar
Thomas Witkowski committed
162
163
164
165
      }
    }
  }

166
167
168
169
170
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

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

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

179
180
181
182
183
184
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return sqrt(nrm);
185
186
  }

187

188
189
190
191
192
  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

193
    checkFeSpace(this->feSpace, vec);
194
195
    
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
196
197
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
198
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
199
200
      nrm += (*vecIterator) * (*vecIterator);

201
202
203
204
205
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

206
207
208
    return nrm;
  }

209

210
211
212
  template<typename T>
  T DOFVector<T>::asum() const
  {
213
    FUNCNAME("DOFVector<T>::asum()");
214

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

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

223
224
225
226
227
228
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
229
230
  }

231

232
233
234
  template<typename T>
  T DOFVector<T>::sum() const
  {
235
    FUNCNAME("DOFVector<T>::sum()");
236

237
    checkFeSpace(this->feSpace, vec);
238

239
    double nrm = 0.0;    
Thomas Witkowski's avatar
Thomas Witkowski committed
240
241
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
242
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
243
244
      nrm += *vecIterator;

245
246
247
248
249
250
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
251
252
  }

253

254
255
256
  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
303
304
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMin = m;
    MPI::COMM_WORLD.Allreduce(&localMin, &m, 1, MPI_DOUBLE, MPI_MIN);
#endif

305
306
307
    return m;
  }

308

309
310
311
  template<typename T>
  T DOFVector<T>::max() const 
  {
312
313
    FUNCNAME("DOFVector<T>::max()");
    
314
    checkFeSpace(this->feSpace, vec);
315

316
    T m;    
317
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
318
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
319
      m = std::max(m, *vecIterator);
320

321
322
323
324
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMax = m;
    MPI::COMM_WORLD.Allreduce(&localMax, &m, 1, MPI_DOUBLE, MPI_MAX);
#endif
325

326
    return m;
327
328
  }

329

Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
332
333
334
335
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351

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

352
353
354
355
356
357
358
#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

359
360
361
362
    return m / count;
  }


363
364
365
  template<typename T>
  void DOFVector<T>::print() const
  {
366
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
367

368
    const DOFAdmin *admin = NULL;
369
    const char *format;
370

371
372
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
373

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

414

415
416
417
418
419
420
421
422
423
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
424

425

426
427
428
429
  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda, 
			     DegreeOfFreedom* dof_indices)
  {
430
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFeSpace()->getBasisFcts());
431
    int nBasisFcts = phi->getNumber();
432
433
    T val = 0.0;

434
    for (int i = 0; i < nBasisFcts; i++)
435
436
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

437
438
439
440
441
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localVal = val;
    MPI::COMM_WORLD.Allreduce(&localVal, &val, 1, MPI_DOUBLE, MPI_SUM);
#endif

442
443
444
    return val;
  }

445

446
447
448
  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
449
    FUNCNAME("DOFVector::interpol()");
450

Thomas Witkowski's avatar
Thomas Witkowski committed
451
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
452
453

    interFct = fct;
454

455
    if (!this->getFeSpace()) {
456
457
458
459
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
460

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

473
474
475
476
477
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
478
479
480

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
481
    TraverseStack stack;
482
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
483
484
485
486
487
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
488
489
  }

490

491
  template<typename T>
492
  void DOFVector<T>::interpolFct(ElInfo* elinfo)
493
  {
494
495
    const BasisFunction *basFct = traverseVector->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFeSpace()->getAdmin();
496
497
498
    const T *inter_val = 
      const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
						   traverseVector->interFct, NULL);
499

500
    int nBasFcts = basFct->getNumber();
501
502
503
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
    basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), 
			    admin, myLocalIndices);
504
    for (int i = 0; i < nBasFcts; i++)
505
      (*traverseVector)[myLocalIndices[i]] = inter_val[i];
506
507
  }

508

509
510
511
  template<typename T>
  double DOFVector<T>::Int(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
512
    FUNCNAME("DOFVector::Int()");
513
  
514
    Mesh* mesh = this->feSpace->getMesh();
515

Thomas Witkowski's avatar
Thomas Witkowski committed
516
    if (!q) {
517
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
518
519
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
520

521
    FastQuadrature *quadFast = 
522
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
523

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

538
539
      elInfo = stack.traverseNext(elInfo);
    }
540

541
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
542
    mpi::globalAdd(result);
543
544
545
#endif
    
    return result;
546
547
  }

548
549
550
551
  
  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
552
  {
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
    FUNCNAME("integrate_Coords()");
    
    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(), feSpace->getBasisFcts()->getDegree());
    Quadrature* quad = 
      Quadrature::provideQuadrature(feSpace->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);

    WorldVector<double> coords;

    TOut value; nullify(value);
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
    ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, traverseFlag);
    while (elInfo) {
    TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
 	tmp += fastQuad->getWeight(iq) * (*fct)(coords);
      }

      value += tmp * elInfo->getDet();
      
      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif

    return value;
  }

  
  template<typename TOut, typename T>
  TOut integrate_Vec(const DOFVector<T> &vec,
		  AbstractFunction<TOut, T> *fct)
  {
    FUNCNAME("integrate_Vec()");
594
595
596
597
598
599
600
601
602
603

    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(),
		       2 * vec.getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad =
      Quadrature::provideQuadrature(vec.getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vec.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

604
605
606
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);
607
608
609
610
611
612

    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(vec.getFeSpace()->getMesh(), -1, traverseFlag);
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);
613
      TOut tmp; nullify(tmp);
614
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
615
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq]);
616
617
618
619
620
621
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

622
623
624
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
625
626
627

    return value;
  }
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
  
  
  template<typename TOut, typename T1, typename T2>
  TOut integrate_Vec2(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct)
  {
    if (vec1.getFeSpace()->getMesh() == vec2.getFeSpace()->getMesh())
      return int_Vec2_SingleMesh(vec1, vec2, fct);
    else
      return int_Vec2_DualMesh(vec1, vec2, fct);
  }


  template<typename TOut, typename T1, typename T2>
  TOut int_Vec2_SingleMesh(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct)
  {
    FUNCNAME("intDualmesh()");
    
    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad = 
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<T1> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<T2> qp2(fastQuad->getNumPoints());

    TOut value; nullify(value);
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
    while (elInfo) {
      vec1.getVecAtQPs(elInfo, quad, fastQuad, qp1);
      vec2.getVecAtQPs(elInfo, quad, fastQuad, qp2);
Praetorius, Simon's avatar
Praetorius, Simon committed
668

669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);     
      value += tmp * elInfo->getDet();
      
      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif
    
    return value;
  }


  template<typename TOut, typename T1, typename T2>
  TOut int_Vec2_DualMesh(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct)
  {
    FUNCNAME("intDualmesh()");
    
    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad = 
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<T1> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<T2> qp2(fastQuad->getNumPoints());

    TOut value; nullify(value);
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    DualTraverse dualTraverse;    
    DualElInfo dualElInfo;
    bool cont = dualTraverse.traverseFirst(vec1.getFeSpace()->getMesh(),
					   vec2.getFeSpace()->getMesh(),
					   -1, -1, traverseFlag, traverseFlag,
					   dualElInfo);
    while (cont) {
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);

      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < quad->getNumPoints(); iq++)
 	tmp += quad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);      
      value += tmp * dualElInfo.smallElInfo->getDet();
      
      cont = dualTraverse.traverseNext(dualElInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif

    return value;
  }
  
  
Praetorius, Simon's avatar
Praetorius, Simon committed
732
  template<typename T1, typename T2>
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
  typename ProductType<T1,T2>::type integrate_VecTimesCoords(const DOFVector<T1> &vec,
		   AbstractFunction<T2, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate_VecTimesCoords()");
    
    TEST_EXIT(fct)("No function defined!\n");
    
    typedef typename ProductType<T1,T2>::type TOut;

    const FiniteElemSpace *feSpace = vec.getFeSpace();
    Mesh *mesh = feSpace->getMesh();

    int deg = std::max(fct->getDegree(), 2 * feSpace->getBasisFcts()->getDegree());
    Quadrature* quad = Quadrature::provideQuadrature(mesh->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<T1> qp(fastQuad->getNumPoints());
    WorldVector<double> coords;

    TOut value; nullify(value);
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);

      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
 	tmp += fastQuad->getWeight(iq) * (qp[iq] * (*fct)(coords));
      }

      value += tmp * elInfo->getDet();
      
      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif

    return value;
  }

  
  template<typename TOut, typename T>
  TOut integrate_VecAndCoords(const DOFVector<T> &vec,
			   BinaryAbstractFunction<T, T, WorldVector<double> > *fct)
Praetorius, Simon's avatar
Praetorius, Simon committed
782
783
784
785
786
787
788
789
790
791
792
793
  {
    FUNCNAME("integrate_VecAndCoords()");

    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(),
		       2 * vec.getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad =
      Quadrature::provideQuadrature(vec.getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vec.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

794
795
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());
    WorldVector<double> coordsAtQP;
Praetorius, Simon's avatar
Praetorius, Simon committed
796

797
798
    TOut value; nullify(value);
    
Praetorius, Simon's avatar
Praetorius, Simon committed
799
800
801
802
803
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(vec.getFeSpace()->getMesh(), -1, traverseFlag);
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);
804
      TOut tmp; nullify(tmp);
Praetorius, Simon's avatar
Praetorius, Simon committed
805
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
806
807
	elInfo->coordToWorld(fastQuad->getLambda(iq), coordsAtQP);
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq], coordsAtQP);
Praetorius, Simon's avatar
Praetorius, Simon committed
808
809
810
811
812
813
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

814
815
816
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif
817

818
819
    return value;
  }
820

821
  
822
823
824
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
825
    FUNCNAME("DOFVector::L1Norm()");
826
  
827
828
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
829
    if (!q) {
830
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
831
832
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
833

834
835
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
836

837
838
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
839
    mtl::dense_vector<T> uh_vec(nPoints);
840
841
842
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
843
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
844
845
846
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
847
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
848
849
850
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
851

852
853
      elInfo = stack.traverseNext(elInfo);
    }
854

855
856
857
858
859
860
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
861
862
  }

863

864
865
866
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
867
    FUNCNAME("DOFVector::L2NormSquare()");
868

869
870
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
871
    if (!q) {
872
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
873
874
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
875

876
877
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
878

879
880
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
881
    mtl::dense_vector<T> uh_vec(nPoints);
882
883
884
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
885
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
886
887
888
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
889
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
890
891
892
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;
893

894
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
895
    }
896

897
898
899
900
901
902
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
903
904
  }

905

906
  template<typename T>  
907
908
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
909
    FUNCNAME("DOFVector::H1NormSquare()");
910

911
912
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
913
    if (!q) {
914
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
915
916
917
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

918
919
920
921
922
923
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
924
    mtl::dense_vector<WorldVector<T> > grduh_vec(nPoints);
925
926
927
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
928
929
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
930
931
932
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
933
      this->getGrdAtQPs(elInfo, NULL, quadFast, grduh_vec);
934
935
936
937
938
939
940

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

942
      result += det * normT;
943

944
945
      elInfo = stack.traverseNext(elInfo);
    }
946

947
948
949
950
951
952
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
953
954
  }

955

956
957

  template<typename T>
Praetorius, Simon's avatar
Praetorius, Simon committed
958
  const bool DOFVector<T>::getDOFidxAtPoint(WorldVector<double> &p, DegreeOfFreedom &idx, ElInfo *oldElInfo, bool useOldElInfo) const
959
960
  { FUNCNAME("DOFVector::getDOFidxAtPoint()");

Praetorius, Simon's avatar
Praetorius, Simon committed
961
962
    Mesh *mesh = this->feSpace->getMesh();
    const BasisFunction *basFcts = this->feSpace->getBasisFcts();
963
964
965
966
967
968
969
970
971

    int dim = mesh->getDim();
    int numBasFcts = basFcts->getNumber();

    std::vector<DegreeOfFreedom> localIndices(numBasFcts);
    DimVec<double> lambda(dim, NO_INIT);

    ElInfo *elInfo = mesh->createNewElInfo();
    idx = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
972
    bool inside = false;
973
974
975
976
977
978
979
980
981
982
983
984
985
986

    if (oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
      delete oldElInfo;
    } else {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
    }
    
    if (oldElInfo)
      oldElInfo = elInfo;
    
    if (!inside)
      return false;
    
Praetorius, Simon's avatar
Praetorius, Simon committed
987
    basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices);
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
    FixVec<WorldVector<double>, VERTEX> coords = elInfo->getCoords();
    int minIdx = -1;
    double minDist = 1.e15;
    
    for (int i = 0; i < coords.size(); i++) {
      WorldVector<double> dist = coords[i] - p;
      double newDist = norm(dist);
      if (newDist < minDist) {
	minIdx = i;
	minDist = newDist;
      }
    }
    
    if (minIdx >= 0)
      idx = localIndices[minIdx];
    
    if(!oldElInfo) delete elInfo;
    return inside;
  };


1009
1010
  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
1011
					std::vector<DegreeOfFreedom> &newDOF)
1012
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1013
1014
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
1015
	vec[newDOF[i]] = vec[i];
1016
1017
  }

1018

1019
1020
1021
1022
  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1023
1024
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
1025
	 op != this->operators.end(); ++op)
1026
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
1027

1028
1029
1030
    return fillFlag;
  }