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
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
    const BasisFunction *basFct = this->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = this->getFeSpace()->getAdmin();
    int nBasFcts = basFct->getNumber();
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
478
    mtl::dense_vector<T> fctInterpolValues(nBasFcts);
479

Thomas Witkowski's avatar
Thomas Witkowski committed
480
    TraverseStack stack;
481
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Thomas Witkowski's avatar
Thomas Witkowski committed
482
483
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
484
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
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
496
  double DOFVector<T>::Int(int meshLevel, Quadrature* q) const
497
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
498
    FUNCNAME("DOFVector::Int()");
Thomas Witkowski's avatar
Thomas Witkowski committed
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

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

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

529
530
      elInfo = stack.traverseNext(elInfo);
    }
531

532
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
533
    mpi::globalAdd(result);
534
535
536
#endif
    
    return result;
537
538
  }

539
540
541
542
  
  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
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
580
581
582
583
584
    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()");
585
586
587
588
589
590
591
592
593
594

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

595
596
597
    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);
598
599
600
601
602
603

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

      elInfo = stack.traverseNext(elInfo);
    }

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

    return value;
  }
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
  
  
  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
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
718
719
720
721
722
      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
723
  template<typename T1, typename T2>
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
768
769
770
771
  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,
772
			   BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct)
Praetorius, Simon's avatar
Praetorius, Simon committed
773
774
775
776
777
778
779
780
781
782
783
784
  {
    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);

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

809
810
    return value;
  }
811

812
  
813
814
815
  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
816
    FUNCNAME("DOFVector::L1Norm()");
817
  
818
819
    Mesh* mesh = this->feSpace->getMesh();

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

825
826
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
827

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

843
844
      elInfo = stack.traverseNext(elInfo);
    }
845

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

854

855
856
857
  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
858
    FUNCNAME("DOFVector::L2NormSquare()");
859

860
861
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
862
    if (!q) {
863
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
864
865
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
866

867
868
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);
869

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

885
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
886
    }
887

888
889
890
891
892
893
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
894
895
  }

896

897
  template<typename T>  
898
899
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
900
    FUNCNAME("DOFVector::H1NormSquare()");
901

902
903
    Mesh* mesh = this->feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
904
    if (!q) {
905
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
906
907
908
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

909
910
911
912
913
914
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

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

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

933
      result += det * normT;
934

935
936
      elInfo = stack.traverseNext(elInfo);
    }
937

938
939
940
941
942
943
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
944
945
  }

946

947
948

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
949
950
951
952
953
954
  const bool DOFVector<T>::getDofIdxAtPoint(WorldVector<double> &p, 
					    DegreeOfFreedom &idx, 
					    ElInfo *oldElInfo, 
					    bool useOldElInfo) const
  { 
    FUNCNAME("DOFVector::getDofIdxAtPoint()");
955

Praetorius, Simon's avatar
Praetorius, Simon committed
956
957
    Mesh *mesh = this->feSpace->getMesh();
    const BasisFunction *basFcts = this->feSpace->getBasisFcts();
958
959
960
961
962
963
964
965
966

    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
967
    bool inside = false;
968
969
970
971
972
973
974
975
976
977
978
979
980
981

    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
982
    basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices);
983
984
    
    WorldVector<double> coord;
985
986
987
    int minIdx = -1;
    double minDist = 1.e15;
    
988
989
990
    for (int i = 0; i < numBasFcts; i++) {
      elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
      WorldVector<double> dist = coord - p;
991
992
993
994
      double newDist = norm(dist);
      if (newDist < minDist) {
	minIdx = i;
	minDist = newDist;
995
996
	if (minDist < 1.e-10)
	  break;
997
998
999
1000
1001
1002
1003
1004
      }
    }
    
    if (minIdx >= 0)
      idx = localIndices[minIdx];
    
    if(!oldElInfo) delete elInfo;
    return inside;
Thomas Witkowski's avatar
Thomas Witkowski committed
1005
  }
1006
1007


1008
1009
  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
1010
					std::vector<DegreeOfFreedom> &newDOF)
1011
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1012
1013
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
1014
	vec[newDOF[i]] = vec[i];
1015
1016
  }

1017

1018
1019
1020
1021
  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1022
1023
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
1024
	 op != this->operators.end(); ++op)
1025
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
1026