DOFVector.hh 56.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include <list>
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include <algorithm>
15
#include <math.h>
16

17
18
19
20
21
#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>

22
23
24
25
26
27
28
29
30
31
32
33
#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"
34
#include "Operator.h"
35
#include "Initfile.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
36
#include "Traverse.h"
37
38
39
#include "DualTraverse.h"

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
40
#include "parallel/MpiHelper.h"
41
#endif
42

43
44
45
46
47
48
49
50
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
// 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



84
85
86
namespace AMDiS {

  template<typename T>
87
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
88
89
    : feSpace(f),
      name(n),
90
      elementVector(f->getBasisFcts()->getNumber()),
Thomas Witkowski's avatar
Thomas Witkowski committed
91
      boundaryManager(NULL)
92
  {    
Thomas Witkowski's avatar
Thomas Witkowski committed
93
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
94
    dim = feSpace->getMesh()->getDim();
95
96
97
98
99
  }
  

  template<typename T>
  DOFVectorBase<T>::~DOFVectorBase()
100
  {}
101

102
 
Thomas Witkowski's avatar
Thomas Witkowski committed
103
  template<typename T>
104
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
105
    : DOFVectorBase<T>(f, n),
106
      coarsenOperation(COARSE_INTERPOL)
107
  {
108
    vec.resize(0);
109
110
111
    init(f, n);
  } 

112

113
  template<typename T>
114
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
115
  {
116
117
118
119
120
    this->name = n;
    this->feSpace = f;
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
    this->boundaryManager = new BoundaryManager(f);
121
122
123
124
125
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
126
127
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);
128
    
Thomas Witkowski's avatar
Thomas Witkowski committed
129
    if (this->boundaryManager)
130
      delete this->boundaryManager;
131
    
132
    vec.clear();
133
134
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
135
136
137
138
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
139
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
142
143
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

144
145
146
    std::vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
147

148
    for (DegreeOfFreedom i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
149
150
151
      BoundaryCondition *condition = 
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;

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

155
	if (add)
Thomas Witkowski's avatar
Thomas Witkowski committed
156
	  (*this)[irow] += factor * elVec[i];
157
158
	else  
	  (*this)[irow] = factor * elVec[i];	
Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
161
162
      }
    }
  }

163
164
165
166
167
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

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

170
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
173
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
174
175
      nrm += (*vecIterator) * (*vecIterator);

176
177
178
179
180
181
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return sqrt(nrm);
182
183
  }

184

185
186
187
188
189
  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

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

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

203
204
205
    return nrm;
  }

206

207
208
209
  template<typename T>
  T DOFVector<T>::asum() const
  {
210
    FUNCNAME("DOFVector<T>::asum()");
211

212
    checkFeSpace(this->feSpace, vec);
213

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

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

    return nrm;
226
227
  }

228

229
230
231
  template<typename T>
  T DOFVector<T>::sum() const
  {
232
    FUNCNAME("DOFVector<T>::sum()");
233

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
    FUNCNAME("DOFVector<T>::set()");
255

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

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


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
267
    FUNCNAME("DOFVector<T>::copy()");
Thomas Witkowski's avatar
Thomas Witkowski committed
268

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

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


  template<typename T>
  T DOFVector<T>::min() const
  {
288
289
    FUNCNAME("DOFVector<T>::min()");
    
290
    checkFeSpace(this->feSpace, vec);
291
292

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

297
298
299
300
301
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMin = m;
    MPI::COMM_WORLD.Allreduce(&localMin, &m, 1, MPI_DOUBLE, MPI_MIN);
#endif

302
303
304
    return m;
  }

305

306
307
308
  template<typename T>
  T DOFVector<T>::max() const 
  {
309
310
    FUNCNAME("DOFVector<T>::max()");
    
311
    checkFeSpace(this->feSpace, vec);
312

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

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

323
    return m;
324
325
  }

326

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

333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348

  template<typename T>
  T DOFVector<T>::average() const
  {
    FUNCNAME("DOFVector<T>::average()");
    
    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
413
414
415
416
417
418
419
420
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    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
469
470
471
472
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
473

474
475
476
477
478
    const BasisFunction *basFct = this->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = this->getFeSpace()->getAdmin();
    int nBasFcts = basFct->getNumber();
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
    mtl::dense_vector<double> 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
485
486
487
488
489
      basFct->interpol(elInfo, 0, NULL, fct, fctInterpolValues);
      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
496
497
  template<typename T>
  double DOFVector<T>::Int(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
498
    FUNCNAME("DOFVector::Int()");
499
  
500
    Mesh* mesh = this->feSpace->getMesh();
501

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

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

510
    T result; nullify(result);
511
    int nPoints = quadFast->getNumPoints();
512
    mtl::dense_vector<T> uh_vec(nPoints);
513
    TraverseStack stack;
514
515
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
      Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
516
517
    while (elInfo) {
      double det = elInfo->getDet();
518
      T normT; nullify(normT);
519
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
520
521
522
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;
523

524
525
      elInfo = stack.traverseNext(elInfo);
    }
526

527
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
528
    mpi::globalAdd(result);
529
530
531
#endif
    
    return result;
532
533
  }

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

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

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif

    return value;
  }

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

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

590
591
592
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);
593
594
595
596
597
598

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

      elInfo = stack.traverseNext(elInfo);
    }

