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
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
525
				   myLocalIndices);
526
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
527

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

555
556
557
558
559
560
561
562
563
564
565

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

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

	v->getLocalVector(el, vLocalCoeffs);

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

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


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

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

616
617
618
619
620
621
622
623
624
625
626
    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;

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

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

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

    int numNodes = 0;
    int numDOFs = 0;

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

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

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

677
    while (elInfo) {
678
679
680
681
682
      const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
      const double *localUh = getLocalVector(elInfo->getElement(), NULL);
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

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

687
	  if (!visited[dofIndex]) {
688
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
689

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


706
707
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
708
  {}
709

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

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

726
    TEST_EXIT_DBG(vec)("no vector\n");
727

728
    int dow = Global::getGeo(WORLD);
729
730
    static WorldVector<DOFVector<double>*> *result = NULL;

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

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

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

    return r;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
749
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
  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];
777
778
      for (int i = 0; i < nPoints; i++)
	localvec[i] = 0.0;      
Thomas Witkowski's avatar
Thomas Witkowski committed
779
780
781
782
783
784
      result = localvec;
    }
      
    double *localVec = localVectors[omp_get_thread_num()];
    getLocalVector(largeElInfo->getElement(), localVec);

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

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

812
    set_to_zero(this->elementVector);
Thomas Witkowski's avatar
Thomas Witkowski committed
813
814
815
816
817
818
819
820
821

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

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

  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;

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

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

864
    addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
865
866
867
  }


868
869
}