DOFVector.hh 57.3 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
#endif
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
84
85
86
87
88
89
90
91
92
// 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



93
94
95
namespace AMDiS {

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

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

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

122

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

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

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

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

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

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

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

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

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

    return sqrt(nrm);
189
190
  }

191

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

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

208
209
210
    return nrm;
  }

211

212
213
214
  template<typename T>
  T DOFVector<T>::asum() const
  {
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
    checkFeSpace(this->feSpace, vec);
236

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

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

    return nrm;
249
250
  }

251

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

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


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

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

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


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

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

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

299
300
301
    return m;
  }

302

303
304
  template<typename T>
  T DOFVector<T>::max() const 
305
  {    
306
    checkFeSpace(this->feSpace, vec);
307

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

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

318
    return m;
319
320
  }

321

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

328
329
330

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

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

349
350
351
352
    return m / count;
  }


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

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

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

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

404

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

    return result;
  }
414

415

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

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

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

432
433
434
    return val;
  }

435

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
473
    TraverseStack stack;
474
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
475
476
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
477
      basFct->interpol(elInfo, 0, nullptr, fct, fctInterpolValues);
478
479
480
481
482
      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
483
484
      elInfo = stack.traverseNext(elInfo);
    }
485
486
  }

487

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

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

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

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

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

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

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

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

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

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

    return value;
  }

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

    return value;
  }
610
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
  
  
  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
650

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

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

    return value;
  }
  
  
Praetorius, Simon's avatar
Praetorius, Simon committed
714
  template<typename T1, typename T2>
715
716
717
718
719
720
721
722
723
724
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
  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
754
    Parallel::mpi::globalAdd(value);
755
756
757
758
759
760
761
762
#endif

    return value;
  }

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

796
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
797
    Parallel::mpi::globalAdd(value);
798
#endif
799

800
801
    return value;
  }
802

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

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

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

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

832
833
      elInfo = stack.traverseNext(elInfo);
    }
834

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

843

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

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

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

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

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

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

883

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

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

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

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

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

918
      result += det * normT;
919

920
921
      elInfo = stack.traverseNext(elInfo);
    }
922

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

931

932
933

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

    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
950
    bool inside = false;
951
952

    if (oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
953
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), nullptr, nullptr);
954
955
      delete oldElInfo;
    } else {
956
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
957
958
959
960
961
    }
    
    if (oldElInfo)
      oldElInfo = elInfo;
    
Nestler, Michael's avatar
Nestler, Michael committed
962
963
    if (!inside){
      delete elInfo; 
964
      return false;
Nestler, Michael's avatar
Nestler, Michael committed
965
    }
966
    
Praetorius, Simon's avatar
Praetorius, Simon committed
967
    basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices);
968
969
    
    WorldVector<double> coord;
970
971
972
    int minIdx = -1;
    double minDist = 1.e15;
    
973
974
975
    for (int i = 0; i < numBasFcts; i++) {
      elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
      WorldVector<double> dist = coord - p;
976
977
978
979
      double newDist = norm(dist);