DOFVector.cc 28.6 KB
Newer Older
1
#include <boost/numeric/mtl/mtl.hpp>
2
3
4
5
6
7
8
9
10
11
#include "DOFVector.h"
#include "Traverse.h"
#include "DualTraverse.h"
#include "FixVec.h"

namespace AMDiS {

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

14
    switch (coarsenOperation) {
15
16
17
18
19
20
21
    case NO_OPERATION:
      return;
      break;
    case COARSE_RESTRICT:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
      break;
    case COARSE_INTERPOL:
22
23
      TEST_EXIT_DBG(feSpace)("Should not happen!\n");
      TEST_EXIT_DBG(feSpace->getBasisFcts())("Shoud not happen!\n");
24

25
26
27
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
      break;
    default:
28
29
      WARNING("Invalid coarsen operation \"%d\" in vector \"%s\"\n", 
	      coarsenOperation, this->name.c_str());
30
31
32
    }
  }

33

34
35
36
37
38
39
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList& list, int n)
  {
    (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->refineInter(this, &list, n);
  }

40

41
42
43
  template<>
  void DOFVector<WorldVector<double> >::refineInterpol(RCNeighbourList& list, int n)
  {
44
45
46
    if (n < 1) 
      return;

47
48
    Element *el = list.getElement(0);
    int n0 = feSpace->getAdmin()->getNumberOfPreDOFs(VERTEX);
49
50
51
    DegreeOfFreedom dof0 = el->getDof(0, n0);
    DegreeOfFreedom dof1 = el->getDof(1, n0);
    DegreeOfFreedom dof_new = el->getChild(0)->getDof(feSpace->getMesh()->getDim(), n0);
52
53
54
55
56
    vec[dof_new] = vec[dof0];
    vec[dof_new] += vec[dof1];
    vec[dof_new] *= 0.5;
  }

57

58
59
60
61
62
63
64
65
66
  template<>
  DOFVector<WorldVector<double> >*
  DOFVector<double>::getGradient(DOFVector<WorldVector<double> > *grad) const 
  {
    FUNCNAME("DOFVector<double>::getGradient()");

    // define result vector
    static DOFVector<WorldVector<double> > *result = NULL;

Thomas Witkowski's avatar
Thomas Witkowski committed
67
    if (grad) {
68
69
      result = grad;
    } else {
70
      if(result && result->getFeSpace() != feSpace) {
71
72
	delete result;
	result = new DOFVector<WorldVector<double> >(feSpace, "gradient");
73
74
75
76
77
78
79
80
81
      }
    }

    Mesh *mesh = feSpace->getMesh();
    int dim = mesh->getDim();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    DOFAdmin *admin = feSpace->getAdmin();

    // count number of nodes and dofs per node
82
83
84
    std::vector<int> numNodeDOFs;
    std::vector<int> numNodePreDOFs;
    std::vector<DimVec<double>*> bary;
85

86
87
    int nNodes = 0;
    int nDofs = 0;
88

Thomas Witkowski's avatar
Thomas Witkowski committed
89
    for (int i = 0; i < dim + 1; i++) {
90
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
91
      int nPositions = mesh->getGeo(geoIndex);
92
      int numPreDOFs = admin->getNumberOfPreDOFs(i);
93
      for (int j = 0; j < nPositions; j++) {
94
95
	int dofs = basFcts->getNumberOfDOFs(geoIndex);
	numNodeDOFs.push_back(dofs);
96
	nDofs += dofs;
97
98
	numNodePreDOFs.push_back(numPreDOFs);
      }
99
      nNodes += nPositions;
100
101
    }

102
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
103
      ("number of dofs != number of basis functions\n");
104
    
105
    for (int i = 0; i < nDofs; i++)
106
107
      bary.push_back(basFcts->getCoords(i));    

108
    ElementVector localUh(basFcts->getNumber());
109
110

    // traverse mesh
111
    std::vector<bool> visited(getUsedSize(), false);
112
    TraverseStack stack;
113
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
114
115
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

Thomas Witkowski's avatar
Thomas Witkowski committed
116
    while (elInfo) {
117
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
118
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
119
      getLocalVector(elInfo->getElement(), localUh);
120
121

      int localDOFNr = 0;
122
      for (int i = 0; i < nNodes; i++) { // for all nodes
Thomas Witkowski's avatar
Thomas Witkowski committed
123
	for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
124
	  DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[i] + j];
Thomas Witkowski's avatar
Thomas Witkowski committed
125
126
127
	  if (!visited[dofIndex]) {	  
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, 
			       localUh, &((*result)[dofIndex]));
128
129
130
131
132
133
134
135
136
137
138
139
140

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }

141

142
143
144
145
146
147
  template<>
  DOFVector<WorldVector<double> >*
  DOFVector<double>::getRecoveryGradient(DOFVector<WorldVector<double> > *grad) const 
  {
    FUNCNAME("DOFVector<double>::getRecoveryGradient()");

148
149
150
151
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif

152
153
154
155
156
    // define result vector
    static DOFVector<WorldVector<double> > *vec = NULL;

    DOFVector<WorldVector<double> > *result = grad;

Thomas Witkowski's avatar
Thomas Witkowski committed
157
    if (!result) {
158
      if (vec && vec->getFeSpace() != feSpace) {
159
	delete vec;
160
161
	vec = NULL;
      }
162
      if (!vec)
163
	vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
164
      
165
166
167
168
169
170
171
172
173
      result = vec;
    }

    result->set(WorldVector<double>(DEFAULT_VALUE, 0.0));

    DOFVector<double> volume(feSpace, "volume");
    volume.set(0.0);

    const BasisFunction *basFcts = feSpace->getBasisFcts();
Thomas Witkowski's avatar
Thomas Witkowski committed
174
    int nPreDOFs = feSpace->getAdmin()->getNumberOfPreDOFs(0);
175
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
Thomas Witkowski's avatar
Thomas Witkowski committed
176
    WorldVector<double> grd;
177
178

    // traverse mesh
Thomas Witkowski's avatar
Thomas Witkowski committed
179
    Mesh *mesh = feSpace->getMesh();
180
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
181
182
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET | 
183
					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
184

185
    ElementVector localUh(basFcts->getNumber());
186

Thomas Witkowski's avatar
Thomas Witkowski committed
187
    while (elInfo) {
188
      double det = elInfo->getDet();
189
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
190
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
191
      getLocalVector(elInfo->getElement(), localUh);
Thomas Witkowski's avatar
Thomas Witkowski committed
192
      basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
193

Thomas Witkowski's avatar
Thomas Witkowski committed
194
      for (int i = 0; i < dim + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
195
	DegreeOfFreedom dofIndex = dof[i][nPreDOFs];
196
197
198
199
200
201
202
203
204
205
	(*result)[dofIndex] += grd * det;
	volume[dofIndex] += det;
      }

      elInfo = stack.traverseNext(elInfo);
    }

    DOFVector<double>::Iterator volIt(&volume, USED_DOFS);
    DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS);

206
207
    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
      if (*volIt != 0.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
208
	*grdIt *= 1.0 / (*volIt);
209
210
211
212

    return result;
  }

213

214
  template<>
215
216
  const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo, 
								const Quadrature *quad,
217
								const FastQuadrature *quadFast,
218
								WorldVector<double> *grdAtQPs) const
219
220
221
  {
    FUNCNAME("DOFVector<double>::getGrdAtQPs()");
  
222
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
223

224
225
226
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
227
228
    }

229
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
230
231
232
      ("invalid basis functions");

    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
233
234
    int myRank = omp_get_thread_num();
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
235
236
237
238
239
240
    static WorldVector<double> *grd = NULL;
    WorldVector<double> *result;

    if (grdAtQPs) {
      result = grdAtQPs;
    } else {
241
242
243
      if (grd) 
	delete [] grd;

Thomas Witkowski's avatar
Thomas Witkowski committed
244
      grd = new WorldVector<double>[nPoints];
245
      for (int i = 0; i < nPoints; i++)
246
	grd[i].set(0.0);
247
      
248
249
250
      result = grd;
    }
  
251
252
    const ElementVector& localVec = localVectors[myRank];
    getLocalVector(elInfo->getElement(), const_cast<ElementVector&>(localVec));
253

Thomas Witkowski's avatar
Thomas Witkowski committed
254
    DimVec<double> &grd1 = *grdTmp[myRank];
255
256
257
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

258
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
259
      for (int i = 0; i < nPoints; i++) {
260
	for (int j = 0; j < dim + 1; j++)
261
262
	  grd1[j] = 0.0;

263
264
	for (int j = 0; j < nBasFcts; j++)
	  for (int k = 0; k < parts; k++)
265
	    grd1[k] += quadFast->getGradient(i, j, k) * localVec[j];
266
	  
267
268
	for (int l=0; l < dow; l++) {
	  result[i][l] = 0.0;
269
	  for (int k = 0; k < parts; k++)
270
271
272
	    result[i][l] += grdLambda[k][l] * grd1[k];
	}
      }
273
274

    } else {      
Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
      const BasisFunction *basFcts = feSpace->getBasisFcts();
      DimVec<double>* grdPhi = grdPhis[myRank];
277

Thomas Witkowski's avatar
Thomas Witkowski committed
278
      for (int i = 0; i < nPoints; i++) {
279
	for (int j = 0; j < dim + 1; j++)
280
	  grd1[j] = 0.0;
281
282

	for (int j = 0; j < nBasFcts; j++) {
283
	  (*(basFcts->getGrdPhi(j)))(quad->getLambda(i), *grdPhi);
284
	  for (int k = 0; k < parts; k++)
285
	    grd1[k] += (*grdPhi)[k] * localVec[j];
286
287
	}

288
289
	for (int l = 0; l < dow; l++) {
	  result[i][l] = 0.0;
290
	  for (int k = 0; k < parts; k++)
291
292
293
294
	    result[i][l] += grdLambda[k][l] * grd1[k];
	}
      }
    }
295

296
297
298
    return const_cast<const WorldVector<double>*>(result);
  }

299

300
  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
301
302
303
304
305
306
307
308
309
310
  const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *smallElInfo,
								const ElInfo *largeElInfo,
								const Quadrature *quad,
								const FastQuadrature *quadFast,
								WorldVector<double> *grdAtQPs) const
  {
    FUNCNAME("DOFVector<double>::getGrdAtQPs()");
  
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");

311
    if (quad && quadFast) {
312
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
313
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330

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

    int dow = Global::getGeo(WORLD);
    int myRank = omp_get_thread_num();
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
    static WorldVector<double> *grd = NULL;
    WorldVector<double> *result;

    if (grdAtQPs) {
      result = grdAtQPs;
    } else {
      if (grd) 
	delete [] grd;

      grd = new WorldVector<double>[nPoints];
331
      for (int i = 0; i < nPoints; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
332
	grd[i].set(0.0);
333

Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
336
      result = grd;
    }
  
337
338
    const ElementVector& localVec = localVectors[myRank];
    getLocalVector(largeElInfo->getElement(), const_cast<ElementVector&>(localVec));
339
340

    const BasisFunction *basFcts = feSpace->getBasisFcts();    
341
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
344
345
346
347
348
349
350

    DimVec<double> &grd1 = *grdTmp[myRank];
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();

    DimVec<double>* grdPhi = grdPhis[myRank];
    DimVec<double> tmp(dim, DEFAULT_VALUE, 0.0);
    
    for (int i = 0; i < nPoints; i++) {
351
      for (int j = 0; j < dim + 1; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
354
355
356
357
358
	grd1[j] = 0.0;
      
      for (int j = 0; j < nBasFcts; j++) {

	grdPhi->set(0.0);
	for (int k = 0; k < nBasFcts; k++) {
	  (*(basFcts->getGrdPhi(k)))(quad->getLambda(i), tmp);
359
	  tmp *= m[j][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
362
	  *grdPhi += tmp;
	}

363
	for (int k = 0; k < parts; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
366
367
368
	  grd1[k] += (*grdPhi)[k] * localVec[j];
      }
      
      for (int l = 0; l < dow; l++) {
	result[i][l] = 0.0;
369
	for (int k = 0; k < parts; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
372
373
374
375
376
	  result[i][l] += grdLambda[k][l] * grd1[k];
      }
    }    
  
    return const_cast<const WorldVector<double>*>(result);
  }

377

Thomas Witkowski's avatar
Thomas Witkowski committed
378
379
380
  template<>
  const WorldMatrix<double> *DOFVectorBase<double>::getD2AtQPs(const ElInfo *elInfo, 
							       const Quadrature *quad,
381
382
383
384
385
							       const FastQuadrature *quadFast,
							       WorldMatrix<double>  *d2AtQPs) const
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
386
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
387

388
389
390
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
391
392
    }

393
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
394
395
      ("invalid basis functions");

396
397
    Element *el = elInfo->getElement(); 

Thomas Witkowski's avatar
Thomas Witkowski committed
398
    int myRank = omp_get_thread_num();
399
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
400
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
401
402
403
404
405
406
407
408

    static WorldMatrix<double> *vec = NULL;
    WorldMatrix<double> *result;

    if (d2AtQPs) {
      result = d2AtQPs;
    } else {
      if(vec) delete [] vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
409
410
      vec = new WorldMatrix<double>[nPoints];
      for (int i = 0; i < nPoints; i++) {
411
412
413
414
415
	vec[i].set(0.0);
      }
      result = vec;
    }
  
416
417
    const ElementVector& localVec = localVectors[myRank];
    getLocalVector(el, const_cast<ElementVector&>(localVec));
418
419
420
421
422

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

423
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
424
425
426
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
427
428
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
429
430
431
432
	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);
433
434
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
435
436
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
437
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
438
439
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
440
441
442
443
		result[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
444
      const BasisFunction *basFcts = feSpace->getBasisFcts();
445
446
      DimMat<double>* D2Phi = D2Phis[omp_get_thread_num()];

Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
449
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
450
451
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
458
	      D2Tmp[k][l] += localVec[i] * (*D2Phi)[k][l];
459
460
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
463
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
464
465
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
466
		result[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
467
468
469
	  }
      }
    }
470

471
472
473
    return const_cast<const WorldMatrix<double>*>(result);
  }

474

475
476
477
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
478
    FUNCNAME("DOFVector<double>::interpol()");
479

480
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
481
482
483
484
485
486
487
488
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

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

    this->set(0.0);

Thomas Witkowski's avatar
Thomas Witkowski committed
489
    DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
490
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
491

492
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
493
494
495
496
497
498
      DimVec<double> *coords = NULL;
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

499
      while (elInfo) {
500
501
	Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
502
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
503
504
	source->getLocalVector(el, sourceLocalCoeffs);

505
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
506
	  if (vec[myLocalIndices[i]] == 0.0) {
507
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
508
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
	  }
	}
	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);
527
      while (nextTraverse) {     
528
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
529
				   myLocalIndices);
530
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
531

532
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
533
	  if (vec[myLocalIndices[i]] == 0.0) {
534
            elInfo1->coordToWorld(*(basisFcts->getCoords(i)), worldVec);
535
536
537
	    elInfo2->worldToCoord(worldVec, &coords2);
	  
	    bool isPositive = true;
538
	    for (int j = 0; j < coords2.size(); j++) {
539
	      if (coords2[j] < -0.00001) {
540
541
542
543
544
		isPositive = false;
		break;
	      }
	    }
	  
545
546
547
	    if (isPositive)
	      vec[myLocalIndices[i]] = 
		sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	   
548
549
550
	  }
	}
      
551
552
	nextTraverse = 
	  dualStack.traverseNext(&elInfo1, &elInfo2, &elInfoSmall, &elInfoLarge);
553
554
555
556
      }
    }
  }

557
558

  template<>
559
560
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
561
562
563
564
565
566
  {
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

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

569
    if (feSpace == vFeSpace)
570
571
572
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
573
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
574
575
576
577

    int numBasFcts = basFcts->getNumber();
    int vNumBasFcts = vBasFcts->getNumber();

578
    if (feSpace->getMesh() == vFeSpace->getMesh()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
579
      DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
580
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
581
582
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
583
584
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
585
586
587

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

	for (int i = 0; i < numBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
592
	  if (vec[myLocalIndices[i]] == nul) {
593
	    coords = basFcts->getCoords(i);
594
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs, NULL) * factor;
595
596
597
598
599
600
601
602
603
604
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


605
606
607
608
609
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
    FUNCNAME("DOFVector<double>::getGradient()");

610
611
612
613
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif

614
615
616
617
618
619
620
621
622
623
624
    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;

625
    if (grad) {
626
627
      result = grad;
    } else {
628
      if (!result) {
629
	result = new WorldVector<DOFVector<double>*>;
630

631
632
	result->set(NULL);
      }
633
      for (int i = 0; i < dow; i++) {
634
	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {
635
636
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
637
638
639
640
641
	}
      }
    }

    // count number of nodes and dofs per node
642
643
644
    std::vector<int> numNodeDOFs;
    std::vector<int> numNodePreDOFs;
    std::vector<DimVec<double>*> bary;
645

646
647
    int nNodes = 0;
    int nDofs = 0;
648

649
    for (int i = 0; i < dim + 1; i++) {
650
651
652
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
      int numPreDOFs = admin->getNumberOfPreDOFs(i);
653
      for (int j = 0; j < numPositionNodes; j++) {
654
655
	int dofs = basFcts->getNumberOfDOFs(geoIndex);
	numNodeDOFs.push_back(dofs);
656
	nDofs += dofs;
657
658
	numNodePreDOFs.push_back(numPreDOFs);
      }
659
      nNodes += numPositionNodes;
660
661
    }

662
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
663
      ("number of dofs != number of basis functions\n");
664
    
665
    for (int i = 0; i < nDofs; i++)
666
      bary.push_back(basFcts->getCoords(i));    
667
668

    // traverse mesh
669
    std::vector<bool> visited(getUsedSize(), false);
670
    TraverseStack stack;
671
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
672
673
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
674
    ElementVector localUh(basFcts->getNumber());
675

676
    while (elInfo) {
677
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
678
      getLocalVector(elInfo->getElement(), localUh);
679
680
681
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
682
      for (int i = 0; i < nNodes; i++) { // for all nodes
683
	for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
684
	  DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[i] + j];
685

686
	  if (!visited[dofIndex]) {
687
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
688

689
690
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
691
692
693
694
695
696
697
698
699
700
701
702
703
704

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


705
706
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
707
  {}
708

709

710
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
711
712
  {}

713
714
715
716

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

719
    TEST_EXIT_DBG(vec)("no vector\n");
720

721
    int dow = Global::getGeo(WORLD);
722
723
    static WorldVector<DOFVector<double>*> *result = NULL;

724
    if (!res && !result) {
725
726
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
727
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
728
729
730
731
732
    }

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

    int vecSize = vec->getSize();
733
734
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
735
736
737
738
739
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

740

Thomas Witkowski's avatar
Thomas Witkowski committed
741
  template<>
742
743
744
745
746
  void DOFVectorBase<double>::getVecAtQPs(const ElInfo *smallElInfo, 
					  const ElInfo *largeElInfo,
					  const Quadrature *quad,
					  const FastQuadrature *quadFast,
					  mtl::dense_vector<double>& vecAtQPs) const
Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
  {
    FUNCNAME("DOFVector<double>::getVecAtQPs()");
749
 
Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");

752
    if (quad && quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
753
754
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
755
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
756
757
758
759

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

760
    if (smallElInfo->getMesh() == feSpace->getMesh())
761
      return getVecAtQPs(smallElInfo, quad, quadFast, vecAtQPs);    
762

Thomas Witkowski's avatar
Thomas Witkowski committed
763
    const BasisFunction *basFcts = feSpace->getBasisFcts();
764
765
    int nPoints = 
      quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
766

767
768
769
770
    vecAtQPs.change_dim(nPoints);
    
    const ElementVector &localVec = localVectors[omp_get_thread_num()];
    getLocalVector(largeElInfo->getElement(), const_cast<ElementVector&>(localVec));
Thomas Witkowski's avatar
Thomas Witkowski committed
771

772
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
773
774

    for (int i = 0; i < nPoints; i++) {
775
      vecAtQPs[i] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
776
777
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
778
	for (int k = 0; k < nBasFcts; k++)
779
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
780
781

	vecAtQPs[i] += localVec[j] * val;
Thomas Witkowski's avatar
Thomas Witkowski committed
782
783
784
785
      }
    }
  }

786

Thomas Witkowski's avatar
Thomas Witkowski committed
787
788
789
790
791
792
793
794
795
796
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    FUNCNAME("DOFVector::assemble()");

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

797
    set_to_zero(this->elementVector);
798
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
799
800
801

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
802
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
803
804
805
806
807
808
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	   it != this->operators.end(); 
809
	   ++it, ++factorIt)
810
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
811
	  (*it)->getElementVector(elInfo, this->elementVector, 
812
				  *factorIt ? **factorIt : 1.0);      
813
814
	  addVector = true;
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
815
816
    }

817
818
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
819
820
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
821

Thomas Witkowski's avatar
Thomas Witkowski committed
822
823
824
825
826
827
828
829
830
831
832
833
  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;

834
    set_to_zero(this->elementVector);
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
    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
850
851
      }
  
852
853
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
854
855
856
  }


857
858
859
860
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct)
  {
861
    if (vec1.getFeSpace()->getMesh() == vec2.getFeSpace()->getMesh())
862
863
864
865
866
867
868
869
870
871
872
873
874
875
      return intSingleMesh(vec1, vec2, fct);
    else
      return intDualMesh(vec1, vec2, fct);
  }


  double intSingleMesh(const DOFVector<double> &vec1,
		       const DOFVector<double> &vec2,
		       BinaryAbstractFunction<double, double, double> *fct)
  {
    FUNCNAME("intDualmesh()");
    
    TEST_EXIT(fct)("No function defined!\n");

876
877
    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
878
    Quadrature* quad = 
879
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
880
    FastQuadrature *fastQuad =
881
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
882

883
884
    mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
885
886
887
888

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
889
    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
890
    while (elInfo) {
891
892
893
894
      vec1.getVecAtQPs(elInfo, quad, fastQuad, qp1);
      vec2.getVecAtQPs(elInfo, quad, fastQuad, qp2);

      double tmp = 0.0;
895
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
896
897
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);     
      value += tmp * elInfo->getDet();
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
      
      elInfo = stack.traverseNext(elInfo);
    }

    return value;
  }


  double intDualMesh(const DOFVector<double> &vec1,
		     const DOFVector<double> &vec2,
		     BinaryAbstractFunction<double, double, double> *fct)
  {
    FUNCNAME("intDualmesh()");
    
    TEST_EXIT(fct)("No function defined!\n");

914
915
    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
916
    Quadrature* quad = 
917
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
918
    FastQuadrature *fastQuad =
919
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
920

921
922
    mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
923
924
925
926
927

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    DualTraverse dualTraverse;    
    DualElInfo dualElInfo;
928
929
    bool cont = dualTraverse.traverseFirst(vec1.getFeSpace()->getMesh(),
					   vec2.getFeSpace()->getMesh(),
930
931
932
					   -1, -1, traverseFlag, traverseFlag,
					   dualElInfo);
    while (cont) {
933
934
935
936
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);

      double tmp = 0.0;
937
      for (int iq = 0; iq < quad->getNumPoints(); iq++)
938
939
 	tmp += quad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);      
      value += tmp * dualElInfo.smallElInfo->getDet();
940
941
942
943
944
945
946
      
      cont = dualTraverse.traverseNext(dualElInfo);
    }

    return value;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
  
  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate()");
    
    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(), 
		       2 * vec.getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad = 
      Quadrature::provideQuadrature(vec.getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(<