DOFVector.hh 58.1 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, bool addToSynch)
115
    : DOFVectorBase<T>(f, n),
116
117
      coarsenOperation(COARSE_INTERPOL),
      refineOperation(REFINE_INTERPOL)
118
  {
119
    vec.resize(0);
120
    init(f, n, addToSynch);
121
122
  } 

123

124
  template<typename T>
125
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n, bool addToSynch)
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
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if ( Parallel::MeshDistributor::globalMeshDistributor != NULL) 
      Parallel::MeshDistributor::globalMeshDistributor->addInterchangeVector(this);
#endif
136
137
138
139
140
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
141
142
143
144
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if ( Parallel::MeshDistributor::globalMeshDistributor != NULL) 
      Parallel::MeshDistributor::globalMeshDistributor->removeInterchangeVector(this);
#endif
145
146
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
147
    
Thomas Witkowski's avatar
Thomas Witkowski committed
148
    if (this->boundaryManager)
149
      delete this->boundaryManager;
150
    
151
    vec.clear();
152
153
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
154
155
156
157
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
158
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
					  bool add)
  {
161
    std::vector<DegreeOfFreedom> indices(nBasFcts);
162
163
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), 
					     feSpace->getAdmin(),
164
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
165

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

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

173
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
174
	  (*this)[irow] += factor * elVec[i];
175
	else  
176
	  (*this)[irow] = factor * elVec[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
177
178
179
180
      }
    }
  }

181
182
183
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
184
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
185

186
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
189
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
190
191
      nrm += (*vecIterator) * (*vecIterator);

192
193
194
195
196
197
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return sqrt(nrm);
198
199
  }

200

201
202
203
  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
204
    checkFeSpace(this->feSpace, vec);
205
206
    
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
207
208
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
209
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
210
211
      nrm += (*vecIterator) * (*vecIterator);

212
213
214
215
216
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

217
218
219
    return nrm;
  }

220

221
222
223
  template<typename T>
  T DOFVector<T>::asum() const
  {
224
    checkFeSpace(this->feSpace, vec);
225

226
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
229
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
230
231
      nrm += abs(*vecIterator);

232
233
234
235
236
237
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
238
239
  }

240

241
242
243
  template<typename T>
  T DOFVector<T>::sum() const
  {
244
    checkFeSpace(this->feSpace, vec);
245

246
    double nrm = 0.0;    
Thomas Witkowski's avatar
Thomas Witkowski committed
247
248
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
249
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
250
251
      nrm += *vecIterator;

252
253
254
255
256
257
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
258
259
  }

260

261
262
263
  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
264
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
265

266
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
267
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
268
269
270
271
272
273
274
      *vecIterator = alpha ; 
  }


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

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

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


  template<typename T>
  T DOFVector<T>::min() const
295
  {    
296
    checkFeSpace(this->feSpace, vec);
297
298

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

303
304
305
306
307
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMin = m;
    MPI::COMM_WORLD.Allreduce(&localMin, &m, 1, MPI_DOUBLE, MPI_MIN);
#endif

308
309
310
    return m;
  }

311

312
313
  template<typename T>
  T DOFVector<T>::max() const 
314
  {    
315
    checkFeSpace(this->feSpace, vec);
316

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

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

327
    return m;
328
329
  }

330

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

337
338
339

  template<typename T>
  T DOFVector<T>::average() const
340
  {    
341
342
343
344
345
346
347
348
349
350
    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++;
    }

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

358
359
360
361
    return m / count;
  }


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

367
    const DOFAdmin *admin = nullptr;
368
    const char *format;
369

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

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

413

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

    return result;
  }
423

424

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

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

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

441
442
443
    return val;
  }

444

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

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

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

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

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

476
477
478
479
    const BasisFunction *basFct = this->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = this->getFeSpace()->getAdmin();
    int nBasFcts = basFct->getNumber();
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
480
    mtl::dense_vector<T> fctInterpolValues(nBasFcts);
481

Thomas Witkowski's avatar
Thomas Witkowski committed
482
    TraverseStack stack;
483
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
484
485
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
486
      basFct->interpol(elInfo, 0, nullptr, fct, fctInterpolValues);
487
488
489
490
491
      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
492
493
      elInfo = stack.traverseNext(elInfo);
    }
494
495
  }

496

497
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
498
  double DOFVector<T>::Int(int meshLevel, Quadrature* q) const
499
  {
500
    Mesh* mesh = this->feSpace->getMesh();
501

Thomas Witkowski's avatar
Thomas Witkowski committed
502
    if (!q) {
503
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
506

507
    FastQuadrature *quadFast = 
508
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
509

Thomas Witkowski's avatar
Thomas Witkowski committed
510
511
512
513
514
515
    Flag flag = Mesh::FILL_COORDS | Mesh::FILL_DET;
    if (meshLevel == -1)
      flag |= Mesh::CALL_LEAF_EL;
    else
      flag |= Mesh::CALL_EL_LEVEL;

516
    T result; nullify(result);
517
    int nPoints = quadFast->getNumPoints();
518
    mtl::dense_vector<T> uh_vec(nPoints);
519
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
520
    ElInfo *elInfo =  stack.traverseFirst(mesh, meshLevel, flag);
521
522
    while (elInfo) {
      double det = elInfo->getDet();
523
      T normT; nullify(normT);
524
      this->getVecAtQPs(elInfo, nullptr, quadFast, uh_vec);
525
526
527
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;
528

529
530
      elInfo = stack.traverseNext(elInfo);
    }
531

532
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
533
    Parallel::mpi::globalAdd(result);
534
535
536
#endif
    
    return result;
537
538
  }

539
540
541
542
  
  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
543
  {
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
    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);
564
 	tmp += (*fct)(coords) * fastQuad->getWeight(iq);
565
566
567
568
569
570
571
572
      }

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

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
573
    Parallel::mpi::globalAdd(value);
