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
    
    #pragma omp critical
147 148
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
149

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

153
    vec.clear();
154 155
  }

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

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

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

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

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

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

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

199
    return std::sqrt(nrm);
200 201
  }

202

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

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

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

219 220 221
    return nrm;
  }

222

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

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

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

    return nrm;
240 241
  }

242

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

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

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

    return nrm;
260 261
  }

262

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

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


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

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

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

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


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

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

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

310 311 312
    return m;
  }

313

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

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

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

329
    return m;
330 331
  }

332

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

339 340 341

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

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

360 361 362 363
    return m / count;
  }


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

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

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

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

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

415

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

    return result;
  }
425

426

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

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

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

443 444 445
    return val;
  }

446

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

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

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

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

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

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

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

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

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

498

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

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

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

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

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

531 532
      elInfo = stack.traverseNext(elInfo);
    }
533

534
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
535
    Parallel::mpi::globalAdd(result);
536
#endif
537

538
    return result;
539 540
  }

541

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

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

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

    WorldVector<double> coords;

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

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

571 572 573 574
      elInfo = stack.traverseNext(elInfo);
    }

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

    return value;
  }

581

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

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

597 598 599
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);
600 601 602 603 604 605

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

      elInfo = stack.traverseNext(elInfo);
    }

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

    return value;
  }
621 622


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

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

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

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

667 668 669 670
      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
671
    Parallel::mpi::globalAdd(value);
672
#endif
673

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

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

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

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

714 715 716 717
      cont = dualTraverse.traverseNext(dualElInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
718
    Parallel::mpi::globalAdd(value);
719 720 721 722
#endif

    return value;
  }
723 724


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

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

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

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

762 763 764 765
      elInfo = stack.traverseNext(elInfo);
    }

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

    return value;
  }

772

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

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

791
    TOut value; nullify(value);
792

Praetorius, Simon's avatar
Praetorius, Simon committed
793 794 795 796 797
    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);