DOFVector.hh 60.8 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * 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.
18
 *
19
 ******************************************************************************/
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
// Defining the interface for MTL4
namespace mtl {

  // Let MTL4 know that DOFVector it is a column vector
  namespace traits {
    template <typename T>
59
    struct category< AMDiS::DOFVector<T> >
60
61
62
63
64
65
66
    {
      typedef tag::dense_col_vector type;
    };
  }

  namespace ashape {
    template <typename T>
67
    struct ashape< AMDiS::DOFVector<T> >
68
69
70
71
72
73
74
    {
      typedef cvec<typename ashape<T>::type> type;
    };
  }

  // Modelling Collection and MutableCollection
  template <typename T>
75
  struct Collection< AMDiS::DOFVector<T> >
76
77
78
79
80
81
82
  {
    typedef T               value_type;
    typedef const T&        const_reference;
    typedef std::size_t     size_type;
  };

  template <typename T>
83
84
  struct MutableCollection< AMDiS::DOFVector<T> >
    : public Collection< AMDiS::DOFVector<T> >
85
86
87
88
89
90
91
92
93
  {
    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(NULL)
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
    vec.resize(0);
118
    init(f, n, addToSynch);
119
  }
120

121

122
  template<typename T>
123
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n, bool addToSynch)
124
  {
125
126
    this->name = n;
    this->feSpace = f;
127
    
128
129
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
130
    
131
    this->boundaryManager = new BoundaryManager(f);
132
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
133
    if (addToSynch && Parallel::MeshDistributor::globalMeshDistributor != NULL)
134
135
      Parallel::MeshDistributor::globalMeshDistributor->addInterchangeVector(this);
#endif
136
137
138
139
140
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
141
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
142
    if ( Parallel::MeshDistributor::globalMeshDistributor != NULL)
143
144
      Parallel::MeshDistributor::globalMeshDistributor->removeInterchangeVector(this);
#endif
145
    
146
    #ifdef _OPENMP
147
    #pragma omp critical
148
    #endif
149
150
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
151

Thomas Witkowski's avatar
Thomas Witkowski committed
152
    if (this->boundaryManager)
153
      delete this->boundaryManager;
154

155
    vec.clear();
156
157
  }

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

170
    for (int i = 0; i < nBasFcts; i++) {
171
      BoundaryCondition *condition =
172
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
173

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

177
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
178
	  (*this)[irow] += factor * elVec[i];
179
	else
180
	  (*this)[irow] = factor * elVec[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
181
182
183
184
      }
    }
  }

185
186
187
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
188
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
189

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

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

201
    return std::sqrt(nrm);
202
203
  }

204

205
206
207
  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
208
    checkFeSpace(this->feSpace, vec);
209

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

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

221
222
223
    return nrm;
  }

224

225
226
227
  template<typename T>
  T DOFVector<T>::asum() const
  {
228
    checkFeSpace(this->feSpace, vec);
229

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

236
237
238
239
240
241
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
242
243
  }

244

245
246
247
  template<typename T>
  T DOFVector<T>::sum() const
  {
248
    checkFeSpace(this->feSpace, vec);
249

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

256
257
258
259
260
261
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
262
263
  }

264

265
266
267
  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
268
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
269

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


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

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

283
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >=
284
		  this->feSpace->getAdmin()->getUsedSize())
285
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(),
286
       this->feSpace->getAdmin()->getUsedSize());
287

288
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
289
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)),
Thomas Witkowski's avatar
Thomas Witkowski committed
290
291
		       USED_DOFS);
    for (vecIterator.reset(), xIterator.reset(); !vecIterator.end();
292
293
	 ++vecIterator, ++xIterator)
      *vecIterator = *xIterator;
294
295
296
297
298
  }


  template<typename T>
  T DOFVector<T>::min() const
299
  {
300
    checkFeSpace(this->feSpace, vec);
301

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

307
308
309
310
311
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMin = m;
    MPI::COMM_WORLD.Allreduce(&localMin, &m, 1, MPI_DOUBLE, MPI_MIN);
#endif

312
313
314
    return m;
  }

