DOFVector.hh 56.5 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
					  bool add)
  {
142
    std::vector<DegreeOfFreedom> indices(nBasFcts);
143
144
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), 
					     feSpace->getAdmin(),
145
					     indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
146

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

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

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

162
163
164
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
165
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
166

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

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

    return sqrt(nrm);
179
180
  }

181

182
183
184
  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
185
    checkFeSpace(this->feSpace, vec);
186
187
    
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
188
189
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
190
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
191
192
      nrm += (*vecIterator) * (*vecIterator);

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

198
199
200
    return nrm;
  }

201

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

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

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

    return nrm;
219
220
  }

221

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

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

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

    return nrm;
239
240
  }

241

242
243
244
  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
245
    checkFeSpace(this->feSpace, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
246

247
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
248
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
249
250
251
252
253
254
255
      *vecIterator = alpha ; 
  }


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

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

260
261
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >= 
		  this->feSpace->getAdmin()->getUsedSize())
262
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
263
       this->feSpace->getAdmin()->getUsedSize());
264
    
265
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
268
269
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), 
		       USED_DOFS);
    for (vecIterator.reset(), xIterator.reset(); !vecIterator.end();
	 ++vecIterator, ++xIterator)      
270
      *vecIterator = *xIterator; 
271
272
273
274
275
  }


  template<typename T>
  T DOFVector<T>::min() const
276
  {    
277
    checkFeSpace(this->feSpace, vec);
278
279

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

284
285
286
287
288
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMin = m;
    MPI::COMM_WORLD.Allreduce(&localMin, &m, 1, MPI_DOUBLE, MPI_MIN);
#endif

289
290
291
    return m;
  }

292

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

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

303
304
305
306
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMax = m;
    MPI::COMM_WORLD.Allreduce(&localMax, &m, 1, MPI_DOUBLE, MPI_MAX);
#endif
307

308
    return m;
309
310
  }

311

Thomas Witkowski's avatar
Thomas Witkowski committed
312
313
314
315
316
317
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

318
319
320

  template<typename T>
  T DOFVector<T>::average() const
321
  {    
322
323
324
325
326
327
328
329
330
331
    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++;
    }

332
333
334
335
336
337
338
#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

339
340
341
342
    return m / count;
  }


343
344
345
  template<typename T>
  void DOFVector<T>::print() const
  {
346
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
347

348
    const DOFAdmin *admin = NULL;
349
    const char *format;
350

351
352
    if (this->feSpace) 
      admin = this->feSpace->getAdmin();
353

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

394

395
396
397
398
399
400
401
402
403
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
404

405

406
407
408
409
  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda, 
			     DegreeOfFreedom* dof_indices)
  {
410
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFeSpace()->getBasisFcts());
411
    int nBasisFcts = phi->getNumber();
412
413
    T val = 0.0;

414
    for (int i = 0; i < nBasisFcts; i++)
415
416
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

417
418
419
420
421
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localVal = val;
    MPI::COMM_WORLD.Allreduce(&localVal, &val, 1, MPI_DOUBLE, MPI_SUM);
#endif

422
423
424
    return val;
  }

425

426
427
428
  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
429
    FUNCNAME("DOFVector::interpol()");
430

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

433
    if (!this->getFeSpace()) {
434
435
436
437
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
438

439
    if (!(this->getFeSpace()->getAdmin())) {
440
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
441
	  this->getFeSpace()->getName().c_str());
442
443
      return;
    }
444
  
445
    if (!(this->getFeSpace()->getBasisFcts())) {
446
447
448
449
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
450

451
452
453
454
455
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
456

457
458
459
460
    const BasisFunction *basFct = this->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = this->getFeSpace()->getAdmin();
    int nBasFcts = basFct->getNumber();
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
461
    mtl::dense_vector<T> fctInterpolValues(nBasFcts);
462

Thomas Witkowski's avatar
Thomas Witkowski committed
463
    TraverseStack stack;
464
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
465
466
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
467
468
469
470
471
472
      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
473
474
      elInfo = stack.traverseNext(elInfo);
    }
