DOFVector.hh 57.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
#include <list>
Thomas Witkowski's avatar
Thomas Witkowski committed
23
#include <algorithm>
24
#include <math.h>
25

26
27
28
29
30
#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>

31
32
33
34
35
36
37
38
39
40
41
42
#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"
43
#include "Operator.h"
44
#include "Initfile.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
45
#include "Traverse.h"
46
47
48
#include "DualTraverse.h"

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
49
#include "parallel/MpiHelper.h"
50
#include "parallel/MeshDistributor.h"
51
#endif
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
84
85
86
87
88
89
90
91
92
93
// 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



94
95
96
namespace AMDiS {

  template<typename T>
97
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
98
99
    : feSpace(f),
      name(n),
100
      elementVector(f->getBasisFcts()->getNumber()),
101
      boundaryManager(nullptr)
102
  {    
Thomas Witkowski's avatar
Thomas Witkowski committed
103
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
104
    dim = feSpace->getMesh()->getDim();
105
106
107
108
109
  }
  

  template<typename T>
  DOFVectorBase<T>::~DOFVectorBase()
110
  {}
111

112
 
Thomas Witkowski's avatar
Thomas Witkowski committed
113
  template<typename T>
114
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
115
    : DOFVectorBase<T>(f, n),
116
117
      coarsenOperation(COARSE_INTERPOL),
      refineOperation(REFINE_INTERPOL)
118
  {
119
    vec.resize(0);
120
121
122
    init(f, n);
  } 

123

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

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
137
138
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
139
    
Thomas Witkowski's avatar
Thomas Witkowski committed
140
    if (this->boundaryManager)
141
      delete this->boundaryManager;
142
    
143
    vec.clear();
144
145
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
146
147
148
149
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
150
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
					  bool add)
  {
153
    std::vector<DegreeOfFreedom> indices(nBasFcts);
154
155
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), 
					     feSpace->getAdmin(),
156
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
157

158
    for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
159
      BoundaryCondition *condition = 
160
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : nullptr;
Thomas Witkowski's avatar
Thomas Witkowski committed
161

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

165
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
166
	  (*this)[irow] += factor * elVec[i];
167
	else  
168
	  (*this)[irow] = factor * elVec[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
169
170
171
172
      }
    }
  }

173
174
175
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
176
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
177

178
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
181
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
182
183
      nrm += (*vecIterator) * (*vecIterator);

184
185
186
187
188
189
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return sqrt(nrm);
190
191
  }

192

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

204
205
206
207
208
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

209
210
211
    return nrm;
  }

212

213
214
215
  template<typename T>
  T DOFVector<T>::asum() const
  {
216
    checkFeSpace(this->feSpace, vec);
217

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

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

    return nrm;
230
231
  }

232

233
234
235
  template<typename T>
  T DOFVector<T>::sum() const
  {
236
    checkFeSpace(this->feSpace, vec);
237

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

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

    return nrm;
250
251
  }

252

253
254
255
  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
256
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
257

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


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
267
    FUNCNAME_DBG("DOFVector<T>::copy()");
Thomas Witkowski's avatar
Thomas Witkowski committed
268

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

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


  template<typename T>
  T DOFVector<T>::min() const
287
  {    
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
  template<typename T>
  T DOFVector<T>::max() const 
306
  {    
307
    checkFeSpace(this->feSpace, vec);
308

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

314
315
316
317
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMax = m;
    MPI::COMM_WORLD.Allreduce(&localMax, &m, 1, MPI_DOUBLE, MPI_MAX);
#endif
318

319
    return m;
320
321
  }

322

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

329
330
331

  template<typename T>
  T DOFVector<T>::average() const
332
  {    
333
334
335
336
337
338
339
340
341
342
    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++;
    }

343
344
345
346
347
348
349
#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

350
351
352
353
    return m / count;
  }


