DOFVector.cc 25.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
#include <boost/numeric/mtl/mtl.hpp>
23
24
25
26
27
28
29
30
31
32
#include "DOFVector.h"
#include "Traverse.h"
#include "DualTraverse.h"
#include "FixVec.h"

namespace AMDiS {

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

35
    switch (coarsenOperation) {
36
37
38
39
40
41
42
    case NO_OPERATION:
      return;
      break;
    case COARSE_RESTRICT:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
      break;
    case COARSE_INTERPOL:
43
44
      TEST_EXIT_DBG(feSpace)("Should not happen!\n");
      TEST_EXIT_DBG(feSpace->getBasisFcts())("Shoud not happen!\n");
45

46
47
48
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
      break;
    default:
49
50
      WARNING("Invalid coarsen operation \"%d\" in vector \"%s\"\n", 
	      coarsenOperation, this->name.c_str());
51
52
53
    }
  }

54

55
56
57
58
59
60
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList& list, int n)
  {
    (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->refineInter(this, &list, n);
  }

61

62
63
64
  template<>
  void DOFVector<WorldVector<double> >::refineInterpol(RCNeighbourList& list, int n)
  {
65
66
67
    if (n < 1) 
      return;

68
    Element *el = list.getElement(0);
69
    int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
70
71
72
    DegreeOfFreedom dof0 = el->getDof(0, n0);
    DegreeOfFreedom dof1 = el->getDof(1, n0);
    DegreeOfFreedom dof_new = el->getChild(0)->getDof(feSpace->getMesh()->getDim(), n0);
73
74
75
76
77
    vec[dof_new] = vec[dof0];
    vec[dof_new] += vec[dof1];
    vec[dof_new] *= 0.5;
  }

78

79
  template<>
80
81
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p,
					ElInfo *oldElInfo) const
82
  {  
83
84
85
86
87
88
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
89
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
90
91
92
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();
93
    double value = 0.0;
94
95
96
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
97
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), nullptr, nullptr);
98
99
      delete oldElInfo;
    } else
100
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
101
102
103
104
105

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Praetorius, Simon's avatar
Praetorius, Simon committed
106
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
107
108
      ElementVector uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
109
110
        uh[i] = operator[](localIndices[i]);
      value = basFcts->evalUh(lambda, uh);
111
112
113
114
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      return 0.0;
#else
115
      ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
116
117
118
#endif
    }

119

120
    if (oldElInfo == nullptr)
121
122
123
      delete elInfo;

    return value;
124
  }
125
126
127


  template<>
128
129
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const
130
  {  
131
132
133
134
135
136
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
137
138
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
    DimVec<double> lambda(dim, NO_INIT);  
139
    ElInfo *elInfo = mesh->createNewElInfo();
140
    WorldVector<double> value(DEFAULT_VALUE, 0.0);
141
142
143
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
144
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), nullptr, nullptr);
145
146
      delete oldElInfo;
    } else
147
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
148
149
150
151
152

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Praetorius, Simon's avatar
Praetorius, Simon committed
153
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
154
155
      mtl::dense_vector<WorldVector<double> > uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
156
        uh[i] = operator[](localIndices[i]);
157
      value = basFcts->evalUh(lambda, uh);
158
159
160
    } else {
      ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
    }
161

162
    if (oldElInfo == nullptr)
163
164
      delete elInfo;

165
166
    return value;
  }
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
  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;
218
          this->getVecAtQPs(elInfo, quadSurfaces[face], nullptr, uh_vec);
219
220
221
222
223
224
225
226
227
          for (int iq = 0; iq < nPoints; iq++)
            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
          result += det * normT;
        }
      }
  
      elInfo = stack.traverseNext(elInfo);
    }
  
228
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
229
230
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
231
#endif
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
    
    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);
289
          this->getVecAtQPs(elInfo, quadSurfaces[face], nullptr, uh_vec);
290
291
292
293
294
295
296
297
298
299
          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);
    }
  
300
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
301
302
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
303
#endif
304
305
306
307
308
    
    return result;
  }


309
  template<>
310
311
312
313
  void DOFVectorBase<double>::getD2AtQPs( const ElInfo *elInfo,
					  const Quadrature *quad,
					  const FastQuadrature *quadFast,
					  mtl::dense_vector<D2Type<double>::type> &d2AtQPs) const
314
315
316
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
317
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
318

319
320
321
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
322
323
    }

324
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
325
326
      ("invalid basis functions");

327
328
    Element *el = elInfo->getElement(); 

329
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
330
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
331
  
332
    mtl::dense_vector<double> localVec(nBasFcts);
333
    getLocalVector(el, localVec);
334
335
336
337
338

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

339
    d2AtQPs.change_dim(nPoints);
340
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
341
342
343
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
344
345
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
346
347
348
349
	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);
350
351
	}

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

Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
366
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
367
368
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
375
	      D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
