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
  {  
74
75
76
77
78
79
    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
80
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
81
82
83
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();
84
    double value = 0.0;
85
86
87
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
Praetorius, Simon's avatar
Praetorius, Simon committed
88
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
89
90
91
92
93
94
95
96
      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
97
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
98
99
      ElementVector uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
100
101
        uh[i] = operator[](localIndices[i]);
      value = basFcts->evalUh(lambda, uh);
102
103
104
105
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      return 0.0;
#else
106
      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
107
108
109
#endif
    }

110
111
112
113
114

    if (oldElInfo == NULL)
      delete elInfo;

    return value;
115
  }
116
117
118


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

    if (oldElInfo == NULL)
      delete elInfo;

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


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

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

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

317
318
    Element *el = elInfo->getElement(); 

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

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

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

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

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

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

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

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

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

379

380
381
382
383
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {

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

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

    this->set(0.0);

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

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

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

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

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

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

461
462

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

    this->set(nul);

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

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

    const BasisFunction *basFcts = feSpace->getBasisFcts();
478
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
479

480
    int nBasFcts = basFcts->getNumber();
481
482
    int vNumBasFcts = vBasFcts->getNumber();

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

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

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


510
511
512
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
513
    FUNCNAME_DBG("DOFVector<double>::getGradient()");
514
515
516
517
518
519
520
521
522
523
524
525

    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;

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

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

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

547
548
    int nNodes = 0;
    int nDofs = 0;
549

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

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

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


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

610

611
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
612
613
  {}

614
615
616
617

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

620
    TEST_EXIT_DBG(vec)("no vector\n");
621

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

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

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

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

    return r;
  }

641
642


Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
645
646
647
648
649
650
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    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
}