354
355
356
  template<typename T>
  void DOFVector<T>::print() const
  {
357
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
358

359
    const DOFAdmin *admin = nullptr;
360
    const char *format;
361

362
363
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
364

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

405

406
  template<typename T>
407
  size_t DOFVector<T>::calcMemoryUsage() const
408
  {
409
    size_t result = 0;
410
411
412
413
414
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
415

416

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

425
    for (int i = 0; i < nBasisFcts; i++)
426
427
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

428
429
430
431
432
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localVal = val;
    MPI::COMM_WORLD.Allreduce(&localVal, &val, 1, MPI_DOUBLE, MPI_SUM);
#endif

433
434
435
    return val;
  }

436

437
438
439
  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
440
    FUNCNAME("DOFVector::interpol()");
441

Thomas Witkowski's avatar
Thomas Witkowski committed
442
    TEST_EXIT_DBG(fct)("No function to interpolate!\n");
443

444
    if (!this->getFeSpace()) {
445
446
447
448
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
449

450
    if (!(this->getFeSpace()->getAdmin())) {
451
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
452
	  this->getFeSpace()->getName().c_str());
453
454
      return;
    }
455
  
456
    if (!(this->getFeSpace()->getBasisFcts())) {
457
458
459
460
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
461

462
    if (!(fct)) {
463
      MSG("function that should be interpolated only pointer to nullptr, ");
464
465
466
      Msg::print("skipping interpolation\n");
      return;
    }
467

468
469
470
471
    const BasisFunction *basFct = this->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = this->getFeSpace()->getAdmin();
    int nBasFcts = basFct->getNumber();
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
472
    mtl::dense_vector<T> fctInterpolValues(nBasFcts);
473

Thomas Witkowski's avatar
Thomas Witkowski committed
474
    TraverseStack stack;
475
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
478
      basFct->interpol(elInfo, 0, nullptr, fct, fctInterpolValues);
479
480
481
482
483
      basFct->getLocalIndices(const_cast<Element*>(elInfo->getElement()), 
			      admin, myLocalIndices);
      for (int i = 0; i < nBasFcts; i++)
	vec[myLocalIndices[i]] = fctInterpolValues[i];

Thomas Witkowski's avatar
Thomas Witkowski committed
484
485
      elInfo = stack.traverseNext(elInfo);
    }
486
487
  }

488

489
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
490
  double DOFVector<T>::Int(int meshLevel, Quadrature* q) const
491
  {
492
    Mesh* mesh = this->feSpace->getMesh();
493

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

499
    FastQuadrature *quadFast = 
500
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
501

Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
504
505
506
507
    Flag flag = Mesh::FILL_COORDS | Mesh::FILL_DET;
    if (meshLevel == -1)
      flag |= Mesh::CALL_LEAF_EL;
    else
      flag |= Mesh::CALL_EL_LEVEL;

508
    T result; nullify(result);
509
    int nPoints = quadFast->getNumPoints();
510
    mtl::dense_vector<T> uh_vec(nPoints);
511
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
512
    ElInfo *elInfo =  stack.traverseFirst(mesh, meshLevel, flag);
513
514
    while (elInfo) {
      double det = elInfo->getDet();
515
      T normT; nullify(normT);
516
      this->getVecAtQPs(elInfo, nullptr, quadFast, uh_vec);
517
518
519
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;
520

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

524
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
525
    Parallel::mpi::globalAdd(result);
526
527
528
#endif
    
    return result;
529
530
  }

531
532
533
534
  
  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
535
  {
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
    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);
556
 	tmp += (*fct)(coords) * fastQuad->getWeight(iq);
557
558
559
560
561
562
563
564
      }

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

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
565
    Parallel::mpi::globalAdd(value);
566
567
568
569
570
571
572
573
574
575
576
#endif

    return value;
  }

  
  template<typename TOut, typename T>
  TOut integrate_Vec(const DOFVector<T> &vec,
		  AbstractFunction<TOut, T> *fct)
  {
    FUNCNAME("integrate_Vec()");
577
578
579
580
581
582
583
584
585
586

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

587
588
589
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);
590
591
592
593
594
595

    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);
596
      TOut tmp; nullify(tmp);
597
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
598
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq]);
599
600
601
602
603
604
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

605
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
606
    Parallel::mpi::globalAdd(value);
607
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
608
609
610

    return value;
  }
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
  
  
  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
651

652
653
654
655
656
657
658
659
660
      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
661
    Parallel::mpi::globalAdd(value);
662
663
664
665
666
667
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
#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) {
696
697
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, nullptr, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, nullptr, qp2);
698
699
700
701
702
703
704
705
706
707

      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
708
    Parallel::mpi::globalAdd(value);
709
710
711
712
713
714
#endif

    return value;
  }
  
  
Praetorius, Simon's avatar
Praetorius, Simon committed
715
  template<typename T1, typename T2>
