DOFVector.cc 24.7 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
    switch (coarsenOperation) {
13
14
15
16
17
18
19
    case NO_OPERATION:
      return;
      break;
    case COARSE_RESTRICT:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
      break;
    case COARSE_INTERPOL:
20
      if (feSpace == NULL || !feSpace)
21
	ERROR_EXIT("ERR1\n");
22
23

      if (feSpace->getBasisFcts() == NULL || !(feSpace->getBasisFcts()))
24
	ERROR_EXIT("ERR2\n");
25

26
27
28
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
      break;
    default:
29
30
31
      std::cout << "ERROR WITH PTR = " << this << std::endl;
      ERROR_EXIT("Invalid coarsen operation \"%d\" in vector \"%s\"\n", 
		 coarsenOperation, this->name.c_str());    
32
33
34
35
36
37
38
39
40
41
42
43
    }
  }

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

  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
57
58
59
60
61
62
63
64
65
    vec[dof_new] = vec[dof0];
    vec[dof_new] += vec[dof1];
    vec[dof_new] *= 0.5;
  }

  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
66
    if (grad) {
67
68
69
      result = grad;
    } else {
      if(result && result->getFESpace() != feSpace) {
70
71
	delete result;
	result = new DOFVector<WorldVector<double> >(feSpace, "gradient");
72
73
74
75
76
77
78
79
80
      }
    }

    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
81
82
83
    std::vector<int> numNodeDOFs;
    std::vector<int> numNodePreDOFs;
    std::vector<DimVec<double>*> bary;
84
85
86
87

    int numNodes = 0;
    int numDOFs = 0;

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

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

    double *localUh = new double[basFcts->getNumber()];
108
109

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

137
138
    delete [] localUh;

139
140
141
142
143
144
145
146
147
    return result;
  }

  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
158
    if (!result) {
      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
183
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET | 
					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
184

185
186
    double *localUh = new double[basFcts->getNumber()];

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
	(*result)[dofIndex] += grd * det;
	volume[dofIndex] += det;
      }

      elInfo = stack.traverseNext(elInfo);
    }

203
204
    delete [] localUh;

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

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

    return result;
  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    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];
330
      for (int i = 0; i < nPoints; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
331
	grd[i].set(0.0);
332

Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
335
336
337
      result = grd;
    }
  
    double *localVec = localVectors[myRank];
    getLocalVector(largeElInfo->getElement(), localVec);
338
339
340
341

    const BasisFunction *basFcts = feSpace->getBasisFcts();    
    mtl::dense2D<double> m(nBasFcts, nBasFcts);
    smallElInfo->getSubElemCoordsMat(m, 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
377
378
379
	  result[i][l] += grdLambda[k][l] * grd1[k];
      }
    }    
  
    return const_cast<const WorldVector<double>*>(result);
  }

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

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

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

    Element *el = elInfo->getElement(); 

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

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

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

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

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

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

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

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

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

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

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

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

  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
476
    FUNCNAME("DOFVector<double>::interpol()");
477
478
479
480
481
482
483
484
485
486

    const FiniteElemSpace *sourceFeSpace = source->getFESpace();
    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
487
    DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
488
    double *sourceLocalCoeffs = new double[nSourceBasisFcts];
489

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

497
      while (elInfo) {
498
499
	Element *el = elInfo->getElement();

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

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

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

557
558
559
560
561
562
563
564
565
566
567

  template<>
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, double factor) 
  {
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

    DimVec<double> *coords = NULL;
    const FiniteElemSpace *vFESpace = v->getFESpace();

568
    if (feSpace == vFESpace)
569
570
571
572
573
574
575
576
577
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    const BasisFunction *vBasFcts = vFESpace->getBasisFcts();

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
589
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
590
591
592
593

	v->getLocalVector(el, vLocalCoeffs);

	for (int i = 0; i < numBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
594
	  if (vec[myLocalIndices[i]] == nul) {
595
	    coords = basFcts->getCoords(i);
596
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs, NULL) * factor;
597
598
599
600
601
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }

602
      delete [] vLocalCoeffs;
603
604
605
606
607
608
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


609
610
611
612
613
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
    FUNCNAME("DOFVector<double>::getGradient()");

614
615
616
617
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif

618
619
620
621
622
623
624
625
626
627
628
    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;

629
    if (grad) {
630
631
      result = grad;
    } else {
632
      if (!result) {
633
	result = new WorldVector<DOFVector<double>*>;
634

635
636
	result->set(NULL);
      }
637
      for (int i = 0; i < dow; i++) {
638
	if ((*result)[i] && (*result)[i]->getFESpace() != feSpace) {
639
640
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
641
642
643
644
645
	}
      }
    }

    // count number of nodes and dofs per node
646
647
648
    std::vector<int> numNodeDOFs;
    std::vector<int> numNodePreDOFs;
    std::vector<DimVec<double>*> bary;
649
650
651
652

    int numNodes = 0;
    int numDOFs = 0;

653
    for (int i = 0; i < dim + 1; i++) {
654
655
656
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
      int numPreDOFs = admin->getNumberOfPreDOFs(i);
657
      for (int j = 0; j < numPositionNodes; j++) {
658
659
660
661
662
663
664
665
	int dofs = basFcts->getNumberOfDOFs(geoIndex);
	numNodeDOFs.push_back(dofs);
	numDOFs += dofs;
	numNodePreDOFs.push_back(numPreDOFs);
      }
      numNodes += numPositionNodes;
    }

666
    TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
667
      ("number of dofs != number of basis functions\n");
668
    
669
670
    for (int i = 0; i < numDOFs; i++)
      bary.push_back(basFcts->getCoords(i));    
671
672

    // traverse mesh
673
    std::vector<bool> visited(getUsedSize(), false);
674
    TraverseStack stack;
675
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
676
677
678
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;

679
    while (elInfo) {
680
681
682
683
684
      const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
      const double *localUh = getLocalVector(elInfo->getElement(), NULL);
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
685
686
      for (int i = 0; i < numNodes; i++) { // for all nodes
	for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
687
	  DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
688

689
	  if (!visited[dofIndex]) {
690
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
691

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


708
709
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
710
  {}
711

712
713
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
  {
714
    /*
715
716
    std::vector<DegreeOfFreedom>::iterator it;
    std::vector<DegreeOfFreedom>::iterator end = vec.end();
717
    DegreeOfFreedom pos = 0;
718
    for (it = vec.begin(); it != end; ++it, ++pos)
719
      if (*it == dof) *it = pos;
720
    */
Thomas Witkowski's avatar
Thomas Witkowski committed
721
  }
722
723
724
725

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

728
    TEST_EXIT_DBG(vec)("no vector\n");
729

730
    int dow = Global::getGeo(WORLD);
731
732
    static WorldVector<DOFVector<double>*> *result = NULL;

733
    if (!res && !result) {
734
735
736
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
	(*result)[i] = new DOFVector<double>(vec->getFESpace(), "transform");
737
738
739
740
741
    }

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

    int vecSize = vec->getSize();
742
743
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
744
745
746
747
748
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
  template<>
  const double *DOFVectorBase<double>::getVecAtQPs(const ElInfo *smallElInfo, 
						   const ElInfo *largeElInfo,
						   const Quadrature *quad,
						   const FastQuadrature *quadFast,
						   double *vecAtQPs) const
  {
    FUNCNAME("DOFVector<double>::getVecAtQPs()");
  
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");

    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
    }

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

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); 
    static double *localvec = NULL;
    double *result;

    if (vecAtQPs) {
      result = vecAtQPs;
    } else {
      if (localvec) 
	delete [] localvec;
      localvec = new double[nPoints];
779
780
      for (int i = 0; i < nPoints; i++)
	localvec[i] = 0.0;      
Thomas Witkowski's avatar
Thomas Witkowski committed
781
782
783
784
785
786
      result = localvec;
    }
      
    double *localVec = localVectors[omp_get_thread_num()];
    getLocalVector(largeElInfo->getElement(), localVec);

787
788
    mtl::dense2D<double> m(nBasFcts, nBasFcts);    
    smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
789
790
791
792
793

    for (int i = 0; i < nPoints; i++) {
      result[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
794
	for (int k = 0; k < nBasFcts; k++)
795
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
796
	
Thomas Witkowski's avatar
Thomas Witkowski committed
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
	result[i] += localVec[j] * val;
      }
    }

    return const_cast<const double*>(result);
  }

  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    FUNCNAME("DOFVector::assemble()");

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

814
    set_to_zero(this->elementVector);
Thomas Witkowski's avatar
Thomas Witkowski committed
815
816
817
818
819
820
821
822
823

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	   it != this->operators.end(); 
824
825
	   ++it, ++factorIt)
	if ((*it)->getNeedDualTraverse() == false)
Thomas Witkowski's avatar
Thomas Witkowski committed
826
	  (*it)->getElementVector(elInfo, this->elementVector, 
827
				  *factorIt ? **factorIt : 1.0);      
Thomas Witkowski's avatar
Thomas Witkowski committed
828
829
    }

830
    addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
831
832
833
834
835
836
837
838
839
840
841
842
843
844
  }

  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;

845
    set_to_zero(this->elementVector);
Thomas Witkowski's avatar
Thomas Witkowski committed
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865

    if (op) {
      ERROR_EXIT("TODO");
      //      op->getElementVector(mainElInfo, this->elementVector);
    } else {
      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);
	}
      }
    }
  

866
    addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
867
868
869
  }


870
871
}