475
476
  }

477

478
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
479
  double DOFVector<T>::Int(int meshLevel, Quadrature* q) const
480
  {
481
    Mesh* mesh = this->feSpace->getMesh();
482

Thomas Witkowski's avatar
Thomas Witkowski committed
483
    if (!q) {
484
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
487

488
    FastQuadrature *quadFast = 
489
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
490

Thomas Witkowski's avatar
Thomas Witkowski committed
491
492
493
494
495
496
    Flag flag = Mesh::FILL_COORDS | Mesh::FILL_DET;
    if (meshLevel == -1)
      flag |= Mesh::CALL_LEAF_EL;
    else
      flag |= Mesh::CALL_EL_LEVEL;

497
    T result; nullify(result);
498
    int nPoints = quadFast->getNumPoints();
499
    mtl::dense_vector<T> uh_vec(nPoints);
500
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
501
    ElInfo *elInfo =  stack.traverseFirst(mesh, meshLevel, flag);
502
503
    while (elInfo) {
      double det = elInfo->getDet();
504
      T normT; nullify(normT);
505
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
506
507
508
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;
509

510
511
      elInfo = stack.traverseNext(elInfo);
    }
512

513
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
514
    mpi::globalAdd(result);
515
516
517
#endif
    
    return result;
518
519
  }

520
521
522
523
  
  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
524
  {
525
526
527
528
529
530
531
532
533
534
535
536
537
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
    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()");
566
567
568
569
570
571
572
573
574
575

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

576
577
578
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);
579
580
581
582
583
584

    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);
585
      TOut tmp; nullify(tmp);
586
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
587
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq]);
588
589
590
591
592
593
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

594
595
596
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
597
598
599

    return value;
  }
600
601
602
603
604
605
606
607
608
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
  
  
  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
640

641
642
643
644
645
646
647
648
649
650
651
652
653
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
      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
704
  template<typename T1, typename T2>
705
706
707
708
709
710
711
712
713
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
    mpi::globalAdd(value);
#endif

    return value;
  }

  
  template<typename TOut, typename T>
  TOut integrate_VecAndCoords(const DOFVector<T> &vec,
753
			   BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct)
Praetorius, Simon's avatar
Praetorius, Simon committed
754
755
756
757
758
759
760
761
762
763
764
765
  {
    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);

766
767
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());
    WorldVector<double> coordsAtQP;
Praetorius, Simon's avatar
Praetorius, Simon committed
768

769
770
    TOut value; nullify(value);
    
Praetorius, Simon's avatar
Praetorius, Simon committed
771
772
773
774
775
    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);
776
      TOut tmp; nullify(tmp);
Praetorius, Simon's avatar
Praetorius, Simon committed
777
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
778
779
	elInfo->coordToWorld(fastQuad->getLambda(iq), coordsAtQP);
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq], coordsAtQP);
Praetorius, Simon's avatar
Praetorius, Simon committed
780
781
782
783
784
785
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

786
787
788
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    mpi::globalAdd(value);
#endif
789

790
791
    return value;
  }
792

793
  
794
795
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
796
  {  
797
798
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
799
    if (!q) {
800
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
801
802
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
803

804
805
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
806

807
808
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
809
    mtl::dense_vector<T> uh_vec(nPoints);
810
811
812
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
813
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
814
815
816
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
817
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
818
819
820
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;
821

822
823
      elInfo = stack.traverseNext(elInfo);
    }
824

825
826
827
828
829
830
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
831
832
  }

833

834
835
836
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
837
838
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
839
    if (!q) {
840
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
841
842
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
843

844
845
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
846

847
848
    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
849
    mtl::dense_vector<T> uh_vec(nPoints);
850
851
852
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
853
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
854
855
856
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
857
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
858
859
860
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;
861

862
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
863
    }
864

865
866
867
868
869
870
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
871
872
  }

873

874
  template<typename T>  
