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

Praetorius, Simon's avatar
Praetorius, Simon committed
106
    delete [] localIndices;
107
108
109
    if (oldElInfo == NULL)
      delete elInfo;

Praetorius, Simon's avatar
Praetorius, Simon committed
110
111
    if (values != NULL)
      *values = value;
112
    return value;
Praetorius, Simon's avatar
Praetorius, Simon committed
113
  };
114
115
116


  template<>
Praetorius, Simon's avatar
Praetorius, Simon committed
117
  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
118
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
119
    FUNCNAME("DOFVector<double>::evalAtCoords()");
120
121
122
123
124
125
126
  
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Praetorius, Simon's avatar
Praetorius, Simon committed
127
    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
128
129
130
131
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();

Praetorius, Simon's avatar
Praetorius, Simon committed
132
133
134
    static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
    WorldVector<double> *val = (NULL != values) ? values : &Values;

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
159
160
    return ((*val));
  };
161
162


163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
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
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
  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);
    }
  
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
    #endif
    
    return result;
  }


  template<>
  double DOFVector<WorldVector<double> >::IntOnBoundaryNormal(
    BoundaryType boundaryType, Quadrature* q) const
  {
    FUNCNAME("DOFVector<T>::IntOnBoundary()");
 
    Mesh* mesh = this->feSpace->getMesh();
    int dim = mesh->getDim();
  
    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(dim - 1, deg);
    } else {
      TEST_EXIT(q->getDim() == dim-1)
        ("Provided quadrature has wrong dimension for surfaceQuadrature!\n");
    }
  
    // create barycentric coords for each vertex of each side
    const Element *refElement = Global::getReferenceElement(dim);
    VectorOfFixVecs<DimVec<double> >**coords;
    coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
  
    // for all element sides
    for (int i = 0; i < dim + 1; i++) {
      coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
        DimVec<double>(dim, DEFAULT_VALUE, 0.0));
      // for each vertex of the side
      for (int k = 0; k < dim; k++) {
        int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), 
                                                    i, k);
        (*coords[i])[k][index] = 1.0;
      }
    }
  
    std::vector<SurfaceQuadrature*> quadSurfaces(dim + 1);
    for (int i = 0; i < dim + 1; i++)
      quadSurfaces[i] = new SurfaceQuadrature(q, *coords[i]);
  
    double result = 0.0;
    WorldVector<double> normal;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
      Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
    while (elInfo) {
      for (int face = 0; face < dim + 1; face++) {
        if (elInfo->getBoundary(face) == boundaryType) {
          elInfo->getNormal(face, normal);
          double det = elInfo->calcSurfaceDet(*coords[face]);
  
          int nPoints = quadSurfaces[face]->getNumPoints();
          mtl::dense_vector<WorldVector<double> > uh_vec(nPoints);
          WorldVector<double> normT; normT.set(0.0);
          this->getVecAtQPs(elInfo, quadSurfaces[face], NULL, uh_vec);
          for (int iq = 0; iq < nPoints; iq++)
            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
          // scalar product between vector-valued solution and normal vector
          result += det * (normT * normal); 
        }
      }
  
      elInfo = stack.traverseNext(elInfo);
    }
  
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
    #endif
    
    return result;
  }


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

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

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

321
322
    Element *el = elInfo->getElement(); 

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

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

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

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

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

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

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

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

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

383

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

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

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

    this->set(0.0);

398
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
399
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
400

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

408
      while (elInfo) {
409
410
	Element *el = elInfo->getElement();

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

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

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

466
467

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

    this->set(nul);

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

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

    const BasisFunction *basFcts = feSpace->getBasisFcts();
482
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
483

484
    int nBasFcts = basFcts->getNumber();
485
486
    int vNumBasFcts = vBasFcts->getNumber();

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

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

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


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

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

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

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

551
552
    int nNodes = 0;
    int nDofs = 0;
553

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

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

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


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

614

615
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
616
617
  {}

618
619
620
621

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

624
    TEST_EXIT_DBG(vec)("no vector\n");
625

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

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

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

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

    return r;
  }

645
646


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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
681

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

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


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

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

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

    std::vector<mtl::dense_vector<double> > qp(vecs.size());
    std::vector<mtl::dense_vector<WorldVector<double> > > qpGrd(vecs.size());
    std::vector<double> qp_local(vecs.size());
    std::vector<WorldVector<double> > grd_local(grds.size());
    for (size_t i = 0; i < vecs.size(); i++)
      qp[i].change_dim(fastQuad->getNumPoints());
    for (size_t i = 0; i < grds.size(); i++)
      qpGrd[i].change_dim(fastQuad->getNumPoints());

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

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

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

818
819
    return value;
  }
820
821
}