DOFVector.cc 24.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
    case NO_OPERATION:
      return;
      break;
    case COARSE_RESTRICT:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
      break;
    case COARSE_INTERPOL:
19
20
21
22
23
24
      if (feSpace == NULL || !feSpace) {
	ERROR_EXIT("ERR1\n");
      }
      if (feSpace->getBasisFcts() == NULL || !(feSpace->getBasisFcts())) {
	ERROR_EXIT("ERR2\n");
      }
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
      (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)
  {
41
42
43
    if (n < 1) 
      return;

44
45
46
47
    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
48
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(feSpace->getMesh()->getDim(), n0);
49
50
51
52
53
54
55
56
57
58
59
60
61
62
    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
63
    if (grad) {
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
      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
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;
    }

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

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

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

    TraverseStack stack;

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

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

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

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

	    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
157
158
    if (!result) {
      if (vec && vec->getFESpace() != feSpace) {
159
160
161
	DELETE vec;
	vec = NULL;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
162
      if (!vec) {
163
164
165
166
167
168
169
170
171
172
173
	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
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

Thomas Witkowski's avatar
Thomas Witkowski committed
185
    while (elInfo) {
186
187
188
189
      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
190
      basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
191

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

    return result;
  }

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

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

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

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

    if (grdAtQPs) {
      result = grdAtQPs;
    } else {
240
241
242
      if (grd) 
	delete [] grd;

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

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

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

263
264
	for (int j = 0; j < nBasFcts; j++) {
	  for (int k = 0; k < parts; k++) {
265
266
267
268
	    grd1[k] += quadFast->getGradient(i, j, k) * localVec[j];
	  }
	}

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

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

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

292
293
294
	for (int l = 0; l < dow; l++) {
	  result[i][l] = 0.0;
	  for (int k = 0; k < parts; k++) {
295
296
297
298
299
300
301
302
303
304
	    result[i][l] += grdLambda[k][l] * grd1[k];
	  }
	}
      }
    }
  
    return const_cast<const WorldVector<double>*>(result);
  }

  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
  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");

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

    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];
      for (int i = 0; i < nPoints; i++) {
	grd[i].set(0.0);
      }
      result = grd;
    }
  
    double *localVec = localVectors[myRank];
    getLocalVector(largeElInfo->getElement(), localVec);
    DimMat<double> *m = smallElInfo->getSubElemCoordsMat();

    DimVec<double> &grd1 = *grdTmp[myRank];
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();

    
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    DimVec<double>* grdPhi = grdPhis[myRank];
    DimVec<double> tmp(dim, DEFAULT_VALUE, 0.0);
    
    for (int i = 0; i < nPoints; i++) {
      for (int j = 0; j < dim + 1; j++) {
	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);
	  tmp *= (*m)[j][k];
	  *grdPhi += tmp;
	}

	for (int k = 0; k < parts; k++) {
	  grd1[k] += (*grdPhi)[k] * localVec[j];
	}
      }
      
      for (int l = 0; l < dow; l++) {
	result[i][l] = 0.0;
	for (int k = 0; k < parts; k++) {
	  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,
388
389
390
391
392
							       const FastQuadrature *quadFast,
							       WorldMatrix<double>  *d2AtQPs) const
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
393
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
394

395
396
397
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
398
399
    }

400
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
401
402
403
404
      ("invalid basis functions");

    Element *el = elInfo->getElement(); 

Thomas Witkowski's avatar
Thomas Witkowski committed
405
    int myRank = omp_get_thread_num();
406
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
407
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
408
409
410
411
412
413
414
415

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

    if (d2AtQPs) {
      result = d2AtQPs;
    } else {
      if(vec) delete [] vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
      vec = new WorldMatrix<double>[nPoints];
      for (int i = 0; i < nPoints; i++) {
418
419
420
421
422
	vec[i].set(0.0);
      }
      result = vec;
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
423
    double *localVec = localVectors[myRank];
424
    getLocalVector(el, localVec);
425
426
427
428
429

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

430
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
433
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
434
435
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438
439
	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);
440
441
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
442
443
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
444
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
445
446
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
447
448
449
450
		result[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
451
      const BasisFunction *basFcts = feSpace->getBasisFcts();
452
453
      DimMat<double>* D2Phi = D2Phis[omp_get_thread_num()];

Thomas Witkowski's avatar
Thomas Witkowski committed
454
455
456
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
457
458
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
465
	      D2Tmp[k][l] += localVec[i] * (*D2Phi)[k][l];
466
467
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
470
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
471
472
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
473
		result[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
474
475
476
	  }
      }
    }
477

478
479
480
481
482
483
    return const_cast<const WorldMatrix<double>*>(result);
  }

  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
484
    FUNCNAME("DOFVector<double>::interpol()");
485
486
487
488
489
490
491
492
493
494

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

498
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
499
500
501
502
503
504
      DimVec<double> *coords = NULL;
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
508
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
509
510
	source->getLocalVector(el, sourceLocalCoeffs);

511
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
512
	  if (vec[myLocalIndices[i]] == 0.0) {
513
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
514
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
	  }
	}
	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);
