DOFVector.cc 31.1 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 <boost/numeric/mtl/mtl.hpp>
14
15
16
17
18
19
20
21
22
23
#include "DOFVector.h"
#include "Traverse.h"
#include "DualTraverse.h"
#include "FixVec.h"

namespace AMDiS {

  template<>
  void DOFVector<double>::coarseRestrict(RCNeighbourList& list, int n)
  {
24
25
    FUNCNAME("DOFVector<double>::coarseRestrict()");

26
    switch (coarsenOperation) {
27
28
29
30
31
32
33
    case NO_OPERATION:
      return;
      break;
    case COARSE_RESTRICT:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
      break;
    case COARSE_INTERPOL:
34
35
      TEST_EXIT_DBG(feSpace)("Should not happen!\n");
      TEST_EXIT_DBG(feSpace->getBasisFcts())("Shoud not happen!\n");
36

37
38
39
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
      break;
    default:
40
41
      WARNING("Invalid coarsen operation \"%d\" in vector \"%s\"\n", 
	      coarsenOperation, this->name.c_str());
42
43
44
    }
  }

45

46
47
48
49
50
51
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList& list, int n)
  {
    (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->refineInter(this, &list, n);
  }

52

53
54
55
  template<>
  void DOFVector<WorldVector<double> >::refineInterpol(RCNeighbourList& list, int n)
  {
56
57
58
    if (n < 1) 
      return;

59
    Element *el = list.getElement(0);
60
    int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
61
62
63
    DegreeOfFreedom dof0 = el->getDof(0, n0);
    DegreeOfFreedom dof1 = el->getDof(1, n0);
    DegreeOfFreedom dof_new = el->getChild(0)->getDof(feSpace->getMesh()->getDim(), n0);
64
65
66
67
68
    vec[dof_new] = vec[dof0];
    vec[dof_new] += vec[dof1];
    vec[dof_new] *= 0.5;
  }

69

70
  template<>
Praetorius, Simon's avatar
Praetorius, Simon committed
71
  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* values) const
72
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
73
    FUNCNAME("DOFVector<double>::evalAtCoords()");
74
75
76
77
78
79
80
  
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Praetorius, Simon's avatar
Praetorius, Simon committed
81
    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
82
83
84
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();
85
    static double value = 0.0;
86
87
88
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
Praetorius, Simon's avatar
Praetorius, Simon committed
89
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
90
91
92
93
94
95
96
97
      delete oldElInfo;
    } else
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Praetorius, Simon's avatar
Praetorius, Simon committed
98
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
99
100
101
102
103
104
105
      ElementVector uh(lambda.size());
      for (int i = 0; i < lambda.size(); i++)
        uh[i] = operator[](localIndices[i]);
      value = basFcts->evalUh(lambda, uh);
    } else
      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));

Praetorius, Simon's avatar
Praetorius, Simon committed
106
    delete [] localIndices;
107
108
109
    if (oldElInfo == NULL)
      delete elInfo;

Praetorius, Simon's avatar
Praetorius, Simon committed
110
111
    if (values != NULL)
      *values = value;
112
    return value;
Praetorius, Simon's avatar
Praetorius, Simon committed
113
  };
114
115
116


  template<>
Praetorius, Simon's avatar
Praetorius, Simon committed
117
  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
118
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
119
    FUNCNAME("DOFVector<double>::evalAtCoords()");
120
121
122
123
124
125
126
  
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Praetorius, Simon's avatar
Praetorius, Simon committed
127
    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
128
129
130
131
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();

Praetorius, Simon's avatar
Praetorius, Simon committed
132
133
134
    static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
    WorldVector<double> *val = (NULL != values) ? values : &Values;

135
136
137
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
Praetorius, Simon's avatar
Praetorius, Simon committed
138
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
139
140
141
142
143
144
145
146
      delete oldElInfo;
    } else
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Praetorius, Simon's avatar
Praetorius, Simon committed
147
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
148
149
150
      mtl::dense_vector<WorldVector<double> > uh(lambda.size());
      for (int i = 0; i < lambda.size(); i++)
        uh[i] = operator[](localIndices[i]);
