DOFVector.cc 24.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
    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
29
30
31
32
33
34
35
36
37
38
39
40
41
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
      break;
    default:
      ERROR_EXIT("invalid coarsen operation\n");
    }
  }

  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)
  {
42
43
44
    if (n < 1) 
      return;

45
46
47
48
    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
49
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(feSpace->getMesh()->getDim(), n0);
50
51
52
53
54
55
56
57
58
59
60
61
62
63
    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
64
    if (grad) {
65
66
67
      result = grad;
    } else {
      if(result && result->getFESpace() != feSpace) {
68
69
	delete result;
	result = new DOFVector<WorldVector<double> >(feSpace, "gradient");
70
71
72
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
88

    int numNodes = 0;
    int numDOFs = 0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
102
103
    //    TEST_EXIT_DBG(numNodes == mesh->getNumberOfNodes())
    //      ("invalid number of nodes\n");
104

105
    TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
106
      ("number of dofs != number of basis functions\n");
107
    
Thomas Witkowski's avatar
Thomas Witkowski committed
108
    for (int i = 0; i < numDOFs; i++) {
109
110
111
112
      bary.push_back(basFcts->getCoords(i));
    }

    // traverse mesh
113
    std::vector<bool> visited(getUsedSize(), false);
114
115
116

    TraverseStack stack;

Thomas Witkowski's avatar
Thomas Witkowski committed
117
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
118
119
120

    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

Thomas Witkowski's avatar
Thomas Witkowski committed
121
    while (elInfo) {
122
123
124
125
126
127

      const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
      const double *localUh = getLocalVector(elInfo->getElement(), NULL);
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
128
129
      for (int i = 0; i < numNodes; i++) { // for all nodes
	for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
130
	  DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
Thomas Witkowski's avatar
Thomas Witkowski committed
131
132
133
	  if (!visited[dofIndex]) {	  
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, 
			       localUh, &((*result)[dofIndex]));
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }

  template<>
  DOFVector<WorldVector<double> >*
  DOFVector<double>::getRecoveryGradient(DOFVector<WorldVector<double> > *grad) const 
  {
    FUNCNAME("DOFVector<double>::getRecoveryGradient()");

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
158
159
    if (!result) {
      if (vec && vec->getFESpace() != feSpace) {
160
	delete vec;
161
162
	vec = NULL;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
163
      if (!vec) {
164
	vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
165
166
167
168
169
170
171
172
173
174
      }
      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
175
    int nPreDOFs = feSpace->getAdmin()->getNumberOfPreDOFs(0);
176
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
Thomas Witkowski's avatar
Thomas Witkowski committed
177
    WorldVector<double> grd;
178
179

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
193
      for (int i = 0; i < dim + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
194
	DegreeOfFreedom dofIndex = dof[i][nPreDOFs];
195
196
197
198
199
200
201
202
203
204
	(*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);

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

    return result;
  }

  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
245
      grd = new WorldVector<double>[nPoints];
      for (int i = 0; i < nPoints; i++) {
246
247
248
249
250
	grd[i].set(0.0);
      }
      result = grd;
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
251
252
    double *localVec = localVectors[myRank];
    getLocalVector(elInfo->getElement(), 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
299
    return const_cast<const WorldVector<double>*>(result);
  }

  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
302
303
304
305
306
307
308
309
  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");

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

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

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

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

    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++) {
350
      for (int j = 0; j < dim + 1; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353
354
355
356
357
	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);
358
	  tmp *= m[j][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
361
	  *grdPhi += tmp;
	}

362
	for (int k = 0; k < parts; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
363
364
365
366
367
	  grd1[k] += (*grdPhi)[k] * localVec[j];
      }
      
      for (int l = 0; l < dow; l++) {
	result[i][l] = 0.0;
368
	for (int k = 0; k < parts; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
371
372
373
374
375
376
377
378
	  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,
379
380
381
382
383
							       const FastQuadrature *quadFast,
							       WorldMatrix<double>  *d2AtQPs) const
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
384
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
385

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

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

    Element *el = elInfo->getElement(); 

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

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

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

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

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

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

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

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

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

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

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

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

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

    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
486
    DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
487
488
    double *sourceLocalCoeffs = GET_MEMORY(double, nSourceBasisFcts);

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

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

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

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

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

558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580

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

    if (feSpace == vFESpace) {
      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
581
      DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
582
      WorldVector<double> *vLocalCoeffs = new WorldVector<double>[vNumBasFcts];
583
584
585
586
587
588
589
590
591
      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
592
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
593
594
595
596

	v->getLocalVector(el, vLocalCoeffs);

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

605
      delete [] vLocalCoeffs;
606
607
608
609
610
611
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
    FUNCNAME("DOFVector<double>::getGradient()");

    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;

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

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

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

    int numNodes = 0;
    int numDOFs = 0;

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

665
666
667
//     TEST_EXIT_DBG(numNodes == mesh->getNumberOfNodes())
//       ("invalid number of nodes\n");

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

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

682
    while (elInfo) {
683
684
685
686
687
688

      const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
      const double *localUh = getLocalVector(elInfo->getElement(), NULL);
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
689
690
      for (int i = 0; i < numNodes; i++) { // for all nodes
	for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
691
	  DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
692

693
	  if (!visited[dofIndex]) {
694
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
695

696
	    for (int k = 0; k < dow; k++) {
697
	      (*(*result)[k])[dofIndex] = grd[k];
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
	    }

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


713
714
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
715
  {}
716

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

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

731
    TEST_EXIT_DBG(vec)("no vector\n");
732

733
    int dow = Global::getGeo(WORLD);
734
735
    static WorldVector<DOFVector<double>*> *result = NULL;

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

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

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

    return r;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
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
779
780
781
782
783
784
785
786
787
788
789
790
  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];
      for (int i = 0; i < nPoints; i++) {
	localvec[i] = 0.0;
      }
      result = localvec;
    }
      
    double *localVec = localVectors[omp_get_thread_num()];
    getLocalVector(largeElInfo->getElement(), localVec);

791
792
    mtl::dense2D<double> m(nBasFcts, nBasFcts);    
    smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
793
794
795
796
797
798

    for (int i = 0; i < nPoints; i++) {
      result[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
	for (int k = 0; k < nBasFcts; k++) {
799
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
Thomas Witkowski's avatar
Thomas Witkowski committed
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
	}
	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;

818
    set_to_zero(this->elementVector);
Thomas Witkowski's avatar
Thomas Witkowski committed
819
820
821
822
823
824
825
826
827

    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(); 
828
829
	   ++it, ++factorIt)
	if ((*it)->getNeedDualTraverse() == false)
Thomas Witkowski's avatar
Thomas Witkowski committed
830
	  (*it)->getElementVector(elInfo, this->elementVector, 
831
				  *factorIt ? **factorIt : 1.0);      
Thomas Witkowski's avatar
Thomas Witkowski committed
832
833
    }

834
    addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
835
836
837
838
839
840
841
842
843
844
845
846
847
848
  }

  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;

849
    set_to_zero(this->elementVector);
Thomas Witkowski's avatar
Thomas Witkowski committed
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869

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

870
    addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
871
872
873
  }


874
875
}