315

316
  template<typename T>
317
318
  T DOFVector<T>::max() const
  {
319
    checkFeSpace(this->feSpace, vec);
320

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

326
327
328
329
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMax = m;
    MPI::COMM_WORLD.Allreduce(&localMax, &m, 1, MPI_DOUBLE, MPI_MAX);
#endif
330

331
    return m;
332
333
  }

334

Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
337
338
339
340
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

341
342
343

  template<typename T>
  T DOFVector<T>::average() const
344
  {
345
346
347
348
349
350
351
352
353
354
    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++;
    }

355
356
357
358
359
360
361
#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

362
363
364
365
    return m / count;
  }


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

371
    const DOFAdmin *admin = NULL;
372
    const char *format;
373

374
    if (this->feSpace)
375
      admin = this->feSpace->getAdmin();
376

377
    MSG("Vec `%s':\n", this->name.c_str());
378
    int j = 0;
379
380
381
382
383
384
385
386
387
    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
388
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
389
	if ((j % 3) == 0) {
390
	  if (j)
391
	    Msg::print("\n");
392
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
393
	} else {
394
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
395
	}
396
	j++;
397
      }
398
      Msg::print("\n");
399
    } else {
400
	MSG("no DOFAdmin, print whole vector.\n");
401

402
403
	for (int i = 0; i < static_cast<int>( vec.size()); i++) {
	  if ((j % 3) == 0) {
404
	    if (j)
405
406
407
	      Msg::print("\n");
	    MSG("(%d,%10.5e)",i,vec[i]);
	  } else {
408
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
409
	  }
410
411
412
413
414
415
416
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }

417

418
  template<typename T>
419
  size_t DOFVector<T>::calcMemoryUsage() const
420
  {
421
    size_t result = 0;
422
423
424
425
426
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
427

428

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

437
    for (int i = 0; i < nBasisFcts; i++)
438
439
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

440
441
442
443
444
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localVal = val;
    MPI::COMM_WORLD.Allreduce(&localVal, &val, 1, MPI_DOUBLE, MPI_SUM);
#endif

445
446
447
    return val;
  }

448

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

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

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

462
    if (!(this->getFeSpace()->getAdmin())) {
463
      MSG("no dof admin in feSpace %s, skipping interpolation\n",
464
	  this->getFeSpace()->getName().c_str());
465
466
      return;
    }
467

468
    if (!(this->getFeSpace()->getBasisFcts())) {
469
470
471
472
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
473

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

480
481
482
483
    const BasisFunction *basFct = this->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = this->getFeSpace()->getAdmin();
    int nBasFcts = basFct->getNumber();
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
484
    mtl::dense_vector<T> fctInterpolValues(nBasFcts);
485

Thomas Witkowski's avatar
Thomas Witkowski committed
486
    TraverseStack stack;
487
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
488
489
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
490
      basFct->interpol(elInfo, 0, NULL, fct, fctInterpolValues);
491
      basFct->getLocalIndices(const_cast<Element*>(elInfo->getElement()),
492
493
494
495
			      admin, myLocalIndices);
      for (int i = 0; i < nBasFcts; i++)
	vec[myLocalIndices[i]] = fctInterpolValues[i];

Thomas Witkowski's avatar
Thomas Witkowski committed
496
497
      elInfo = stack.traverseNext(elInfo);
    }
498
499
  }

500

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

Thomas Witkowski's avatar
Thomas Witkowski committed
506
    if (!q) {
507
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
508
509
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
510

511
    FastQuadrature *quadFast =
512
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
513

Thomas Witkowski's avatar
Thomas Witkowski committed
514
515
516
517
518
519
    Flag flag = Mesh::FILL_COORDS | Mesh::FILL_DET;
    if (meshLevel == -1)
      flag |= Mesh::CALL_LEAF_EL;
    else
      flag |= Mesh::CALL_EL_LEVEL;

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

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

536
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
537
    Parallel::mpi::globalAdd(result);
538
#endif
539

540
    return result;
541
542
  }