533
      while (nextTraverse) {     
534
535
	basisFcts->getLocalIndices(elInfo1->getElement(), 
				   feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
536
				   myLocalIndices);
537
538
539
	source->getLocalVector(elInfo2->getElement(), 
			       sourceLocalCoeffs);

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

567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589

  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
590
      DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
591
592
593
594
595
596
597
598
599
600
      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
601
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
602
603
604
605

	v->getLocalVector(el, vLocalCoeffs);

	for (int i = 0; i < numBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
606
	  if (vec[myLocalIndices[i]] == nul) {
607
	    coords = basFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
608
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs,NULL) * factor;
609
610
611
612
613
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }

614
      DELETE [] vLocalCoeffs;
615
616
617
618
619
620
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
  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;

637
    if (grad) {
638
639
      result = grad;
    } else {
640
      if (!result) {
641
	result = NEW WorldVector<DOFVector<double>*>;
642

643
644
	result->set(NULL);
      }
645
      for (int i = 0; i < dow; i++) {
646
	if ((*result)[i] && (*result)[i]->getFESpace() != feSpace) {
647
648
649
650
651
652
653
	  DELETE (*result)[i];
	  (*result)[i] = NEW DOFVector<double>(feSpace, "gradient");
	}
      }
    }

    // count number of nodes and dofs per node
654
655
656
    std::vector<int> numNodeDOFs;
    std::vector<int> numNodePreDOFs;
    std::vector<DimVec<double>*> bary;
657
658
659
660

    int numNodes = 0;
    int numDOFs = 0;

661
    for (int i = 0; i < dim + 1; i++) {
662
663
664
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
      int numPreDOFs = admin->getNumberOfPreDOFs(i);
665
      for (int j = 0; j < numPositionNodes; j++) {
666
667
668
669
670
671
672
673
	int dofs = basFcts->getNumberOfDOFs(geoIndex);
	numNodeDOFs.push_back(dofs);
	numDOFs += dofs;
	numNodePreDOFs.push_back(numPreDOFs);
      }
      numNodes += numPositionNodes;
    }

674
675
676
//     TEST_EXIT_DBG(numNodes == mesh->getNumberOfNodes())
//       ("invalid number of nodes\n");

677
    TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
678
      ("number of dofs != number of basis functions\n");
679
    
680
    for (int i = 0; i < numDOFs; i++) {
681
682
683
684
      bary.push_back(basFcts->getCoords(i));
    }

    // traverse mesh
685
    std::vector<bool> visited(getUsedSize(), false);
686
    TraverseStack stack;
687
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
688
689
690
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;

691
    while (elInfo) {
692
693
694
695
696
697

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

      int localDOFNr = 0;
698
699
      for (int i = 0; i < numNodes; i++) { // for all nodes
	for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
700
	  DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
701

702
	  if (!visited[dofIndex]) {
703
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
704

705
	    for (int k = 0; k < dow; k++) {
706
	      (*(*result)[k])[dofIndex] = grd[k];
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
	    }

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


722
723
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
724
  {}
725

726

727
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) {
728
729
    std::vector<DegreeOfFreedom>::iterator it;
    std::vector<DegreeOfFreedom>::iterator end = vec.end();
730
731
732
    DegreeOfFreedom pos = 0;
    for (it = vec.begin(); it != end; ++it, ++pos) {
      if (*it == dof) *it = pos;
733
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
734
  }
735

736

737
738
739
  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *res)
  {
740
    FUNCNAME("DOFVector<double>::transform()");
741

742
    TEST_EXIT_DBG(vec)("no vector\n");
743

744
    int dow = Global::getGeo(WORLD);
745
746
    static WorldVector<DOFVector<double>*> *result = NULL;

747
    if (!res && !result) {
748
      result = NEW WorldVector<DOFVector<double>*>;
749
      for (int i = 0; i < dow; i++) {
750
751
752
753
754
755
756
	(*result)[i] = NEW DOFVector<double>(vec->getFESpace(), "transform");
      }
    }

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

    int vecSize = vec->getSize();
757
758
    for (int i = 0; i < vecSize; i++) {
      for (int j = 0; j < dow; j++) {
759
760
761
762
763
764
765
	(*((*r)[j]))[i] = (*vec)[i][j];
      }
    }

    return r;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
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
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
  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);

    DimMat<double> *m = smallElInfo->getSubElemCoordsMat();

    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++) {
	  val += (*m)[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
	}
	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()");

828
829
    TEST_EXIT_DBG(this->elementVector)("No ElementVector!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
830
831
832
833
834
    if (!(op || this->operators.size())) 
      return;

    Operator *operat = op ? op : this->operators[0];

835
    operat->getAssembler(omp_get_thread_num())->initElementVector(this->elementVector, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869

    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(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse() == false) {
	  (*it)->getElementVector(elInfo, this->elementVector, 
				  *factorIt ? **factorIt : 1.0);
	}
      }
    }

    addElementVector(factor, *this->elementVector, bound);
  }

  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;

    Operator *operat = op ? op : this->operators[0];

870
871
    operat->getAssembler(omp_get_thread_num())->
	initElementVector(this->elementVector, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895

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

    addElementVector(factor, *this->elementVector, bound);
  }


896
897
}