DOFVector.cc 18.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#include "DOFVector.h"
#include "Traverse.h"
#include "DualTraverse.h"
#include "FixVec.h"

namespace AMDiS {

  template<>
  void DOFVector<double>::coarseRestrict(RCNeighbourList& list, int n)
  {
11
    switch (coarsenOperation) {
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
    case NO_OPERATION:
      return;
      break;
    case COARSE_RESTRICT:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
      break;
    case COARSE_INTERPOL:
      (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)
  {
35
36
37
    if (n < 1) 
      return;

38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
    int dim = feSpace->getMesh()->getDim();
    Element *el = list.getElement(0);
    int n0 = feSpace->getAdmin()->getNumberOfPreDOFs(VERTEX);
    DegreeOfFreedom dof0 = el->getDOF(0, n0);
    DegreeOfFreedom dof1 = el->getDOF(1, n0);
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(dim, n0);
    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
58
    if (grad) {
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
      result = grad;
    } else {
      if(result && result->getFESpace() != feSpace) {
	DELETE result;
	result = NEW DOFVector<WorldVector<double> >(feSpace, "gradient");
      }
    }

    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
76
77
78
    std::vector<int> numNodeDOFs;
    std::vector<int> numNodePreDOFs;
    std::vector<DimVec<double>*> bary;
79
80
81
82

    int numNodes = 0;
    int numDOFs = 0;

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

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

99
    TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
100
      ("number of dofs != number of basis functions\n");
101
    
Thomas Witkowski's avatar
Thomas Witkowski committed
102
    for (int i = 0; i < numDOFs; i++) {
103
104
105
106
      bary.push_back(basFcts->getCoords(i));
    }

    // traverse mesh
107
    std::vector<bool> visited(getUsedSize(), false);
108
109
110

    TraverseStack stack;

Thomas Witkowski's avatar
Thomas Witkowski committed
111
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
112
113
114

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

Thomas Witkowski's avatar
Thomas Witkowski committed
115
    while (elInfo) {
116
117
118
119
120
121

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

	    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
152
153
    if (!result) {
      if (vec && vec->getFESpace() != feSpace) {
154
155
156
	DELETE vec;
	vec = NULL;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
157
      if (!vec) {
158
159
160
161
162
163
164
165
166
167
168
	vec = NEW DOFVector<WorldVector<double> >(feSpace, "gradient");
      }
      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
169
    int nPreDOFs = feSpace->getAdmin()->getNumberOfPreDOFs(0);
170
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
Thomas Witkowski's avatar
Thomas Witkowski committed
171
    WorldVector<double> grd;
172
173

    // traverse mesh
Thomas Witkowski's avatar
Thomas Witkowski committed
174
    Mesh *mesh = feSpace->getMesh();
175
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
176
177
178
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET | 
					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
179

Thomas Witkowski's avatar
Thomas Witkowski committed
180
    while (elInfo) {
181
182
183
184
      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
185
      basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
186

Thomas Witkowski's avatar
Thomas Witkowski committed
187
      for (int i = 0; i < dim + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
188
	DegreeOfFreedom dofIndex = dof[i][nPreDOFs];
189
190
191
192
193
194
195
196
197
198
	(*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
199
200
    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt) {
      if (*volIt != 0.0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
201
	*grdIt *= 1.0 / (*volIt);
202
203
204
205
206
207
208
      }
    }

    return result;
  }

  template<>
209
210
  const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo, 
								const Quadrature *quad,
211
								const FastQuadrature *quadFast,
212
								WorldVector<double> *grdAtQPs) const
213
214
215
  {
    FUNCNAME("DOFVector<double>::getGrdAtQPs()");
  
216
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
217

218
219
220
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
221
222
    }

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

    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
    int myRank = omp_get_thread_num();
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
229
230
231
232
233
234
    static WorldVector<double> *grd = NULL;
    WorldVector<double> *result;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
238
239
      grd = new WorldVector<double>[nPoints];
      for (int i = 0; i < nPoints; i++) {
240
241
242
243
244
	grd[i].set(0.0);
      }
      result = grd;
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
    double *localVec = localVectors[myRank];
    getLocalVector(elInfo->getElement(), localVec);
247

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

252
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
253
      for (int i = 0; i < nPoints; i++) {
254
	for (int j = 0; j < dim + 1; j++) {
255
	  grd1[j] = 0.0;
256
	}
257

258
259
	for (int j = 0; j < nBasFcts; j++) {
	  for (int k = 0; k < parts; k++) {
260
261
262
263
	    grd1[k] += quadFast->getGradient(i, j, k) * localVec[j];
	  }
	}

264
265
266
	for (int l=0; l < dow; l++) {
	  result[i][l] = 0.0;
	  for (int k = 0; k < parts; k++) {
267
268
269
270
271
	    result[i][l] += grdLambda[k][l] * grd1[k];
	  }
	}
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
      const BasisFunction *basFcts = feSpace->getBasisFcts();
      DimVec<double>* grdPhi = grdPhis[myRank];
274

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

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

287
288
289
	for (int l = 0; l < dow; l++) {
	  result[i][l] = 0.0;
	  for (int k = 0; k < parts; k++) {
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
	    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,
							       const FastQuadrature *quadFast,
							       WorldMatrix<double>  *d2AtQPs) const
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
307
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
308

309
310
311
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
312
313
    }

314
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
315
316
317
318
      ("invalid basis functions");

    Element *el = elInfo->getElement(); 

Thomas Witkowski's avatar
Thomas Witkowski committed
319
    int myRank = omp_get_thread_num();
320
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
321
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
322
323
324
325
326
327
328
329

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

    if (d2AtQPs) {
      result = d2AtQPs;
    } else {
      if(vec) delete [] vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
      vec = new WorldMatrix<double>[nPoints];
      for (int i = 0; i < nPoints; i++) {
332
333
334
335
336
	vec[i].set(0.0);
      }
      result = vec;
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
337
    double *localVec = localVectors[myRank];
338
    getLocalVector(el, localVec);
339
340
341
342
343

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

344
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
347
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
348
349
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
350
351
352
353
	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);
354
355
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
358
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
361
362
363
364
		result[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
365
      const BasisFunction *basFcts = feSpace->getBasisFcts();
366
367
      DimMat<double>* D2Phi = D2Phis[omp_get_thread_num()];

Thomas Witkowski's avatar
Thomas Witkowski committed
368
369
370
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
371
372
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
379
	      D2Tmp[k][l] += localVec[i] * (*D2Phi)[k][l];
380
381
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
384
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
387
		result[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
388
389
390
	  }
      }
    }
391

392
393
394
395
396
397
    return const_cast<const WorldMatrix<double>*>(result);
  }

  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
398
    FUNCNAME("DOFVector<double>::interpol()");
399
400
401
402
403
404
405
406
407
408

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

412
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
413
414
415
416
417
418
      DimVec<double> *coords = NULL;
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
422
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
423
424
	source->getLocalVector(el, sourceLocalCoeffs);

425
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
426
	  if (vec[myLocalIndices[i]] == 0.0) {
427
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
428
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      DimVec<double> *coords1 = NULL;
      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);
448
      while (nextTraverse) {     
449
450
	basisFcts->getLocalIndices(elInfo1->getElement(), 
				   feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
451
				   myLocalIndices);
452
453
454
	source->getLocalVector(elInfo2->getElement(), 
			       sourceLocalCoeffs);

455
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
456
	  if (vec[myLocalIndices[i]] == 0.0) {
457
458
459
460
461
	    coords1 = basisFcts->getCoords(i);
	    elInfo1->coordToWorld(*coords1, &worldVec);
	    elInfo2->worldToCoord(worldVec, &coords2);
	  
	    bool isPositive = true;
462
	    for (int j = 0; j < coords2.size(); j++) {
463
464
465
466
467
468
	      if (coords2[j] < 0) {
		isPositive = false;
		break;
	      }
	    }
	  
469
	    if (isPositive) {
Thomas Witkowski's avatar
Thomas Witkowski committed
470
	      vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	    
471
472
473
474
475
476
477
478
479
480
481
482
	    }	
	  }
	}
      
	nextTraverse = dualStack.traverseNext(&elInfo1, &elInfo2,
					      &elInfoSmall, &elInfoLarge);
      }
    }
  
    FREE_MEMORY(sourceLocalCoeffs, double, nSourceBasisFcts);
  }

483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505

  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
506
      DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
507
508
509
510
511
512
513
514
515
516
      WorldVector<double> *vLocalCoeffs = NEW WorldVector<double>[vNumBasFcts];
      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
517
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
518
519
520
521

	v->getLocalVector(el, vLocalCoeffs);

	for (int i = 0; i < numBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
522
	  if (vec[myLocalIndices[i]] == nul) {
523
	    coords = basFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
524
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs,NULL) * factor;
525
526
527
528
529
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }

530
      DELETE [] vLocalCoeffs;
531
532
533
534
535
536
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
  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;

553
    if (grad) {
554
555
      result = grad;
    } else {
556
      if (!result) {
557
558
559
	result = NEW WorldVector<DOFVector<double>*>;
	result->set(NULL);
      }
560
      for (int i = 0; i < dow; i++) {
561
	if ((*result)[i] && (*result)[i]->getFESpace() != feSpace) {
562
563
564
565
566
567
568
	  DELETE (*result)[i];
	  (*result)[i] = NEW DOFVector<double>(feSpace, "gradient");
	}
      }
    }

    // count number of nodes and dofs per node
569
570
571
    std::vector<int> numNodeDOFs;
    std::vector<int> numNodePreDOFs;
    std::vector<DimVec<double>*> bary;
572
573
574
575

    int numNodes = 0;
    int numDOFs = 0;

576
    for (int i = 0; i < dim + 1; i++) {
577
578
579
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
      int numPreDOFs = admin->getNumberOfPreDOFs(i);
580
      for (int j = 0; j < numPositionNodes; j++) {
581
582
583
584
585
586
587
588
	int dofs = basFcts->getNumberOfDOFs(geoIndex);
	numNodeDOFs.push_back(dofs);
	numDOFs += dofs;
	numNodePreDOFs.push_back(numPreDOFs);
      }
      numNodes += numPositionNodes;
    }

589
    TEST_EXIT_DBG(numNodes == mesh->getNumberOfNodes())
590
      ("invalid number of nodes\n");
591
    TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
592
      ("number of dofs != number of basis functions\n");
593
    
594
    for (int i = 0; i < numDOFs; i++) {
595
596
597
598
      bary.push_back(basFcts->getCoords(i));
    }

    // traverse mesh
599
    std::vector<bool> visited(getUsedSize(), false);
600
601
602
603
604
    TraverseStack stack;
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;

605
    while (elInfo) {
606
607
608
609
610
611

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

      int localDOFNr = 0;
612
613
      for (int i = 0; i < numNodes; i++) { // for all nodes
	for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
614
	  DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
615

616
	  if (!visited[dofIndex]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
617
618
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, 
			       localUh, &grd);
619

620
	    for (int k = 0; k < dow; k++) {
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
	      (*result)[k][dofIndex] = grd[k];
	    }

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


637
638
639
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
  {};
640

641

642
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) {
643
644
    std::vector<DegreeOfFreedom>::iterator it;
    std::vector<DegreeOfFreedom>::iterator end = vec.end();
645
646
647
    DegreeOfFreedom pos = 0;
    for (it = vec.begin(); it != end; ++it, ++pos) {
      if (*it == dof) *it = pos;
648
649
650
    }
  };

651

652
653
654
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *res)
  {
655
    FUNCNAME("DOFVector<double>::transform()");
656

657
    TEST_EXIT_DBG(vec)("no vector\n");
658

659
    int dow = Global::getGeo(WORLD);
660
661
    static WorldVector<DOFVector<double>*> *result = NULL;

662
    if (!res && !result) {
663
      result = NEW WorldVector<DOFVector<double>*>;
664
      for (int i = 0; i < dow; i++) {
665
666
667
668
669
670
671
	(*result)[i] = NEW DOFVector<double>(vec->getFESpace(), "transform");
      }
    }

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

    int vecSize = vec->getSize();
672
673
    for (int i = 0; i < vecSize; i++) {
      for (int j = 0; j < dow; j++) {
674
675
676
677
678
679
680
681
682
	(*((*r)[j]))[i] = (*vec)[i][j];
      }
    }

    return r;
  }

}