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

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

    int numNodes = 0;
    int numDOFs = 0;

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

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

    double *localUh = new double[basFcts->getNumber()];
106
107

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

135
136
    delete [] localUh;

137
138
139
140
141
142
143
144
145
    return result;
  }

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

146
147
148
149
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif

150
151
152
153
154
    // define result vector
    static DOFVector<WorldVector<double> > *vec = NULL;

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

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

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

183
184
    double *localUh = new double[basFcts->getNumber()];

Thomas Witkowski's avatar
Thomas Witkowski committed
185
    while (elInfo) {
186
187
188
      double det = elInfo->getDet();
      const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
189
      getLocalVector(elInfo->getElement(), localUh);
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
	(*result)[dofIndex] += grd * det;
	volume[dofIndex] += det;
      }

      elInfo = stack.traverseNext(elInfo);
    }

201
202
    delete [] localUh;

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

206
207
    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
      if (*volIt != 0.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
208
	*grdIt *= 1.0 / (*volIt);
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
      grd = new WorldVector<double>[nPoints];
244
      for (int i = 0; i < nPoints; i++)
245
	grd[i].set(0.0);
246
      
247
248
249
      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
261
	  grd1[j] = 0.0;

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

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

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

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

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

295
296
297
298
    return const_cast<const WorldVector<double>*>(result);
  }

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

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

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

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

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

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

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

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

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

    Element *el = elInfo->getElement(); 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

557
558
559
560
561
562
563
564
565
566
567

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

    this->set(nul);

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

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

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

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

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

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

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

	v->getLocalVector(el, vLocalCoeffs);

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

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


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

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

618
619
620
621
622
623
624
625
626
627
628
    Mesh *mesh = feSpace->getMesh();
    int dim = mesh->getDim();
    int dow = Global::getGeo(WORLD);

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

    DOFAdmin *admin = feSpace->getAdmin();

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

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

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

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

    int numNodes = 0;
    int numDOFs = 0;

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

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

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

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


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

713
714
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
  {
715
    /*
716
717
    std::vector<DegreeOfFreedom>::iterator it;
    std::vector<DegreeOfFreedom>::iterator end = vec.end();
718
    DegreeOfFreedom pos = 0;
719
    for (it = vec.begin(); it != end; ++it, ++pos)
720
      if (*it == dof) *it = pos;
721
    */
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
  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];
780
781
      for (int i = 0; i < nPoints; i++)
	localvec[i] = 0.0;      
Thomas Witkowski's avatar
Thomas Witkowski committed
782
783
784
785
786
787
      result = localvec;
    }
      
    double *localVec = localVectors[omp_get_thread_num()];
    getLocalVector(largeElInfo->getElement(), localVec);

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

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

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

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

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

  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;

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

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

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


871
872
}