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


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

namespace AMDiS {

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

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

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

55

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

70

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

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

90

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

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

131

132
    if (oldElInfo == NULL)
133
134
      delete elInfo;

135
136
137
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(value);
#endif
138
    return value;
139
  }
140
141
142


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

    if (oldElInfo)
      oldElInfo = elInfo;

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

177
    if (oldElInfo == NULL)
178
179
      delete elInfo;

180
181
    return value;
  }
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
  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;
233
          this->getVecAtQPs(elInfo, quadSurfaces[face], NULL, uh_vec);
234
235
236
237
238
239
240
241
242
          for (int iq = 0; iq < nPoints; iq++)
            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
          result += det * normT;
        }
      }
  
      elInfo = stack.traverseNext(elInfo);
    }
  
243
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
244
    Parallel::mpi::globalAdd(result);
245
#endif
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
    
    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);
303
          this->getVecAtQPs(elInfo, quadSurfaces[face], NULL, uh_vec);
304
305
306
307
308
309
310
311
312
313
          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);
    }
  
314
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
315
    Parallel::mpi::globalAdd(result);
316
#endif
317
318
319
320
321
    
    return result;
  }


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

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

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

340
341
    Element *el = elInfo->getElement(); 

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

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

352
    d2AtQPs.change_dim(nPoints);
353
    if (quadFast) {
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
360
361
362
	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);
363
364
	}

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

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

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

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

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

402

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

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

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

    this->set(0.0);

416
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
417
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
418

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

426
      if (basisFcts->isNodal()) {
427
428
	while (elInfo) {
	  Element *el = elInfo->getElement();
429

430
431
	  basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
	  source->getLocalVector(el, sourceLocalCoeffs);
432

433
434
435
436
437
	  for (int i = 0; i < nBasisFcts; i++) {
	    if (vec[myLocalIndices[i]] == 0.0) {
	      coords = basisFcts->getCoords(i);
	      vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
	    }
438
	  }
439
	  elInfo = stack.traverseNext(elInfo);
440
	}
441
      } else {
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
	// nonnodal base:
	ElementFunctionWorld<double> F(source, factor);
	
	while (elInfo) {
	  F.setElInfo(elInfo);
	  
	  Element *el = elInfo->getElement();
	  
	  basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
	  
	  ElementVector rvec(nBasisFcts);	  
	  basisFcts->interpol(elInfo, nBasisFcts, NULL, &F, rvec);
	  
	  for (int i = 0; i < nBasisFcts; i++) {
	    if (vec[myLocalIndices[i]] == 0.0) {
	      vec[myLocalIndices[i]] = rvec[i];
	    }
	  }
	  elInfo = stack.traverseNext(elInfo);
	}	
      }
    } else {
      if (basisFcts->isNodal()) {
465
466
467
468
469
470
471
472
473
474
475
476
477
      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);
478
      while (nextTraverse) {     
479
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
480
				   myLocalIndices);
481
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
482

483
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
484
	  if (vec[myLocalIndices[i]] == 0.0) {
485
            elInfo1->coordToWorld(*(basisFcts->getCoords(i)), worldVec);
486
487
488
	    elInfo2->worldToCoord(worldVec, &coords2);
	  
	    bool isPositive = true;
489
	    for (int j = 0; j < coords2.size(); j++) {
490
	      if (coords2[j] < -0.00001) {
491
492
493
494
495
		isPositive = false;
		break;
	      }
	    }
	  
496
497
498
	    if (isPositive)
	      vec[myLocalIndices[i]] = 
		sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	   
499
500
501
	  }
	}
      
502
503
	nextTraverse = 
	  dualStack.traverseNext(&elInfo1, &elInfo2, &elInfoSmall, &elInfoLarge);
504
      }
505
506
507
508
      } else {
	// nonnodal base
	ERROR_EXIT("not yet implemented\n");
      }
509
510
511
    }
  }

512
513

  template<>