Praetorius, Simon's avatar
Praetorius, Simon committed
151
      *val = basFcts->evalUh(lambda, uh);
152
153
154
    } else
      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));

Praetorius, Simon's avatar
Praetorius, Simon committed
155
    delete [] localIndices;
156
157
158
    if (oldElInfo == NULL)
      delete elInfo;

Praetorius, Simon's avatar
Praetorius, Simon committed
159
160
    return ((*val));
  };
161
162


163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
  template<>
  double DOFVector<double>::IntOnBoundary(BoundaryType boundaryType, 
                                          Quadrature* q) const
  {
    FUNCNAME("DOFVector<T>::IntOnBoundary()");
  
    Mesh* mesh = this->feSpace->getMesh();
    int dim = mesh->getDim();
  
    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(dim - 1, deg);
    } else {
      TEST_EXIT(q->getDim() == dim-1)
        ("Provided quadrature has wrong dimension for surfaceQuadrature!\n");
    }
  
    // create barycentric coords for each vertex of each side
    const Element *refElement = Global::getReferenceElement(dim);
    VectorOfFixVecs<DimVec<double> >**coords;
    coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
  
    // for all element sides
    for (int i = 0; i < dim + 1; i++) {
      coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
        DimVec<double>(dim, DEFAULT_VALUE, 0.0));
      // for each vertex of the side
      for (int k = 0; k < dim; k++) {
        int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), 
                                                    i, k);
        (*coords[i])[k][index] = 1.0;
      }
    }
  
    std::vector<SurfaceQuadrature*> quadSurfaces(dim + 1);
    for (int i = 0; i < dim + 1; i++)
      quadSurfaces[i] = new SurfaceQuadrature(q, *coords[i]);
  
    double result = 0.0;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
      Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
    while (elInfo) {
      for (int face = 0; face < dim + 1; face++) {
        if (elInfo->getBoundary(face) == boundaryType) {
          int nPoints = quadSurfaces[face]->getNumPoints();
          mtl::dense_vector<double> uh_vec(nPoints);
          double det = elInfo->calcSurfaceDet(*coords[face]);
          double normT = 0.0;
          this->getVecAtQPs(elInfo, quadSurfaces[face], NULL, uh_vec);
          for (int iq = 0; iq < nPoints; iq++)
            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
          result += det * normT;
        }
      }
  
      elInfo = stack.traverseNext(elInfo);
    }
  
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
    #endif
    
    return result;
  }


  template<>
  double DOFVector<WorldVector<double> >::IntOnBoundaryNormal(
    BoundaryType boundaryType, Quadrature* q) const
  {
    FUNCNAME("DOFVector<T>::IntOnBoundary()");
 
    Mesh* mesh = this->feSpace->getMesh();
    int dim = mesh->getDim();
  
    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(dim - 1, deg);
    } else {
      TEST_EXIT(q->getDim() == dim-1)
        ("Provided quadrature has wrong dimension for surfaceQuadrature!\n");
    }
  
    // create barycentric coords for each vertex of each side
    const Element *refElement = Global::getReferenceElement(dim);
    VectorOfFixVecs<DimVec<double> >**coords;
    coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
  
    // for all element sides
    for (int i = 0; i < dim + 1; i++) {
      coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
        DimVec<double>(dim, DEFAULT_VALUE, 0.0));
      // for each vertex of the side
      for (int k = 0; k < dim; k++) {
        int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), 
                                                    i, k);
        (*coords[i])[k][index] = 1.0;
      }
    }
  
    std::vector<SurfaceQuadrature*> quadSurfaces(dim + 1);
    for (int i = 0; i < dim + 1; i++)
      quadSurfaces[i] = new SurfaceQuadrature(q, *coords[i]);
  
    double result = 0.0;
    WorldVector<double> normal;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
      Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
    while (elInfo) {
      for (int face = 0; face < dim + 1; face++) {
        if (elInfo->getBoundary(face) == boundaryType) {
          elInfo->getNormal(face, normal);
          double det = elInfo->calcSurfaceDet(*coords[face]);
  
          int nPoints = quadSurfaces[face]->getNumPoints();
          mtl::dense_vector<WorldVector<double> > uh_vec(nPoints);
          WorldVector<double> normT; normT.set(0.0);
          this->getVecAtQPs(elInfo, quadSurfaces[face], NULL, uh_vec);
          for (int iq = 0; iq < nPoints; iq++)
            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
          // scalar product between vector-valued solution and normal vector
          result += det * (normT * normal); 
        }
      }
  
      elInfo = stack.traverseNext(elInfo);
    }
  
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
    #endif
    
    return result;
  }


