DOFVector.cc 25.5 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
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList& list, int n)
  {
58
59
60
61
62
63
64
65
66
    switch (refineOperation) {
    case NO_OPERATION:
      return;
      break;
    case REFINE_INTERPOL:
    default:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->refineInter(this, &list, n);
      break;
    }
67
68
  }

69

70
71
72
  template<>
  void DOFVector<WorldVector<double> >::refineInterpol(RCNeighbourList& list, int n)
  {
73
74
75
    if (refineOperation == NO_OPERATION)
      return;
    
76
77
78
    if (n < 1) 
      return;

79
    Element *el = list.getElement(0);
80
    int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
81
82
83
    DegreeOfFreedom dof0 = el->getDof(0, n0);
    DegreeOfFreedom dof1 = el->getDof(1, n0);
    DegreeOfFreedom dof_new = el->getChild(0)->getDof(feSpace->getMesh()->getDim(), n0);
84
85
86
87
88
    vec[dof_new] = vec[dof0];
    vec[dof_new] += vec[dof1];
    vec[dof_new] *= 0.5;
  }

89

90
  template<>
91
92
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p,
					ElInfo *oldElInfo) const
93
  {  
94
95
96
97
98
99
    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
100
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
101
102
103
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();
104
    double value = 0.0;
105
106
107
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
108
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), nullptr, nullptr);
109
110
      delete oldElInfo;
    } else
111
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
112
113
114
115
116

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Praetorius, Simon's avatar
Praetorius, Simon committed
117
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
118
119
      ElementVector uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
120
121
        uh[i] = operator[](localIndices[i]);
      value = basFcts->evalUh(lambda, uh);
122
123
124
125
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      return 0.0;
#else
126
      ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
127
128
129
#endif
    }

130

131
    if (oldElInfo == nullptr)
132
133
134
      delete elInfo;

    return value;
135
  }
136
137
138


  template<>
139
140
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const
141
  {  
142
143
144
145
146
147
    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
148
149
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
    DimVec<double> lambda(dim, NO_INIT);  
150
    ElInfo *elInfo = mesh->createNewElInfo();
151
    WorldVector<double> value(DEFAULT_VALUE, 0.0);
152
153
154
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
155
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), nullptr, nullptr);
156
157
      delete oldElInfo;
    } else
158
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
159
160
161
162
163

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Praetorius, Simon's avatar
Praetorius, Simon committed
164
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
165
166
      mtl::dense_vector<WorldVector<double> > uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
167
        uh[i] = operator[](localIndices[i]);
168
      value = basFcts->evalUh(lambda, uh);
169
170
171
    } else {
      ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
    }
172

173
    if (oldElInfo == nullptr)
174
175
      delete elInfo;

176
177
    return value;
  }
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
  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;
229
          this->getVecAtQPs(elInfo, quadSurfaces[face], nullptr, uh_vec);
230
231
232
233
234
235
236
237
238
          for (int iq = 0; iq < nPoints; iq++)
            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
          result += det * normT;
        }
      }
  
      elInfo = stack.traverseNext(elInfo);
    }
  
239
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
240
241
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
242
#endif
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
    
    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);
300
          this->getVecAtQPs(elInfo, quadSurfaces[face], nullptr, uh_vec);
301
302
303
304
305
306
307
308
309
310
          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);
    }
  
311
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
312
313
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
314
#endif
315
316
317
318
319
    
    return result;
  }


320
  template<>
321
322
323
324
  void DOFVectorBase<double>::getD2AtQPs( const ElInfo *elInfo,
					  const Quadrature *quad,
					  const FastQuadrature *quadFast,
					  mtl::dense_vector<D2Type<double>::type> &d2AtQPs) const
325
326
327
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
328
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
329

330
331
332
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
333
334
    }

335
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
336
337
      ("invalid basis functions");

338
339
    Element *el = elInfo->getElement(); 

340
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
341
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
342
  
343
    mtl::dense_vector<double> localVec(nBasFcts);
344
    getLocalVector(el, localVec);
345
346
347
348
349

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

350
    d2AtQPs.change_dim(nPoints);
351
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
354
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
355
356
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
357
358
359
360
	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);
361
362
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
363
364
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
365
	    d2AtQPs[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
366
367
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
368
		d2AtQPs[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
369
370
371
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
372
      const BasisFunction *basFcts = feSpace->getBasisFcts();
373
      DimMat<double> D2Phi(dim, NO_INIT);
374

Thomas Witkowski's avatar
Thomas Witkowski committed
375
376
377
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
378
379
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
384
385
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
386
	      D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
387
388
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
391
	    d2AtQPs[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
394
		d2AtQPs[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
395
396
397
398
399
	  }
      }
    }
  }

400

401
402
403
404
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {

405
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
406
407
408
409
410
411
412
413
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

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

    this->set(0.0);

414
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
415
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
416

417
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
418
      DimVec<double> *coords = nullptr;
419
420
421
422
423
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

424
      while (elInfo) {
425
426
	Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
427
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
428
429
	source->getLocalVector(el, sourceLocalCoeffs);

430
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
431
	  if (vec[myLocalIndices[i]] == 0.0) {
432
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
433
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
	  }
	}
	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);
452
      while (nextTraverse) {     
453
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
454
				   myLocalIndices);
455
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
456

457
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
458
	  if (vec[myLocalIndices[i]] == 0.0) {
459
            elInfo1->coordToWorld(*(basisFcts->getCoords(i)), worldVec);
460
461
462
	    elInfo2->worldToCoord(worldVec, &coords2);
	  
	    bool isPositive = true;
463
	    for (int j = 0; j < coords2.size(); j++) {
464
	      if (coords2[j] < -0.00001) {
465
466
467
468
469
		isPositive = false;
		break;
	      }
	    }
	  
470
471
472
	    if (isPositive)
	      vec[myLocalIndices[i]] = 
		sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	   
473
474
475
	  }
	}
      