543

544
545
546
  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
547
  {
548
    FUNCNAME("integrate_Coords()");
549

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

    int deg = std::max(fct->getDegree(), feSpace->getBasisFcts()->getDegree());
553
    Quadrature* quad =
554
555
556
557
558
559
560
      Quadrature::provideQuadrature(feSpace->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);

    WorldVector<double> coords;

    TOut value; nullify(value);
561
562
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
563
564
565
566
567
    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);
568
 	tmp += (*fct)(coords) * fastQuad->getWeight(iq);
569
570
571
      }

      value += tmp * elInfo->getDet();
572

573
574
575
576
      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
577
    Parallel::mpi::globalAdd(value);
578
579
580
581
582
#endif

    return value;
  }

583

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

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

599
600
601
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);
602
603
604
605
606
607

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

      elInfo = stack.traverseNext(elInfo);
    }

617
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
618
    Parallel::mpi::globalAdd(value);
619
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
620
621
622

    return value;
  }
623
624


625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
  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()");
643

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

646
    int deg = std::max(fct->getDegree(),
647
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
648
    Quadrature* quad =
649
650
651
652
653
654
655
656
      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);
657
658
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
659
660
661
662
    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
663

664
665
      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
666
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
667
      value += tmp * elInfo->getDet();
668

669
670
671
672
      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
673
    Parallel::mpi::globalAdd(value);
674
#endif
675

676
677
678
679
680
681
682
683
684
685
    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()");
686

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

689
    int deg = std::max(fct->getDegree(),
690
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
691
    Quadrature* quad =
692
693
694
695
696
697
698
699
      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);
700
701
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    DualTraverse dualTraverse;
702
703
704
705
706
707
    DualElInfo dualElInfo;
    bool cont = dualTraverse.traverseFirst(vec1.getFeSpace()->getMesh(),
					   vec2.getFeSpace()->getMesh(),
					   -1, -1, traverseFlag, traverseFlag,
					   dualElInfo);
    while (cont) {
708
709
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);
710
711
712

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

716
717
718
719
      cont = dualTraverse.traverseNext(dualElInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
720
    Parallel::mpi::globalAdd(value);
721
722
723
724
#endif

    return value;
  }
725
726


Praetorius, Simon's avatar
Praetorius, Simon committed
727
  template<typename T1, typename T2>
728
  typename traits::mult_type<T1,T2>::type
729
730
  integrate_VecTimesCoords(const DOFVector<T1> &vec,
			   AbstractFunction<T2, WorldVector<double> > *fct)
731
732
  {
    FUNCNAME("integrate_VecTimesCoords()");
733

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

736
    typedef typename traits::mult_type<T1,T2>::type TOut;
737
738
739
740
741
742
743
744
745
746
747
748
749

    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);
750
751
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
752
753
754
755
756
757
758
759
760
761
762
    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();
763

764
765
766
767
      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
768
    Parallel::mpi::globalAdd(value);
769
770
771
772
773
#endif

    return value;
  }

774

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

790
791
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());
    WorldVector<double> coordsAtQP;
Praetorius, Simon's avatar
Praetorius, Simon committed
792

793
    TOut value; nullify(value);
794

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

      elInfo = stack.traverseNext(elInfo);
    }

810
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
811
    Parallel::mpi::globalAdd(value);
812
#endif
813

814
815
    return value;
  }
816

817

818
819
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
820
  {
821
822
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
823
    if (!q) {
824
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
825
826
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
827

828
    FastQuadrature *quadFast =
829
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
830

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

846
847
      elInfo = stack.traverseNext(elInfo);
    }
848

849
850
851
852
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
853

854
    return result;
855
856
  }

857

858
859
860
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
861
862
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
863
    if (!q) {
864
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
865
866
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
867

868
    FastQuadrature *quadFast =
869
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
870

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

886
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
887
    }
888

889
890
891
892
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
893

894
    return result;
895
896
  }

897