303
  template<>
304
305
306
307
  void DOFVectorBase<double>::getD2AtQPs( const ElInfo *elInfo,
					  const Quadrature *quad,
					  const FastQuadrature *quadFast,
					  mtl::dense_vector<D2Type<double>::type> &d2AtQPs) const
308
309
310
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
311
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
312

313
314
315
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
316
317
    }

318
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
319
320
      ("invalid basis functions");

321
322
    Element *el = elInfo->getElement(); 

323
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
324
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
325
  
326
    mtl::dense_vector<double> localVec(nBasFcts);
327
    getLocalVector(el, localVec);
328
329
330
331
332

    DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0);
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

333
    d2AtQPs.change_dim(nPoints);
334
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
337
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
338
339
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
340
341
342
343
	for (int i = 0; i < nBasFcts; i++) {
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
	      D2Tmp[k][l] += localVec[i] * quadFast->getSecDer(iq, i, k, l);
344
345
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
346
347
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
348
	    d2AtQPs[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
351
		d2AtQPs[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
352
353
354
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
355
      const BasisFunction *basFcts = feSpace->getBasisFcts();
356
      DimMat<double> D2Phi(dim, NO_INIT);
357

Thomas Witkowski's avatar
Thomas Witkowski committed
358
359
360
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
361
362
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
363
	for (int i = 0; i < nBasFcts; i++) {
364
	  WARNING("not tested after index correction\n");
365
	  (*(basFcts->getD2Phi(i)))(quad->getLambda(iq), D2Phi);
366

Thomas Witkowski's avatar
Thomas Witkowski committed
367
368
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
369
	      D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
370
371
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
374
	    d2AtQPs[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
375
376
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
377
		d2AtQPs[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
378
379
380
381
382
	  }
      }
    }
  }

383

384
385
386
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
387
    FUNCNAME("DOFVector<double>::interpol()");
388

389
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
390
391
392
393
394
395
396
397
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

    int nBasisFcts = basisFcts->getNumber();
    int nSourceBasisFcts = sourceBasisFcts->getNumber();

    this->set(0.0);

398
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
399
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
400

401
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
402
403
404
405
406
407
      DimVec<double> *coords = NULL;
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

408
      while (elInfo) {
409
410
	Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
411
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
412
413
	source->getLocalVector(el, sourceLocalCoeffs);

414
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
415
	  if (vec[myLocalIndices[i]] == 0.0) {
416
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
417
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      DimVec<double> coords2(feSpace->getMesh()->getDim(), NO_INIT);
      DualTraverse dualStack;
      ElInfo *elInfo1, *elInfo2;
      ElInfo *elInfoSmall, *elInfoLarge;
      WorldVector<double> worldVec;

      bool nextTraverse = dualStack.traverseFirst(feSpace->getMesh(),
						  sourceFeSpace->getMesh(),
						  -1, -1,
						  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS,
						  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS,
						  &elInfo1, &elInfo2,
						  &elInfoSmall, &elInfoLarge);
436
      while (nextTraverse) {     
437
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
438
				   myLocalIndices);
439
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
440

441
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
442
	  if (vec[myLocalIndices[i]] == 0.0) {
443
            elInfo1->coordToWorld(*(basisFcts->getCoords(i)), worldVec);
444
445
446
	    elInfo2->worldToCoord(worldVec, &coords2);
	  
	    bool isPositive = true;
447
	    for (int j = 0; j < coords2.size(); j++) {
448
	      if (coords2[j] < -0.00001) {
449
450
451
452
453
		isPositive = false;
		break;
	      }
	    }
	  
454
455
456
	    if (isPositive)
	      vec[myLocalIndices[i]] = 
		sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	   
457
458
459
	  }
	}
      
460
461
	nextTraverse = 
	  dualStack.traverseNext(&elInfo1, &elInfo2, &elInfoSmall, &elInfoLarge);
462
463
464
465
      }
    }
  }