476
477
	nextTraverse = 
	  dualStack.traverseNext(&elInfo1, &elInfo2, &elInfoSmall, &elInfoLarge);
478
479
480
481
      }
    }
  }

482
483

  template<>
484
485
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
486
  {
487
    FUNCNAME("DOFVector<WorldVector<double> >::interpol()");
488
489
490
491
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

492
    DimVec<double> *coords = nullptr;
493
    const FiniteElemSpace *vFeSpace = v->getFeSpace();
494

495
    if (feSpace == vFeSpace)
496
497
498
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
499
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
500

501
    int nBasFcts = basFcts->getNumber();
502
503
    int vNumBasFcts = vBasFcts->getNumber();

504
505
    if (feSpace->getMesh() == vFeSpace->getMesh()) {      
      std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
506
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
507
508
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
509
510
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
511
512
513

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

517
	for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
518
	  if (vec[myLocalIndices[i]] == nul) {
519
	    coords = basFcts->getCoords(i);
520
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs) * factor;
521
522
523
524
525
526
527
528
529
530
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


531
532
533
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
534
    FUNCNAME_DBG("DOFVector<double>::getGradient()");
535
536
537
538
539
540
541
542
543
544

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

547
    if (grad) {
548
549
      result = grad;
    } else {
550
      if (!result) {
551
	result = new WorldVector<DOFVector<double>*>;
552

553
	result->set(nullptr);
554
      }
555
      for (int i = 0; i < dow; i++) {
556
	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {
557
558
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
559
560
561
562
563
	}
      }
    }

    // count number of nodes and dofs per node
564
565
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
566
    std::vector<DimVec<double>*> bary;
567

568
569
    int nNodes = 0;
    int nDofs = 0;
570

571
    for (int i = 0; i < dim + 1; i++) {
572
573
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
574
      int numPreDofs = admin->getNumberOfPreDofs(i);
575
      for (int j = 0; j < numPositionNodes; j++) {
576
	int dofs = basFcts->getNumberOfDofs(geoIndex);
577
	nNodeDOFs.push_back(dofs);
578
	nDofs += dofs;
579
	nNodePreDofs.push_back(numPreDofs);
580
      }
581
      nNodes += numPositionNodes;
582
583
    }

584
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
585
      ("number of dofs != number of basis functions\n");
586
    
587
    for (int i = 0; i < nDofs; i++)
588
      bary.push_back(basFcts->getCoords(i));    
589
590

    // traverse mesh
591
    std::vector<bool> visited(getUsedSize(), false);
592
    TraverseStack stack;
593
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
594
595
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
596
    ElementVector localUh(basFcts->getNumber());
597

598
    while (elInfo) {
599
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
600
      getLocalVector(elInfo->getElement(), localUh);
601
602
603
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
604
      for (int i = 0; i < nNodes; i++) { // for all nodes
605
606
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
607

608
	  if (!visited[dofIndex]) {
609
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, grd);
610

611
612
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
613
614
615
616
617
618
619
620
621
622
623
624
625
626

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


627
628
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
629
  {}
630

631

632
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
633
634
  {}

635
636
637
638

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

641
    TEST_EXIT_DBG(vec)("no vector\n");
642

643
    int dow = Global::getGeo(WORLD);
644
    static WorldVector<DOFVector<double>*> *result = nullptr; // TODO: REMOVE STATIC
645

646
    if (!res && !result) {
647
648
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
649
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
650
651
652
653
654
    }

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

    int vecSize = vec->getSize();
655
656
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
657
658
659
660
661
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

662
663


Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
666
667
668
669
670
671
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    if (!(op || this->operators.size())) 
      return;

672
    set_to_zero(this->elementVector);
673
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
674
675
676

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
677
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
678
679
680
681
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

682
      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
683
	   it != this->operators.end(); 
684
	   ++it, ++factorIt)
685
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
686
	  (*it)->getElementVector(elInfo, this->elementVector, 
687
				  *factorIt ? **factorIt : 1.0);      
688
	  addVector = true;	  
689
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
690
691
    }

692
693
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
694
695
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
696

Thomas Witkowski's avatar
Thomas Witkowski committed
697
698
699
700
701
702
703
704
705
706
707
708
  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;

709
    set_to_zero(this->elementVector);
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
    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
725
726
      }
  
727
728
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
729
730
731
  }


732
733
734
735
736
737
738
739
740
  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
741
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
742
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
743
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
744
    FastQuadrature *fastQuad =
Praetorius, Simon's avatar
Praetorius, Simon committed
745
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
746
747
748
749
750
751
752
753
754

    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
755
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
756
757
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
758
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780

      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
781
				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > > *fct)
782
783
784
785
786
787
788
789
  {
    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
790
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
791
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
792
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
793
    FastQuadrature *fastQuad =
794
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad,  INIT_PHI | INIT_GRD_PHI);
795
796
797
798
799
800
801
802
803
804
805

    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;
806
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;
807
    TraverseStack stack;
Praetorius, Simon's avatar
Praetorius, Simon committed
808
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
809
810
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
811
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
812
      for (size_t i = 0; i < grds.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
813
	grds[i]->getGrdAtQPs(elInfo, quad, fastQuad, qpGrd[i]);
814
815
816
817
818
819
820
821
822
823
824
825
826
827

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

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

833
834
    return value;
  }
835
836
}