DOFVector.hh 57.4 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(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
127
128
129
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
    this->boundaryManager = new BoundaryManager(f);
130
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
131
    if (addToSynch && Parallel::MeshDistributor::globalMeshDistributor != NULL) 
132
133
      Parallel::MeshDistributor::globalMeshDistributor->addInterchangeVector(this);
#endif
134
135
136
137
138
  }

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

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

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

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

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

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

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

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

    return sqrt(nrm);
196
197
  }

198

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

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

215
216
217
    return nrm;
  }

218

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

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

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

    return nrm;
236
237
  }

238

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

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

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

    return nrm;
256
257
  }

258

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

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


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

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

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


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

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

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

306
307
308
    return m;
  }

309

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

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

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

325
    return m;
326
327
  }

328

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

335
336
337

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

349
350
351
352
353
354
355
#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

356
357
358
359
    return m / count;
  }


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

365
    const DOFAdmin *admin = NULL;
366
    const char *format;
367

368
369
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
370

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

411

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

    return result;
  }
421

422

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

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

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

439
440
441
    return val;
  }

442

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

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

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

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

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

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

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

494

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

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

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

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

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

527
528
      elInfo = stack.traverseNext(elInfo);
    }
529

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

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

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

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

    return value;
  }

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

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

593
594
595
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);
596
597
598
599
600
601

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

      elInfo = stack.traverseNext(elInfo);
    }

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

    return value;
  }
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
  
  
  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
657

658
659
660
661
662
663
664
665
666
      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
667
    Parallel::mpi::globalAdd(value);
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
#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) {
702
703
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);
704
705
706
707
708
709
710
711
712
713

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

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

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

    return value;
  }

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

804
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
805
    Parallel::mpi::globalAdd(value);
806
#endif
807

808
809
    return value;
  }
810

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

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

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

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

840
841
      elInfo = stack.traverseNext(elInfo);
    }
842

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

851

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

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

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

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

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

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

891

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

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

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

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

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

926
      result += det * normT;
927

928
929
      elInfo = stack.traverseNext(elInfo);
    }
930

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

939

940
941

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

    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
958
    bool inside = false;
959
960

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