DOFVector.cc 24.7 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
86

    int numNodes = 0;
    int numDOFs = 0;

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

100
    TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
101
      ("number of dofs != number of basis functions\n");
102
    
103
104
105
106
    for (int i = 0; i < numDOFs; i++)
      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;
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
      for (int i = 0; i < numNodes; i++) { // for all nodes
	for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
122
	  DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + 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
311
    if (quad && quadFast)
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328

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

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

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

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

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

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

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

    Element *el = elInfo->getElement(); 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

556
557
558
559
560
561
562
563
564
565
566

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

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

	v->getLocalVector(el, vLocalCoeffs);

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

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


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

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

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

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

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

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

    int numNodes = 0;
    int numDOFs = 0;

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

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

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


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

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

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

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

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

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

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

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

    return r;
  }

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

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

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

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

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

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

  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;

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

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

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


869
870
}