514
515
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
516
  {
517
    FUNCNAME("DOFVector<WorldVector<double> >::interpol()");
518
519
520
521
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

522
    DimVec<double> *coords = NULL;
523
    const FiniteElemSpace *vFeSpace = v->getFeSpace();
524

525
    if (feSpace == vFeSpace)
526
527
528
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
529
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
530

531
    int nBasFcts = basFcts->getNumber();
532
533
    int vNumBasFcts = vBasFcts->getNumber();

534
535
    if (feSpace->getMesh() == vFeSpace->getMesh()) {      
      std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
536
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
537
538
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
539
540
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
541
      if (basFcts->isNodal()){
542
543
      while (elInfo) {
	Element *el = elInfo->getElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
544
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
545
546
	v->getLocalVector(el, vLocalCoeffs);

547
	for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
548
	  if (vec[myLocalIndices[i]] == nul) {
549
	    coords = basFcts->getCoords(i);
550
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs) * factor;
551
552
553
554
555
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
	//nonnodal base:
	ElementFunctionWorld<WorldVector<double> > F(v, factor);
	
	while (elInfo) {
	  F.setElInfo(elInfo);
	  
	  Element *el = elInfo->getElement();
	  
	  basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
	  v->getLocalVector(el, vLocalCoeffs);
	  
	  mtl::dense_vector<WorldVector<double> > rvec(nBasFcts);
	  basFcts->interpol(elInfo, nBasFcts, NULL, &F, rvec);
	  
	  for (int i = 0; i < nBasFcts; i++) {
	    if (vec[myLocalIndices[i]] == nul) {
	      vec[myLocalIndices[i]] = rvec[i];
	    }
	  }
	  elInfo = stack.traverseNext(elInfo);
	}	
      }
    } else {
      ERROR_EXIT("not yet implemented for dual traverse\n");
580
581
582
583
    }
  }


584
585
586
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
587
    FUNCNAME_DBG("DOFVector<double>::getGradient()");
588
589
590
591
592
593
594
595
596
597

    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
598
    static WorldVector<DOFVector<double>*> *result = NULL; // TODO: REMOVE STATIC
599

600
    if (grad) {
601
602
      result = grad;
    } else {
603
      if (!result) {
604
	result = new WorldVector<DOFVector<double>*>;
605

606
	result->set(NULL);
607
      }
608
      for (int i = 0; i < dow; i++) {
609
	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {
610
611
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
612
613
614
615
616
	}
      }
    }

    // count number of nodes and dofs per node
617
618
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
619
    std::vector<DimVec<double>*> bary;
620

621
622
    int nNodes = 0;
    int nDofs = 0;
623

624
    for (int i = 0; i < dim + 1; i++) {
625
626
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
627
      int numPreDofs = admin->getNumberOfPreDofs(i);
628
      for (int j = 0; j < numPositionNodes; j++) {
629
	int dofs = basFcts->getNumberOfDofs(geoIndex);
630
	nNodeDOFs.push_back(dofs);
631
	nDofs += dofs;
632
	nNodePreDofs.push_back(numPreDofs);
633
      }
634
      nNodes += numPositionNodes;
635
636
    }

637
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
638
      ("number of dofs != number of basis functions\n");
639
    
640
    for (int i = 0; i < nDofs; i++)
641
      bary.push_back(basFcts->getCoords(i));    
642
643

    // traverse mesh
644
    std::vector<bool> visited(getUsedSize(), false);
645
    TraverseStack stack;
646
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
647
648
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
649
    ElementVector localUh(basFcts->getNumber());
650

651
    while (elInfo) {
652
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
653
      getLocalVector(elInfo->getElement(), localUh);
654
655
656
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
657
      for (int i = 0; i < nNodes; i++) { // for all nodes
658
659
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
660

661
	  if (!visited[dofIndex]) {
662
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, grd);
663

664
665
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
666
667
668
669
670
671
672
673
674
675
676
677
678
679

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


680
681
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
682
  {}
683

684

685
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
686
687
  {}

688
689
690
691

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

694
    TEST_EXIT_DBG(vec)("no vector\n");
695

696
    int dow = Global::getGeo(WORLD);
697
    static WorldVector<DOFVector<double>*> *result = NULL; // TODO: REMOVE STATIC
698

699
    if (!res && !result) {
700
701
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
702
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
703
704
705
706
707
    }

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

    int vecSize = vec->getSize();
708
709
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
710
711
712
713
714
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

715
716


Thomas Witkowski's avatar
Thomas Witkowski committed
717
718
719
720
721
722
723
724
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    if (!(op || this->operators.size())) 
      return;

725
    set_to_zero(this->elementVector);
726
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
727
728
729

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
730
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
731
732
733
734
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

735
      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
736
	   it != this->operators.end(); 
737
	   ++it, ++factorIt)
738
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
739
	  (*it)->getElementVector(elInfo, this->elementVector, 
740
				  *factorIt ? **factorIt : 1.0);      
741
	  addVector = true;	  
742
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
    }

745
746
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
749

Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
752
753
754
755
756
757
758
759
760
761
  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;

762
    set_to_zero(this->elementVector);
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
    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
778
779
      }
  