466
467

  template<>
468
469
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
470
471
472
473
474
475
  {
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

    DimVec<double> *coords = NULL;
476
    const FiniteElemSpace *vFeSpace = v->getFeSpace();
477

478
    if (feSpace == vFeSpace)
479
480
481
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
482
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
483

484
    int nBasFcts = basFcts->getNumber();
485
486
    int vNumBasFcts = vBasFcts->getNumber();

487
488
    if (feSpace->getMesh() == vFeSpace->getMesh()) {      
      std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
489
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
490
491
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
492
493
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
494
495
496

      while (elInfo) {
	Element *el = elInfo->getElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
497
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
498
499
	v->getLocalVector(el, vLocalCoeffs);

500
	for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
501
	  if (vec[myLocalIndices[i]] == nul) {
502
	    coords = basFcts->getCoords(i);
503
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs) * factor;
504
505
506
507
508
509
510
511
512
513
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
    FUNCNAME("DOFVector<double>::getGradient()");

    Mesh *mesh = feSpace->getMesh();
    int dim = mesh->getDim();
    int dow = Global::getGeo(WORLD);

    const BasisFunction *basFcts = feSpace->getBasisFcts();

    DOFAdmin *admin = feSpace->getAdmin();

    // define result vector
    static WorldVector<DOFVector<double>*> *result = NULL;

530
    if (grad) {
531
532
      result = grad;
    } else {
533
      if (!result) {
534
	result = new WorldVector<DOFVector<double>*>;
535

536
537
	result->set(NULL);
      }
538
      for (int i = 0; i < dow; i++) {
539
	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {
540
541
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
542
543
544
545
546
	}
      }
    }

    // count number of nodes and dofs per node
547
548
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
549
    std::vector<DimVec<double>*> bary;
550

551
552
    int nNodes = 0;
    int nDofs = 0;
553

554
    for (int i = 0; i < dim + 1; i++) {
555
556
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
557
      int numPreDofs = admin->getNumberOfPreDofs(i);
558
      for (int j = 0; j < numPositionNodes; j++) {
559
	int dofs = basFcts->getNumberOfDofs(geoIndex);
560
	nNodeDOFs.push_back(dofs);
561
	nDofs += dofs;
562
	nNodePreDofs.push_back(numPreDofs);
563
      }
564
      nNodes += numPositionNodes;
565
566
    }

567
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
568
      ("number of dofs != number of basis functions\n");
569
    
570
    for (int i = 0; i < nDofs; i++)
571
      bary.push_back(basFcts->getCoords(i));    
572
573

    // traverse mesh
574
    std::vector<bool> visited(getUsedSize(), false);
575
    TraverseStack stack;
576
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
577
578
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
579
    ElementVector localUh(basFcts->getNumber());
580

581
    while (elInfo) {
582
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
583
      getLocalVector(elInfo->getElement(), localUh);
584
585
586
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
587
      for (int i = 0; i < nNodes; i++) { // for all nodes
588
589
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
590

591
	  if (!visited[dofIndex]) {
592
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, grd);
593

594
595
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
596
597
598
599
600
601
602
603
604
605
606
607
608
609

	    visited[dofIndex] = true;
	  }
	  localDOFNr++;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


610
611
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
612
  {}
613

614

615
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
616
617
  {}

618
619
620
621

  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *res)
  {
622
    FUNCNAME("DOFVector<double>::transform()");
623

624
    TEST_EXIT_DBG(vec)("no vector\n");
625

626
    int dow = Global::getGeo(WORLD);
627
628
    static WorldVector<DOFVector<double>*> *result = NULL;

629
    if (!res && !result) {
630
631
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
632
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
633
634
635
636
637
    }

    WorldVector<DOFVector<double>*> *r = res ? res : result;

    int vecSize = vec->getSize();
638
639
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
640
641
642
643
644
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

645
646


Thomas Witkowski's avatar
Thomas Witkowski committed
647
648
649
650
651
652
653
654
655
656
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    FUNCNAME("DOFVector::assemble()");

    if (!(op || this->operators.size())) 
      return;

657
    set_to_zero(this->elementVector);
658
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
659
660
661

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
662
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
663
664
665
666
667
668
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	   it != this->operators.end(); 
669
	   ++it, ++factorIt)
670
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
671
	  (*it)->getElementVector(elInfo, this->elementVector, 
672
				  *factorIt ? **factorIt : 1.0);      
673
674
	  addVector = true;
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
675
676
    }

677
678
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
681

Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
684
685
686
687
688
689
690
691
692
693
  template<>
  void DOFVectorBase<double>::assemble2(double factor, 
					ElInfo *mainElInfo, ElInfo *auxElInfo,
					ElInfo *smallElInfo, ElInfo *largeElInfo,
					const BoundaryType *bound, 
					Operator *op)
  {
    FUNCNAME("DOFVector::assemble2()");

    if (!(op || this->operators.size())) 
      return;

694
    set_to_zero(this->elementVector);
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
    bool addVector = false;

    TEST_EXIT(!op)("Not yet implemented!\n");

    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator factorIt;
    for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	 it != this->operators.end(); 
	 ++it, ++factorIt)
      if ((*it)->getNeedDualTraverse()) {
	(*it)->getElementVector(mainElInfo, auxElInfo,
				smallElInfo, largeElInfo,
				this->elementVector, 
				*factorIt ? **factorIt : 1.0);
	addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
710
711
      }
  
712
713
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
716
  }


717
718
719
720
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct)
  {
721
    if (vec1.getFeSpace()->getMesh() == vec2.getFeSpace()->getMesh())
722
723
724
725
726
727
728
729
730
731
732
733
734
735
      return intSingleMesh(vec1, vec2, fct);
    else
      return intDualMesh(vec1, vec2, fct);
  }


  double intSingleMesh(const DOFVector<double> &vec1,
		       const DOFVector<double> &vec2,
		       BinaryAbstractFunction<double, double, double> *fct)
  {
    FUNCNAME("intDualmesh()");
    
    TEST_EXIT(fct)("No function defined!\n");

736
737
    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
738
    Quadrature* quad = 
739
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
740
    FastQuadrature *fastQuad =
741
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
742

743
744
    mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
745
746
747
748

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
749
    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
750
    while (elInfo) {
751
752
753
754
      vec1.getVecAtQPs(elInfo, quad, fastQuad, qp1);
      vec2.getVecAtQPs(elInfo, quad, fastQuad, qp2);

      double tmp = 0.0;
755
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
756
757
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);     
      value += tmp * elInfo->getDet();
758
759
760
761
      
      elInfo = stack.traverseNext(elInfo);
    }