376
377
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
378
379
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
380
	    d2AtQPs[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
381
382
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
383
		d2AtQPs[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
384
385
386
387
388
	  }
      }
    }
  }

389

390
391
392
393
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {

394
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
395
396
397
398
399
400
401
402
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

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

    this->set(0.0);

403
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
404
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
405

406
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
407
      DimVec<double> *coords = nullptr;
408
409
410
411
412
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

413
      while (elInfo) {
414
415
	Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
416
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
417
418
	source->getLocalVector(el, sourceLocalCoeffs);

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

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

471
472

  template<>
473
474
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
475
  {
476
    FUNCNAME("DOFVector<WorldVector<double> >::interpol()");
477
478
479
480
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

481
    DimVec<double> *coords = nullptr;
482
    const FiniteElemSpace *vFeSpace = v->getFeSpace();
483

484
    if (feSpace == vFeSpace)
485
486
487
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
488
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
489

490
    int nBasFcts = basFcts->getNumber();
491
492
    int vNumBasFcts = vBasFcts->getNumber();

493
494
    if (feSpace->getMesh() == vFeSpace->getMesh()) {      
      std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
495
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
496
497
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
498
499
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
500
501
502

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

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


520
521
522
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
523
    FUNCNAME_DBG("DOFVector<double>::getGradient()");
524
525
526
527
528
529
530
531
532
533

    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
534
    static WorldVector<DOFVector<double>*> *result = nullptr; // TODO: REMOVE STATIC
535

536
    if (grad) {
537
538
      result = grad;
    } else {
539
      if (!result) {
540
	result = new WorldVector<DOFVector<double>*>;
541

542
	result->set(nullptr);
543
      }
544
      for (int i = 0; i < dow; i++) {
545
	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {
546
547
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
548
549
550
551
552
	}
      }
    }

    // count number of nodes and dofs per node
553
554
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
555
    std::vector<DimVec<double>*> bary;
556

557
558
    int nNodes = 0;
    int nDofs = 0;
559

560
    for (int i = 0; i < dim + 1; i++) {
561
562
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
563
      int numPreDofs = admin->getNumberOfPreDofs(i);
564
      for (int j = 0; j < numPositionNodes; j++) {
565
	int dofs = basFcts->getNumberOfDofs(geoIndex);
566
	nNodeDOFs.push_back(dofs);
567
	nDofs += dofs;
568
	nNodePreDofs.push_back(numPreDofs);
569
      }
570
      nNodes += numPositionNodes;
571
572
    }

573
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
574
      ("number of dofs != number of basis functions\n");
575
    
576
    for (int i = 0; i < nDofs; i++)
577
      bary.push_back(basFcts->getCoords(i));    
578
579

    // traverse mesh
580
    std::vector<bool> visited(getUsedSize(), false);
581
    TraverseStack stack;
582
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
583
584
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
585
    ElementVector localUh(basFcts->getNumber());
586

587
    while (elInfo) {
588
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
589
      getLocalVector(elInfo->getElement(), localUh);
590
591
592
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
593
      for (int i = 0; i < nNodes; i++) { // for all nodes
594
595
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
596

597
	  if (!visited[dofIndex]) {
598
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, grd);
599

600
601
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
602
603
604
605
606
607
608
609
610
611
612
613
614
615

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


616
617
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
618
  {}
619

620

621
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
622
623
  {}

624
625
626
627

  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *res)
  {
628
    FUNCNAME_DBG("DOFVector<double>::transform()");
629

630
    TEST_EXIT_DBG(vec)("no vector\n");
631

632
    int dow = Global::getGeo(WORLD);
633
    static WorldVector<DOFVector<double>*> *result = nullptr; // TODO: REMOVE STATIC
634

635
    if (!res && !result) {
636
637
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
638
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
639
640
641
642
643
    }

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

    int vecSize = vec->getSize();
644
645
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
646
647
648
649
650
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

651
652


Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
655
656
657
658
659
660
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    if (!(op || this->operators.size())) 
      return;

661
    set_to_zero(this->elementVector);
662
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
663
664
665

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

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

681
682
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
683
684
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
685

Thomas Witkowski's avatar
Thomas Witkowski committed
686
687
688
689
690
691
692
693
694
695
696
697
  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;

698
    set_to_zero(this->elementVector);
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
    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
714
715
      }
  
716
717
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
720
  }


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

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

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

    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;
795
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;
796
    TraverseStack stack;
Praetorius, Simon's avatar
Praetorius, Simon committed
797
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
798
799
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
800
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
801
      for (size_t i = 0; i < grds.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
802
	grds[i]->getGrdAtQPs(elInfo, quad, fastQuad, qpGrd[i]);
803
804
805
806
807
808
809
810
811
812
813
814
815
816

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

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

822
823
    return value;
  }
824
825
}