780
781
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
782
783
  }

784
785
786
787
788
  
  template<>
  void DOFVectorBase<double>::assembleOperator(Operator &op)
  {
    FUNCNAME("DOFVectorBase::assembleOperator()");
Thomas Witkowski's avatar
Thomas Witkowski committed
789

790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
    TEST_EXIT(op.getRowFeSpace() == feSpace)
      ("Row FE spaces do not fit together!\n");

    Mesh *mesh = feSpace->getMesh();
    mesh->dofCompress();
    const BasisFunction *basisFcts = feSpace->getBasisFcts();

    Flag assembleFlag = getAssembleFlag() | 
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_NEIGH |
      Mesh::FILL_BOUND;

    BoundaryType *bound = new BoundaryType[basisFcts->getNumber()];

    if (getBoundaryManager())
      getBoundaryManager()->initVector(this);


    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
    while (elInfo) {
      basisFcts->getBound(elInfo, bound);

      assemble(1.0, elInfo, bound, &op);

      if (getBoundaryManager())
	getBoundaryManager()->fillBoundaryConditions(elInfo, this);

      elInfo = stack.traverseNext(elInfo);
    }

    finishAssembling();
    getBoundaryManager()->exitVector(this);

    delete [] bound;
  }
  
  
831
832
833
834
835
836
837
838
839
  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
840
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
841
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
842
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
843
    FastQuadrature *fastQuad =
Praetorius, Simon's avatar
Praetorius, Simon committed
844
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
845
846
847
848
849
850
851
852
853

    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
854
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
855
856
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
857
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
858
859
860
861
862
863
864
865
866
867
868
869
870

      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
871
    Parallel::mpi::globalAdd(value);
872
873
874
875
876
877
878
#endif

    return value;
  }

  double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
				  const std::vector<DOFVector<double>*> &grds,
Praetorius, Simon's avatar
Praetorius, Simon committed
879
				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > > *fct)
880
881
882
883
884
885
886
887
  {
    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
888
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
889
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
890
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
891
    FastQuadrature *fastQuad =
892
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad,  INIT_PHI | INIT_GRD_PHI);
893
894
895
896
897
898
899
900
901
902
903

    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;
904
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;
905
    TraverseStack stack;
Praetorius, Simon's avatar
Praetorius, Simon committed
906
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
907
908
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
909
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
910
      for (size_t i = 0; i < grds.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
911
	grds[i]->getGrdAtQPs(elInfo, quad, fastQuad, qpGrd[i]);
912
913
914
915
916
917
918
919
920
921
922
923
924
925

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

926
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
927
    Parallel::mpi::globalAdd(value);
928
929
#endif

930
931
    return value;
  }
932
933
}