DOFVector.cc 24.8 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();
  
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
82
    std::vector<DegreeOfFreedom> localIndices(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
107
108
109
110
        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."));

    if (oldElInfo == NULL)
      delete elInfo;

    return value;
111
  }
112
113
114


  template<>
115
116
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const
117
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
118
    FUNCNAME("DOFVector<double>::evalAtCoords()");
119
120
121
122
123
124
125
  
    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
126
127
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
    DimVec<double> lambda(dim, NO_INIT);  
128
    ElInfo *elInfo = mesh->createNewElInfo();
129
    WorldVector<double> value(DEFAULT_VALUE, 0.0);
130
131
132
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
Praetorius, Simon's avatar
Praetorius, Simon committed
133
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
134
135
136
137
138
139
140
141
      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
142
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
143
144
      mtl::dense_vector<WorldVector<double> > uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
145
        uh[i] = operator[](localIndices[i]);
146
      value = basFcts->evalUh(lambda, uh);
147
148
149
150
151
152
    } else
      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));

    if (oldElInfo == NULL)
      delete elInfo;

153
154
    return value;
  }
155
156


157
158
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
  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);
    }
  
216
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
217
218
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
219
#endif
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
    
    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);
    }
  
288
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
289
290
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
291
#endif
292
293
294
295
296
    
    return result;
  }


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

307
308
309
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
310
311
    }

312
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
313
314
      ("invalid basis functions");

315
316
    Element *el = elInfo->getElement(); 

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
336
337
	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);
338
339
	}

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

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
	for (int i = 0; i < nBasFcts; i++) {
358
	  WARNING("not tested after index correction\n");
359
	  (*(basFcts->getD2Phi(i)))(quad->getLambda(iq), D2Phi);
360

Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
363
	      D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
364
365
	}

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

377

378
379
380
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
381
    FUNCNAME("DOFVector<double>::interpol()");
382

383
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
384
385
386
387
388
389
390
391
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

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

    this->set(0.0);

392
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
393
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
394

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

402
      while (elInfo) {
403
404
	Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
405
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
406
407
	source->getLocalVector(el, sourceLocalCoeffs);

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

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

460
461

  template<>
462
463
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
464
465
466
467
468
469
  {
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

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

472
    if (feSpace == vFeSpace)
473
474
475
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
476
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
477

478
    int nBasFcts = basFcts->getNumber();
479
480
    int vNumBasFcts = vBasFcts->getNumber();

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

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

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


508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
  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;

524
    if (grad) {
525
526
      result = grad;
    } else {
527
      if (!result) {
528
	result = new WorldVector<DOFVector<double>*>;
529

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

    // count number of nodes and dofs per node
541
542
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
543
    std::vector<DimVec<double>*> bary;
544

545
546
    int nNodes = 0;
    int nDofs = 0;
547

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

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

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

575
    while (elInfo) {
576
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
577
      getLocalVector(elInfo->getElement(), localUh);
578
579
580
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
581
      for (int i = 0; i < nNodes; i++) { // for all nodes
582
583
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
584

585
	  if (!visited[dofIndex]) {
586
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, grd);
587

588
589
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
590
591
592
593
594
595
596
597
598
599
600
601
602
603

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


604
605
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
606
  {}
607

608

609
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
610
611
  {}

612
613
614
615

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

618
    TEST_EXIT_DBG(vec)("no vector\n");
619

620
    int dow = Global::getGeo(WORLD);
621
622
    static WorldVector<DOFVector<double>*> *result = NULL;

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

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

    int vecSize = vec->getSize();
632
633
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
634
635
636
637
638
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

639
640


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

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

651
    set_to_zero(this->elementVector);
652
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
653
654
655

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

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

671
672
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
673
674
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
675

Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
678
679
680
681
682
683
684
685
686
687
  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;

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


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

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

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

    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;
785
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;
786
    TraverseStack stack;
Praetorius, Simon's avatar
Praetorius, Simon committed
787
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
788
789
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
790
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
791
      for (size_t i = 0; i < grds.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
792
	grds[i]->getGrdAtQPs(elInfo, quad, fastQuad, qpGrd[i]);
793
794
795
796
797
798
799
800
801
802
803
804
805
806

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

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

812
813
    return value;
  }
814
815
}