762
763
764
765
766
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
767
768
769
770
771
772
773
774
775
776
777
778
    return value;
  }


  double intDualMesh(const DOFVector<double> &vec1,
		     const DOFVector<double> &vec2,
		     BinaryAbstractFunction<double, double, double> *fct)
  {
    FUNCNAME("intDualmesh()");
    
    TEST_EXIT(fct)("No function defined!\n");

779
780
    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
781
    Quadrature* quad = 
782
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
783
    FastQuadrature *fastQuad =
784
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
785

786
787
    mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
788
789
790
791
792

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    DualTraverse dualTraverse;    
    DualElInfo dualElInfo;
793
794
    bool cont = dualTraverse.traverseFirst(vec1.getFeSpace()->getMesh(),
					   vec2.getFeSpace()->getMesh(),
795
796
797
					   -1, -1, traverseFlag, traverseFlag,
					   dualElInfo);
    while (cont) {
798
799
800
801
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);

      double tmp = 0.0;
802
      for (int iq = 0; iq < quad->getNumPoints(); iq++)
803
804
 	tmp += quad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);      
      value += tmp * dualElInfo.smallElInfo->getDet();
805
806
807
808
      
      cont = dualTraverse.traverseNext(dualElInfo);
    }

809
810
811
812
813
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif

