DOFVector.cc 24.9 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<>
71
72
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p,
					ElInfo *oldElInfo) const
73
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
74
    FUNCNAME("DOFVector<double>::evalAtCoords()");
75
76
77
78
79
80
81
  
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Praetorius, Simon's avatar
Praetorius, Simon committed
82
    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
83
84
85
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();
86
    double value = 0.0;
87
88
89
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
Praetorius, Simon's avatar
Praetorius, Simon committed
90
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
91
92
93
94
95
96
97
98
      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
99
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
100
101
      ElementVector uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
102
103
104
105
106
        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
107
    delete [] localIndices;
108
109
110
111
    if (oldElInfo == NULL)
      delete elInfo;

    return value;
112
  }
113
114
115


  template<>
116
117
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) 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();

132
    WorldVector<double> value(DEFAULT_VALUE, 0.0);
133
134
135
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
Praetorius, Simon's avatar
Praetorius, Simon committed
136
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
137
138
139
140
141
142
143
144
      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
145
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
146
147
      mtl::dense_vector<WorldVector<double> > uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
148
        uh[i] = operator[](localIndices[i]);
149
      value = basFcts->evalUh(lambda, uh);
150
151
152
    } else
      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));

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

157
158
    return value;
  }
159
160


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
  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);
    }
  
220
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
221
222
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
223
#endif
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
    
    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;
  }


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

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

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

319
320
    Element *el = elInfo->getElement(); 

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
340
341
	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);
342
343
	}

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

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

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

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

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

381

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

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

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

    this->set(0.0);

396
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
397
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
398

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

406
      while (elInfo) {
407
408
	Element *el = elInfo->getElement();

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

412
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
413
	  if (vec[myLocalIndices[i]] == 0.0) {
414
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
415
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
	  }
	}
	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);
434
      while (nextTraverse) {     
435
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
436
				   myLocalIndices);
437
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
438

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

464
465

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

    this->set(nul);

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

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

    const BasisFunction *basFcts = feSpace->getBasisFcts();
480
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
481

482
    int nBasFcts = basFcts->getNumber();
483
484
    int vNumBasFcts = vBasFcts->getNumber();

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

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

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


512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
  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;

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

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

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

549
550
    int nNodes = 0;
    int nDofs = 0;
551

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

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

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


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

612

613
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
614
615
  {}

616
617
618
619

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

622
    TEST_EXIT_DBG(vec)("no vector\n");
623

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

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

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

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

    return r;
  }

643
644


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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
679

Thomas Witkowski's avatar
Thomas Witkowski committed
680
681
682
683
684
685
686
687
688
689
690
691
  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;

692
    set_to_zero(this->elementVector);
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
    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
708
709
      }
  
710
711
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
712
713
714
  }


715
716
717
718
719
720
721
722
723
  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
724
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
725
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
726
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
727
    FastQuadrature *fastQuad =
Praetorius, Simon's avatar
Praetorius, Simon committed
728
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
729
730
731
732
733
734
735
736
737

    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
738
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
739
740
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
741
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763

      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
764
				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > > *fct)
765
766
767
768
769
770
771
772
  {
    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
773
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
774
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
775
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
776
    FastQuadrature *fastQuad =
Praetorius, Simon's avatar
Praetorius, Simon committed
777
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
778
779
780
781
782
783
784
785
786
787
788
789
790

    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
791
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
792
793
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
794
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
795
      for (size_t i = 0; i < grds.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
796
	grds[i]->getGrdAtQPs(elInfo, quad, fastQuad, qpGrd[i]);
797
798
799
800
801
802
803
804
805
806
807
808
809
810

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

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

816
817
    return value;
  }
818
819
}