716
717
718
  typename traits::mult_type<T1,T2>::type 
  integrate_VecTimesCoords(const DOFVector<T1> &vec,
			   AbstractFunction<T2, WorldVector<double> > *fct)
719
720
721
722
723
  {
    FUNCNAME("integrate_VecTimesCoords()");
    
    TEST_EXIT(fct)("No function defined!\n");
    
724
    typedef typename traits::mult_type<T1,T2>::type TOut;
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755

    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
756
    Parallel::mpi::globalAdd(value);
757
758
759
760
761
762
763
764
#endif

    return value;
  }

  
  template<typename TOut, typename T>
  TOut integrate_VecAndCoords(const DOFVector<T> &vec,
765
			   BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct)
Praetorius, Simon's avatar
Praetorius, Simon committed
766
767
768
769
770
771
772
773
774
775
776
777
  {
    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);

778
779
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());
    WorldVector<double> coordsAtQP;
Praetorius, Simon's avatar
Praetorius, Simon committed
780

781
782
    TOut value; nullify(value);
    
Praetorius, Simon's avatar
Praetorius, Simon committed
783
784
785
786
787
    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);
788
      TOut tmp; nullify(tmp);
Praetorius, Simon's avatar
Praetorius, Simon committed
789
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
790
791
	elInfo->coordToWorld(fastQuad->getLambda(iq), coordsAtQP);
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq], coordsAtQP);
Praetorius, Simon's avatar
Praetorius, Simon committed
792
793
794
795
796
797
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

798
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
799
    Parallel::mpi::globalAdd(value);
800
#endif
801

802
803
    return value;
  }
804

805
  
806
807
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
808
  {  
809
810
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
811
    if (!q) {
812
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
813
814
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
815

816
817
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
818

819
820
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
821
    mtl::dense_vector<T> uh_vec(nPoints);
822
823
824
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
825
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
826
827
828
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
829
      this->getVecAtQPs(elInfo, nullptr, quadFast, uh_vec);
830
831
832
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
833

834
835
      elInfo = stack.traverseNext(elInfo);
    }
836

837
838
839
840
841
842
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
843
844
  }

845

846
847
848
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
849
850
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
851
    if (!q) {
852
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
853
854
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
855

856
857
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
858

859
860
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
861
    mtl::dense_vector<T> uh_vec(nPoints);
862
863
864
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
865
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
866
867
868
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
869
      this->getVecAtQPs(elInfo, nullptr, quadFast, uh_vec);
870
871
872
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;
873

874
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
875
    }
876

877
878
879
880
881
882
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
883
884
  }

885

886
  template<typename T>  
887
888
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
889
890
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
891
    if (!q) {
892
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
893
894
895
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

896
897
898
899
900
901
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
902
    mtl::dense_vector<WorldVector<T> > grduh_vec(nPoints);
903
904
905
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
906
907
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
908
909
910
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
911
      this->getGrdAtQPs(elInfo, nullptr, quadFast, grduh_vec);
912
913
914
915
916
917
918

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

920
      result += det * normT;
921

922
923
      elInfo = stack.traverseNext(elInfo);
    }
924

925
926
927
928
929
930
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
931
932
  }

933

934
935

  template<typename T>
936
  bool DOFVector<T>::getDofIdxAtPoint(WorldVector<double> &p, 
Thomas Witkowski's avatar
Thomas Witkowski committed
937
938
939
940
					    DegreeOfFreedom &idx, 
					    ElInfo *oldElInfo, 
					    bool useOldElInfo) const
  { 
Praetorius, Simon's avatar
Praetorius, Simon committed
941
942
    Mesh *mesh = this->feSpace->getMesh();
    const BasisFunction *basFcts = this->feSpace->getBasisFcts();
943
944
945
946
947
948
949
950
951

    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
952
    bool inside = false;
953
954

    if (oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
955
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), nullptr, nullptr);
956
957
      delete oldElInfo;
    } else {
958
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
959
960
961
962
963
    }
    
    if (oldElInfo)
      oldElInfo = elInfo;
    
Nestler, Michael's avatar
Nestler, Michael committed
964
965
    if (!inside){
      delete elInfo; 
966
      return false;
Nestler, Michael's avatar
Nestler, Michael committed
967
    }
968
    
Praetorius, Simon's avatar
Praetorius, Simon committed
969
    basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices);
970
971
    
    WorldVector<double> coord;