DOFVector.cc 24.5 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
    double *sourceLocalCoeffs = new double[nSourceBasisFcts];
488

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
	    }	
	  }
	}
      
	nextTraverse = dualStack.traverseNext(&elInfo1, &elInfo2,
					      &elInfoSmall, &elInfoLarge);
      }
    }
  
555
    delete [] sourceLocalCoeffs;
556
557
  }

558
559
560
561
562
563
564
565
566
567
568

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

569
    if (feSpace == vFESpace)
570
571
572
573
574
575
576
577
578
      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
579
      DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
580
      WorldVector<double> *vLocalCoeffs = new WorldVector<double>[vNumBasFcts];
581
582
583
584
585
586
587
588
589
      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
590
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
591
592
593
594

	v->getLocalVector(el, vLocalCoeffs);

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

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


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

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

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

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

    int numNodes = 0;
    int numDOFs = 0;

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

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

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

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

680
    while (elInfo) {
681
682
683
684
685
686

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


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

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

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

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

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

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

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

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

    return r;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
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
779
780
781
782
783
784
785
786
787
788
  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);

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

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

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

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

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

  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;

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

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

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


872
873
}