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
    DOFVector<double>* null_ptr = NULL;
600

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

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

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

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

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

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

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


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

685

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

689
690
691
692

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

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

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

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

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

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

    return r;
  }

716
717


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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
750

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

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

785
786
787
788
789
  
  template<>
  void DOFVectorBase<double>::assembleOperator(Operator &op)
  {
    FUNCNAME("DOFVectorBase::assembleOperator()");
Thomas Witkowski's avatar
Thomas Witkowski committed
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
831
    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;
  }
  
  
832
833
834
835
836
837
838
839
840
  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
841
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
842
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
843
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
844
    FastQuadrature *fastQuad =
Praetorius, Simon's avatar
Praetorius, Simon committed
845
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
846
847
848
849
850
851
852
853
854

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

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

    return value;
  }

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

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

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

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

931
932
    return value;
  }
933
934
}