814
815
816
    return value;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
817
818
819
820
821
822
823
824
  
  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate()");
    
    TEST_EXIT(fct)("No function defined!\n");

825
826
827
828
829
    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);
Thomas Witkowski's avatar
Thomas Witkowski committed
830
    FastQuadrature *fastQuad =
831
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);
Thomas Witkowski's avatar
Thomas Witkowski committed
832
833
834
835
836
837
838

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

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
839
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
Thomas Witkowski's avatar
Thomas Witkowski committed
840
841
842
843
844
845
846
847
848
849
850
851
852
853
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);

      double tmp = 0.0;
      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);
    }

854
855
856
857
858
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
859
860
    return value;
  }
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893


  double integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<double, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate()");
    
    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;

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
    ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, traverseFlag);
    while (elInfo) {
      double tmp = 0.0;
      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);
    }

894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return value;
  }



  double integrateGeneral(const std::vector<DOFVector<double>*> &vecs,
		       AbstractFunction<double, std::vector<double> > *fct)
  {
    FUNCNAME("integrateGeneral()");

    TEST_EXIT(fct)("No function defined!\n");
    TEST_EXIT(vecs.size()>0)("No DOFVectors provided!\n");

    int deg = std::max(fct->getDegree(),
Praetorius, Simon's avatar
Praetorius, Simon committed
913
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
914
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
915
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
916
    FastQuadrature *fastQuad =
Praetorius, Simon's avatar
Praetorius, Simon committed
917
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
918
919
920
921
922
923
924
925
926

    std::vector<mtl::dense_vector<double> > qp(vecs.size());
    std::vector<double> qp_local(vecs.size());
    for (size_t i = 0; i < vecs.size(); i++)
      qp[i].change_dim(fastQuad->getNumPoints());

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
Praetorius, Simon's avatar
Praetorius, Simon committed
927
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
928
929
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
930
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952

      double tmp = 0.0;
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	for (size_t i = 0; i < vecs.size(); i++)
	  qp_local[i] = qp[i][iq];
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp_local);
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return value;
  }

  double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
				  const std::vector<DOFVector<double>*> &grds,
Praetorius, Simon's avatar
Praetorius, Simon committed
953
				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > > *fct)
954
955
956
957
958
959
960
961
  {
    FUNCNAME("integrateGeneral()");

    TEST_EXIT(fct)("No function defined!\n");
    TEST_EXIT(vecs.size()>0)("No DOFVectors provided!\n");
    TEST_EXIT(grds.size()>0)("No DOFVectors for gradients provided!\n");

    int deg = std::max(fct->getDegree(),
Praetorius, Simon's avatar
Praetorius, Simon committed
962
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
963
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
964
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
965
    FastQuadrature *fastQuad =
Praetorius, Simon's avatar
Praetorius, Simon committed
966
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
967
968
969
970
971
972
973
974
975
976
977
978
979

    std::vector<mtl::dense_vector<double> > qp(vecs.size());
    std::vector<mtl::dense_vector<WorldVector<double> > > qpGrd(vecs.size());
    std::vector<double> qp_local(vecs.size());
    std::vector<WorldVector<double> > grd_local(grds.size());
    for (size_t i = 0; i < vecs.size(); i++)
      qp[i].change_dim(fastQuad->getNumPoints());
    for (size_t i = 0; i < grds.size(); i++)
      qpGrd[i].change_dim(fastQuad->getNumPoints());

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
Praetorius, Simon's avatar
Praetorius, Simon committed
980
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
981
982
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
983
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
984
      for (size_t i = 0; i < grds.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
985
	grds[i]->getGrdAtQPs(elInfo, quad, fastQuad, qpGrd[i]);
986
987
988
989
990
991
992
993
994
995
996
997
998
999

      double tmp = 0.0;
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	for (size_t i = 0; i < vecs.size(); i++)
	  qp_local[i] = qp[i][iq];
	for (size_t i = 0; i < grds.size(); i++)
	  grd_local[i] = qpGrd[i][iq];
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp_local, grd_local);
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

1000
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS