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
49
50
    Element *el = list.getElement(0);
    int n0 = feSpace->getAdmin()->getNumberOfPreDOFs(VERTEX);
    DegreeOfFreedom dof0 = el->getDOF(0, n0);
    DegreeOfFreedom dof1 = el->getDOF(1, n0);
Thomas Witkowski's avatar
Thomas Witkowski committed
51
    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
118
      const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
      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
189
190
      double det = elInfo->getDet();
      const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
      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
755
756
757
758
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");

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

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

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

766
767
768
769
    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
770

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

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

784

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
819

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

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


855
856
857
858
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct)
  {
859
    if (vec1.getFeSpace()->getMesh() == vec2.getFeSpace()->getMesh())
860
861
862
863
864
865
866
867
868
869
870
871
872
873
      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");

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

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

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

      double tmp = 0.0;
893
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
894
895
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);     
      value += tmp * elInfo->getDet();
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
      
      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");

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

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

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

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

    return value;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
945
946
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
  
  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(vec.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<double> qp(fastQuad->getNumPoints());
    WorldVector<double> coords;

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
    ElInfo *elInfo = stack.traverseFirst(vec.getFeSpace()->getMesh(), -1, traverseFlag);
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);

      double tmp = 0.0;
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
 	tmp += fastQuad->getWeight(iq) * (qp[iq] * (*fct)(coords));
      }

      value += tmp * elInfo->getDet();
      
      elInfo = stack.traverseNext(elInfo);
    }

    return value;
  }
983
984
}