574
575
576
577
578
579
580
581
582
583
584
#endif

    return value;
  }

  
  template<typename TOut, typename T>
  TOut integrate_Vec(const DOFVector<T> &vec,
		  AbstractFunction<TOut, T> *fct)
  {
    FUNCNAME("integrate_Vec()");
585
586
587
588
589
590
591
592
593
594

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

595
596
597
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);
598
599
600
601
602
603

    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);
604
      TOut tmp; nullify(tmp);
605
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
606
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq]);
607
608
609
610
611
612
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

613
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
614
    Parallel::mpi::globalAdd(value);
615
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
616
617
618

    return value;
  }
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
651
652
653
654
655
656
657
658
  
  
  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
659

660
661
662
663
664
665
666
667
668
      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
669
    Parallel::mpi::globalAdd(value);
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
#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) {
704
705
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, nullptr, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, nullptr, qp2);
706
707
708
709
710
711
712
713
714
715

      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
716
    Parallel::mpi::globalAdd(value);
717
718
719
720
721
722
#endif

    return value;
  }
  
  
Praetorius, Simon's avatar
Praetorius, Simon committed
723
  template<typename T1, typename T2>
724
725
726
  typename traits::mult_type<T1,T2>::type 
  integrate_VecTimesCoords(const DOFVector<T1> &vec,
			   AbstractFunction<T2, WorldVector<double> > *fct)
727
728
729
730
731
  {
    FUNCNAME("integrate_VecTimesCoords()");
    
    TEST_EXIT(fct)("No function defined!\n");
    
732
    typedef typename traits::mult_type<T1,T2>::type TOut;
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

    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
764
    Parallel::mpi::globalAdd(value);
765
766
767
768
769
770
771
772
#endif

    return value;
  }

  
  template<typename TOut, typename T>
  TOut integrate_VecAndCoords(const DOFVector<T> &vec,
773
			   BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct)
Praetorius, Simon's avatar
Praetorius, Simon committed
774
775
776
777
778
779
780
781
782
783
784
785
  {
    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);

786
787
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());
    WorldVector<double> coordsAtQP;
Praetorius, Simon's avatar
Praetorius, Simon committed
788

789
790
    TOut value; nullify(value);
    
Praetorius, Simon's avatar
Praetorius, Simon committed
791
792
793
794
795
    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);
796
      TOut tmp; nullify(tmp);
Praetorius, Simon's avatar
Praetorius, Simon committed
797
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
798
799
	elInfo->coordToWorld(fastQuad->getLambda(iq), coordsAtQP);
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq], coordsAtQP);
Praetorius, Simon's avatar
Praetorius, Simon committed
800
801
802
803
804
805
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

806
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
807
    Parallel::mpi::globalAdd(value);
808
#endif
809

810
811
    return value;
  }
812

813
  
814
815
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
816
  {  
817
818
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
819
    if (!q) {
820
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
821
822
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
823

824
825
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
826

827
828
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
829
    mtl::dense_vector<T> uh_vec(nPoints);
830
831
832
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
833
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
834
835
836
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
837
      this->getVecAtQPs(elInfo, nullptr, quadFast, uh_vec);
838
839
840
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
841

842
843
      elInfo = stack.traverseNext(elInfo);
    }
844

845
846
847
848
849
850
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
851
852
  }

853

854
855
856
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
857
858
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
859
    if (!q) {
860
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
861
862
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
863

864
865
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
866

867
868
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
869
    mtl::dense_vector<T> uh_vec(nPoints);
870
871
872
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
873
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
874
875
876
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
877
      this->getVecAtQPs(elInfo, nullptr, quadFast, uh_vec);
878
879
880
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;
881

882
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
883
    }
884

885
886
887
888
889
890
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
891
892
  }

893

894
  template<typename T>  
895
896
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
897
898
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
899
    if (!q) {
900
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
901
902
903
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

904
905
906
907
908
909
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
910
    mtl::dense_vector<WorldVector<T> > grduh_vec(nPoints);
911
912
913
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
914
915
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
916
917
918
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
919
      this->getGrdAtQPs(elInfo, nullptr, quadFast, grduh_vec);
920
921
922
923
924
925
926

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

928
      result += det * normT;
929

930
931
      elInfo = stack.traverseNext(elInfo);
    }
932

933
934
935
936
937
938
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
939
940
  }

941

942
943

  template<typename T>
944
  bool DOFVector<T>::getDofIdxAtPoint(WorldVector<double> &p, 
Thomas Witkowski's avatar
Thomas Witkowski committed
945
946
947
948
					    DegreeOfFreedom &idx, 
					    ElInfo *oldElInfo, 
					    bool useOldElInfo) const
  { 
Praetorius, Simon's avatar
Praetorius, Simon committed
949
950
    Mesh *mesh = this->feSpace->getMesh();
    const BasisFunction *basFcts = this->feSpace->getBasisFcts();
951
952
953
954
955
956
957
958
959

    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
960
    bool inside = false;
961
962

    if (oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
963
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), nullptr, nullptr);
964
965
      delete oldElInfo;
    } else {
966
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
967
968
969
970
971
    }
    
    if (oldElInfo)
      oldElInfo = elInfo;
    
Nestler, Michael's avatar
Nestler, Michael committed
972
973
    if (!inside){
      delete elInfo; 
974
      return false;