DOFVector.hh 57 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
#include <list>
Thomas Witkowski's avatar
Thomas Witkowski committed
23
#include <algorithm>
24
#include <math.h>
25

26
27
28
29
30
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/utility/category.hpp>
#include <boost/numeric/linear_algebra/identity.hpp>

31
32
33
34
35
36
37
38
39
40
41
42
#include "FixVec.h"
#include "Boundary.h"
#include "DOFAdmin.h"
#include "ElInfo.h"
#include "Error.h"
#include "FiniteElemSpace.h"
#include "Global.h"
#include "Mesh.h"
#include "Quadrature.h"
#include "AbstractFunction.h"
#include "BoundaryManager.h"
#include "Assembler.h"
43
#include "Operator.h"
44
#include "Initfile.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
45
#include "Traverse.h"
46
47
48
#include "DualTraverse.h"

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
49
#include "parallel/MpiHelper.h"
50
#endif
51

52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
// Defining the interface for MTL4
namespace mtl {

  // Let MTL4 know that DOFVector it is a column vector
  namespace traits {
    template <typename T>
    struct category< AMDiS::DOFVector<T> > 
    {
      typedef tag::dense_col_vector type;
    };
  }

  namespace ashape {
    template <typename T>
    struct ashape< AMDiS::DOFVector<T> > 
    {
      typedef cvec<typename ashape<T>::type> type;
    };
  }

  // Modelling Collection and MutableCollection
  template <typename T>
  struct Collection< AMDiS::DOFVector<T> > 
  {
    typedef T               value_type;
    typedef const T&        const_reference;
    typedef std::size_t     size_type;
  };

  template <typename T>
  struct MutableCollection< AMDiS::DOFVector<T> > 
    : public Collection< AMDiS::DOFVector<T> > 
  {
    typedef T&              reference;
  };


} // namespace mtl



93
94
95
namespace AMDiS {

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

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

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

121

122
  template<typename T>
123
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
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
131
132
133
134
  }

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

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

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

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

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

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

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

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

    return sqrt(nrm);
188
189
  }

190

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

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

207
208
209
    return nrm;
  }

210

211
212
213
  template<typename T>
  T DOFVector<T>::asum() const
  {
214
    checkFeSpace(this->feSpace, vec);
215

216
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
217
218
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
219
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
220
221
      nrm += abs(*vecIterator);

222
223
224
225
226
227
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
228
229
  }

230

231
232
233
  template<typename T>
  T DOFVector<T>::sum() const
  {
234
    checkFeSpace(this->feSpace, vec);
235

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

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

    return nrm;
248
249
  }

250

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

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


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

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

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


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

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

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

298
299
300
    return m;
  }

301

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

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

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

317
    return m;
318
319
  }

320

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

327
328
329

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

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

348
349
350
351
    return m / count;
  }


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

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

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

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

403

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

    return result;
  }
413

414

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

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

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

431
432
433
    return val;
  }

434

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

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

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

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

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

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

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

486

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

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

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

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

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

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

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

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

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

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

    return value;
  }

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

    return value;
  }
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
  
  
  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
649

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

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

    return value;
  }
  
  
Praetorius, Simon's avatar
Praetorius, Simon committed
713
  template<typename T1, typename T2>
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
  typename ProductType<T1,T2>::type integrate_VecTimesCoords(const DOFVector<T1> &vec,
		   AbstractFunction<T2, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate_VecTimesCoords()");
    
    TEST_EXIT(fct)("No function defined!\n");
    
    typedef typename ProductType<T1,T2>::type TOut;

    const FiniteElemSpace *feSpace = vec.getFeSpace();
    Mesh *mesh = feSpace->getMesh();

    int deg = std::max(fct->getDegree(), 2 * feSpace->getBasisFcts()->getDegree());
    Quadrature* quad = Quadrature::provideQuadrature(mesh->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<T1> qp(fastQuad->getNumPoints());
    WorldVector<double> coords;

    TOut value; nullify(value);
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);

      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
 	tmp += fastQuad->getWeight(iq) * (qp[iq] * (*fct)(coords));
      }

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

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
753
    Parallel::mpi::globalAdd(value);
754
755
756
757
758
759
760
761
#endif

    return value;
  }

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

799
800
    return value;
  }
801

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

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

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

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

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

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

842

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

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

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

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

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

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

882

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

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

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

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

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

917
      result += det * normT;
918

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

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

930

931
932

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

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

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


990
991
  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
992
					std::vector<DegreeOfFreedom> &newDOF)
993
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
994
995
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
996
	vec[newDOF[i]] = vec[i];
997
998
  }

999

1000
1001
1002
1003
  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1004
1005
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
1006
	 op != this->operators.end(); ++op)
1007
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
1008

1009
1010
1011
    return fillFlag;
  }

1012

Thomas Witkowski's avatar
Thomas Witkowski committed
1013
1014
1015
1016
1017
  template<typename T>
  void DOFVectorBase<T>::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = this->operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1018
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
1019
1020
1021
      (*it)->finishAssembling();
  }

1022

1023
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
1024
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
1025
  {
1026
    this->feSpace = rhs.feSpace;
1027
    this->dim = rhs.dim;
1028
    this->nBasFcts = rhs.nBasFcts;
1029
    vec = rhs.vec;
1030
    this->elementVector.change_dim(this->nBasFcts);
1031
1032
1033
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
1034
    
1035
    if (rhs.boundaryManager) {
1036
1037
1038
      if (this->boundaryManager) 
	delete this->boundaryManager; 

1039
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
1040
    } else {
1041
      this->boundaryManager = nullptr;