DOFVector.cc 24.4 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
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
      break;
    default:
29
30
      WARNING("Invalid coarsen operation \"%d\" in vector \"%s\"\n", 
	      coarsenOperation, this->name.c_str());
31
32
33
34
35
36
37
38
39
40
41
42
    }
  }

  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)
  {
43
44
45
    if (n < 1) 
      return;

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

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

84
85
    int nNodes = 0;
    int nDofs = 0;
86

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

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

136
137
    delete [] localUh;

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
186
    while (elInfo) {
187
188
189
      double det = elInfo->getDet();
      const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
190
      getLocalVector(elInfo->getElement(), localUh);
Thomas Witkowski's avatar
Thomas Witkowski committed
191
      basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
192

Thomas Witkowski's avatar
Thomas Witkowski committed
193
      for (int i = 0; i < dim + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
194
	DegreeOfFreedom dofIndex = dof[i][nPreDOFs];
195
196
197
198
199
200
201
	(*result)[dofIndex] += grd * det;
	volume[dofIndex] += det;
      }

      elInfo = stack.traverseNext(elInfo);
    }

202
203
    delete [] localUh;

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

207
208
    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
      if (*volIt != 0.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
209
	*grdIt *= 1.0 / (*volIt);
210
211
212
213
214

    return result;
  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
302
303
304
305
306
307
308
309
  const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *smallElInfo,
								const ElInfo *largeElInfo,
								const Quadrature *quad,
								const FastQuadrature *quadFast,
								WorldVector<double> *grdAtQPs) const
  {
    FUNCNAME("DOFVector<double>::getGrdAtQPs()");
  
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");

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

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

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

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

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

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

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

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

    Element *el = elInfo->getElement(); 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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
	      if (coords2[j] < -0.00001) {
538
539
540
541
542
		isPositive = false;
		break;
	      }
	    }
	  
543
544
545
	    if (isPositive)
	      vec[myLocalIndices[i]] = 
		sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	   
546
547
548
	  }
	}
      
549
550
	nextTraverse = 
	  dualStack.traverseNext(&elInfo1, &elInfo2, &elInfoSmall, &elInfoLarge);
551
552
553
      }
    }
  
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
    int nNodes = 0;
    int nDofs = 0;
652

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
	int dofs = basFcts->getNumberOfDOFs(geoIndex);
	numNodeDOFs.push_back(dofs);
660
	nDofs += dofs;
661
662
	numNodePreDOFs.push_back(numPreDOFs);
      }
663
      nNodes += numPositionNodes;
664
665
    }

666
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
667
      ("number of dofs != number of basis functions\n");
668
    
669
    for (int i = 0; i < nDofs; i++)
670
      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
      const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
      const double *localUh = getLocalVector(elInfo->getElement(), NULL);
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


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

712

713
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
714
715
  {}

716
717
718
719

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

722
    TEST_EXIT_DBG(vec)("no vector\n");
723

724
    int dow = Global::getGeo(WORLD);
725
726
    static WorldVector<DOFVector<double>*> *result = NULL;

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

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

    int vecSize = vec->getSize();
736
737
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
738
739
740
741
742
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

743

Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
746
747
748
749
750
751
752
753
754
  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");

755
    if (quad && quadFast)
Thomas Witkowski's avatar
Thomas Witkowski committed
756
757
758
759
760
761
762
      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();
763
764
    int nPoints = 
      quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
765
766
767
768
769
770
771
772
773
    static double *localvec = NULL;
    double *result;

    if (vecAtQPs) {
      result = vecAtQPs;
    } else {
      if (localvec) 
	delete [] localvec;
      localvec = new double[nPoints];
774
775
      for (int i = 0; i < nPoints; i++)
	localvec[i] = 0.0;      
Thomas Witkowski's avatar
Thomas Witkowski committed
776
777
778
779
780
781
      result = localvec;
    }
      
    double *localVec = localVectors[omp_get_thread_num()];
    getLocalVector(largeElInfo->getElement(), localVec);

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

    for (int i = 0; i < nPoints; i++) {
      result[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
789
	for (int k = 0; k < nBasFcts; k++)
790
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
791
	
Thomas Witkowski's avatar
Thomas Witkowski committed
792
793
794
795
796
797
798
	result[i] += localVec[j] * val;
      }
    }

    return const_cast<const double*>(result);
  }

799

Thomas Witkowski's avatar
Thomas Witkowski committed
800
801
802
803
804
805
806
807
808
809
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    FUNCNAME("DOFVector::assemble()");

    if (!(op || this->operators.size())) 
      return;

810
    set_to_zero(this->elementVector);
811
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
812
813
814

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
815
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
816
817
818
819
820
821
    } 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
	   ++it, ++factorIt)
823
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
824
	  (*it)->getElementVector(elInfo, this->elementVector, 
825
				  *factorIt ? **factorIt : 1.0);      
826
827
	  addVector = true;
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
828
829
    }

830
831
    if (addVector)
      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);
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
    bool addVector = false;

    TEST_EXIT(!op)("Not yet implemented!\n");

    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);
	addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
862
863
      }
  
864
865
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
866
867
868
  }


869
870
}