875
876
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
877
878
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
879
    if (!q) {
880
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
881
882
883
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

884
885
886
887
888
889
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
890
    mtl::dense_vector<WorldVector<T> > grduh_vec(nPoints);
891
892
893
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
894
895
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
896
897
898
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
899
      this->getGrdAtQPs(elInfo, NULL, quadFast, grduh_vec);
900
901
902
903
904
905
906

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

908
      result += det * normT;
909

910
911
      elInfo = stack.traverseNext(elInfo);
    }
912

913
914
915
916
917
918
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
919
920
  }

921

922
923

  template<typename T>
924
  bool DOFVector<T>::getDofIdxAtPoint(WorldVector<double> &p, 
Thomas Witkowski's avatar
Thomas Witkowski committed
925
926
927
928
					    DegreeOfFreedom &idx, 
					    ElInfo *oldElInfo, 
					    bool useOldElInfo) const
  { 
Praetorius, Simon's avatar
Praetorius, Simon committed
929
930
    Mesh *mesh = this->feSpace->getMesh();
    const BasisFunction *basFcts = this->feSpace->getBasisFcts();
931
932
933
934
935
936
937
938
939

    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
940
    bool inside = false;
941
942
943
944
945
946
947
948
949
950
951
952
953
954

    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
955
    basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices);
956
957
    
    WorldVector<double> coord;
958
959
960
    int minIdx = -1;
    double minDist = 1.e15;
    
961
962
963
    for (int i = 0; i < numBasFcts; i++) {
      elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
      WorldVector<double> dist = coord - p;
964
965
966
967
      double newDist = norm(dist);
      if (newDist < minDist) {
	minIdx = i;
	minDist = newDist;
968
969
	if (minDist < 1.e-10)
	  break;
970
971
972
973
974
975
976
977
      }
    }
    
    if (minIdx >= 0)
      idx = localIndices[minIdx];
    
    if(!oldElInfo) delete elInfo;
    return inside;
Thomas Witkowski's avatar
Thomas Witkowski committed
978
  }
979
980


981
982
  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
983
					std::vector<DegreeOfFreedom> &newDOF)
984
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
985
986
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
987
	vec[newDOF[i]] = vec[i];
988
989
  }

990

991
992
993
994
  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
995
996
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
997
	 op != this->operators.end(); ++op)
998
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
999

1000
1001
1002
    return fillFlag;
  }

1003

Thomas Witkowski's avatar
Thomas Witkowski committed
1004
1005
1006
1007
1008
  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
1009
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
1010
1011
1012
      (*it)->finishAssembling();
  }

1013

1014
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
1015
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
1016
  {
1017
    this->feSpace = rhs.feSpace;
1018
    this->dim = rhs.dim;
1019
    this->nBasFcts = rhs.nBasFcts;
1020
    vec = rhs.vec;
1021
    this->elementVector.change_dim(this->nBasFcts);
1022
1023
1024
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
1025
    
1026
    if (rhs.boundaryManager) {
1027
1028
1029
      if (this->boundaryManager) 
	delete this->boundaryManager; 

1030
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
1031
    } else {
1032
      this->boundaryManager = NULL;
1033
1034
    }

1035
1036
1037
    return *this;
  }

1038

1039
1040
1041
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal)
  {
1042
    FUNCNAME_DBG("DOFVector<T>::operator*=(DOFVector<T>& x, T scal)");
1043

1044
1045
    TEST_EXIT_DBG(x.getFeSpace() && x.getFeSpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFeSpace(), x.getFeSpace()->getAdmin());
1046

Thomas Witkowski's avatar
Thomas Witkowski committed
1047
1048
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
1049
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
1050
      (*vecIterator) *= scal; 
1051

1052
1053
1054
1055
1056
1057
1058
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
1059
    FUNCNAME_DBG("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
1060
    
1061
1062
1063
1064
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
1065
      ("no admin or different admins: %8X, %8X\n",
1066
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
1067
1068
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");