608
609
610
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
611
612
613

    return value;
  }
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
649
650
651
652
653
  
  
  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
654

655
656
657
658
659
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
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
      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
    mpi::globalAdd(value);
#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) {
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);

      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
    mpi::globalAdd(value);
#endif

    return value;
  }
  
  
Praetorius, Simon's avatar
Praetorius, Simon committed
718
  template<typename T1, typename T2>
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
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
  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
    mpi::globalAdd(value);
#endif

    return value;
  }

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

780
781
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());
    WorldVector<double> coordsAtQP;
Praetorius, Simon's avatar
Praetorius, Simon committed
782

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

      elInfo = stack.traverseNext(elInfo);
    }

800
801
802
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif
803

804
805
    return value;
  }
806

807
  
808
809
810
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
811
    FUNCNAME("DOFVector::L1Norm()");
812
  
813
814
    Mesh* mesh = this->feSpace->getMesh();

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

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

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

838
839
      elInfo = stack.traverseNext(elInfo);
    }
840

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

849

850
851
852
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
853
    FUNCNAME("DOFVector::L2NormSquare()");
854

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
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
895
    FUNCNAME("DOFVector::H1NormSquare()");
896

897
898
    Mesh* mesh = this->feSpace->getMesh();

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

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

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

      for (int iq = 0; iq < nPoints; iq++) {
	double norm2 = 0.0;
	for (int j = 0; j < dimOfWorld; j++)
	  norm2 += sqr(grduh_vec[iq][j]);
	normT += quadFast->getWeight(iq) * norm2;
      }
927

928
      result += det * normT;
929

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

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

941

942
943

  template<typename T>
Praetorius, Simon's avatar
Praetorius, Simon committed
944
  const bool DOFVector<T>::getDOFidxAtPoint(WorldVector<double> &p, DegreeOfFreedom &idx, ElInfo *oldElInfo, bool useOldElInfo) const
945
946
  { FUNCNAME("DOFVector::getDOFidxAtPoint()");

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
961
962
963
964
965
966
967
968
969
970
971
972

    if (oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
      delete oldElInfo;
    } else {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
    }
    
    if (oldElInfo)
      oldElInfo = elInfo;
    
    if (!inside)
      return false;
    
Praetorius, Simon's avatar
Praetorius, Simon committed
973
    basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices);
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
    FixVec<WorldVector<double>, VERTEX> coords = elInfo->getCoords();
    int minIdx = -1;
    double minDist = 1.e15;
    
    for (int i = 0; i < coords.size(); i++) {
      WorldVector<double> dist = coords[i] - p;
      double newDist = norm(dist);
      if (newDist < minDist) {
	minIdx = i;
	minDist = newDist;
      }
    }
    
    if (minIdx >= 0)
      idx = localIndices[minIdx];
    
    if(!oldElInfo) delete elInfo;
    return inside;
  };


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

1004

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

1014
1015
1016
    return fillFlag;
  }

1017

Thomas Witkowski's avatar
Thomas Witkowski committed
1018
1019
1020
1021
1022
  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
1023
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
1024
1025
1026
      (*it)->finishAssembling();
  }

1027

1028
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
1029
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
1030
  {
1031
    this->feSpace = rhs.feSpace;