DOFVector.cc 28.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
13
    FUNCNAME("DOFVector<double>::coarseRestrict()");

14
    switch (coarsenOperation) {
15
16
17
18
19
20
21
    case NO_OPERATION:
      return;
      break;
    case COARSE_RESTRICT:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
      break;
    case COARSE_INTERPOL:
22
23
      TEST_EXIT_DBG(feSpace)("Should not happen!\n");
      TEST_EXIT_DBG(feSpace->getBasisFcts())("Shoud not happen!\n");
24

25
26
27
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
      break;
    default:
28
29
      WARNING("Invalid coarsen operation \"%d\" in vector \"%s\"\n", 
	      coarsenOperation, this->name.c_str());
30
31
32
    }
  }

33

34
35
36
37
38
39
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList& list, int n)
  {
    (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->refineInter(this, &list, n);
  }

40

41
42
43
  template<>
  void DOFVector<WorldVector<double> >::refineInterpol(RCNeighbourList& list, int n)
  {
44
45
46
    if (n < 1) 
      return;

47
48
    Element *el = list.getElement(0);
    int n0 = feSpace->getAdmin()->getNumberOfPreDOFs(VERTEX);
49
50
51
    DegreeOfFreedom dof0 = el->getDof(0, n0);
    DegreeOfFreedom dof1 = el->getDof(1, n0);
    DegreeOfFreedom dof_new = el->getChild(0)->getDof(feSpace->getMesh()->getDim(), n0);
52
53
54
55
56
    vec[dof_new] = vec[dof0];
    vec[dof_new] += vec[dof1];
    vec[dof_new] *= 0.5;
  }

57

58
59
60
61
62
63
64
65
66
  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
67
    if (grad) {
68
69
      result = grad;
    } else {
70
      if(result && result->getFeSpace() != feSpace) {
71
72
	delete result;
	result = new DOFVector<WorldVector<double> >(feSpace, "gradient");
73
74
75
76
77
78
79
80
81
      }
    }

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

86
87
    int nNodes = 0;
    int nDofs = 0;
88

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

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

108
    ElementVector localUh(basFcts->getNumber());
109
110

    // traverse mesh
111
    std::vector<bool> visited(getUsedSize(), false);
112
    TraverseStack stack;
113
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
114
115
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }

141

142
143
144
145
146
147
  template<>
  DOFVector<WorldVector<double> >*
  DOFVector<double>::getRecoveryGradient(DOFVector<WorldVector<double> > *grad) const 
  {
    FUNCNAME("DOFVector<double>::getRecoveryGradient()");

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
157
    if (!result) {
158
      if (vec && vec->getFeSpace() != feSpace) {
159
	delete vec;
160
161
	vec = NULL;
      }
162
      if (!vec)
163
	vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
164
      
165
166
167
168
169
170
171
172
173
      result = vec;
    }

    result->set(WorldVector<double>(DEFAULT_VALUE, 0.0));

    DOFVector<double> volume(feSpace, "volume");
    volume.set(0.0);

    const BasisFunction *basFcts = feSpace->getBasisFcts();
Thomas Witkowski's avatar
Thomas Witkowski committed
174
    int nPreDOFs = feSpace->getAdmin()->getNumberOfPreDOFs(0);
175

176
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
Thomas Witkowski's avatar
Thomas Witkowski committed
177
    WorldVector<double> grd;
178
179

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

186
    ElementVector localUh(basFcts->getNumber());
187

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

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

    return result;
  }

214

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

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

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

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

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

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

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

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

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

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

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

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

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

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

300

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
337
      result = grd;
    }
  
338
339
    const ElementVector& localVec = localVectors[myRank];
    getLocalVector(largeElInfo->getElement(), const_cast<ElementVector&>(localVec));
340
341

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

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

364
	for (int k = 0; k < parts; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
367
368
369
	  grd1[k] += (*grdPhi)[k] * localVec[j];
      }
      
      for (int l = 0; l < dow; l++) {
	result[i][l] = 0.0;
370
	for (int k = 0; k < parts; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
373
374
375
376
377
	  result[i][l] += grdLambda[k][l] * grd1[k];
      }
    }    
  
    return const_cast<const WorldVector<double>*>(result);
  }

378

Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
381
  template<>
  const WorldMatrix<double> *DOFVectorBase<double>::getD2AtQPs(const ElInfo *elInfo, 
							       const Quadrature *quad,
382
383
384
385
386
							       const FastQuadrature *quadFast,
							       WorldMatrix<double>  *d2AtQPs) const
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
387
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
388

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

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

397
398
    Element *el = elInfo->getElement(); 

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

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

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

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

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

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

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

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

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

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

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

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

475

476
477
478
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
479
    FUNCNAME("DOFVector<double>::interpol()");
480

481
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
482
483
484
485
486
487
488
489
    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
490
    DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
491
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
492

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
503
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
504
505
	source->getLocalVector(el, sourceLocalCoeffs);

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

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

558
559

  template<>
560
561
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
562
563
564
565
566
567
  {
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

    DimVec<double> *coords = NULL;
568
    const FiniteElemSpace *vFeSpace = v->getFeSpace();
569

570
    if (feSpace == vFeSpace)
571
572
573
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
574
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
575
576
577
578

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

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

      while (elInfo) {
	Element *el = elInfo->getElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
589
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
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
601
602
603
604
605
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


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

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

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

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

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

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

647
648
    int nNodes = 0;
    int nDofs = 0;
649

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

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

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

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

      int localDOFNr = 0;
683
      for (int i = 0; i < nNodes; i++) { // for all nodes
684
	for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
685
	  DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[i] + 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
715
716
717

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

720
    TEST_EXIT_DBG(vec)("no vector\n");
721

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

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

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

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

    return r;
  }

741

Thomas Witkowski's avatar
Thomas Witkowski committed
742
  template<>
743
744
745
746
747
  void DOFVectorBase<double>::getVecAtQPs(const ElInfo *smallElInfo, 
					  const ElInfo *largeElInfo,
					  const Quadrature *quad,
					  const FastQuadrature *quadFast,
					  mtl::dense_vector<double>& vecAtQPs) const
Thomas Witkowski's avatar
Thomas Witkowski committed
748
749
  {
    FUNCNAME("DOFVector<double>::getVecAtQPs()");
750
 
Thomas Witkowski's avatar
Thomas Witkowski committed
751
752
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");

753
    if (quad && quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
754
755
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
756
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
757
758
759
760

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

761
    if (smallElInfo->getMesh() == feSpace->getMesh())
762
      return getVecAtQPs(smallElInfo, quad, quadFast, vecAtQPs);    
763

Thomas Witkowski's avatar
Thomas Witkowski committed
764
    const BasisFunction *basFcts = feSpace->getBasisFcts();
765
766
    int nPoints = 
      quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
767

768
769
770
771
    vecAtQPs.change_dim(nPoints);
    
    const ElementVector &localVec = localVectors[omp_get_thread_num()];
    getLocalVector(largeElInfo->getElement(), const_cast<ElementVector&>(localVec));
Thomas Witkowski's avatar
Thomas Witkowski committed
772

773
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
774
775

    for (int i = 0; i < nPoints; i++) {
776
      vecAtQPs[i] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
777
778
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
779
	for (int k = 0; k < nBasFcts; k++)
780
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
781
782

	vecAtQPs[i] += localVec[j] * val;
Thomas Witkowski's avatar
Thomas Witkowski committed
783
784
785
786
      }
    }
  }

787

Thomas Witkowski's avatar
Thomas Witkowski committed
788
789
790
791
792
793
794
795
796
797
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    FUNCNAME("DOFVector::assemble()");

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

798
    set_to_zero(this->elementVector);
799
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
800
801
802

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
803
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
804
805
806
807
808
809
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	   it != this->operators.end(); 
810
	   ++it, ++factorIt)
811
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
812
	  (*it)->getElementVector(elInfo, this->elementVector, 
813
				  *factorIt ? **factorIt : 1.0);      
814
815
	  addVector = true;
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
816
817
    }

818
819
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
820
821
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
822

Thomas Witkowski's avatar
Thomas Witkowski committed
823
824
825
826
827
828
829
830
831
832
833
834
  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;

835
    set_to_zero(this->elementVector);
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
    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
851
852
      }
  
853
854
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
855
856
857
  }


858
859
860
861
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct)
  {
862
    if (vec1.getFeSpace()->getMesh() == vec2.getFeSpace()->getMesh())
863
864
865
866
867
868
869
870
871
872
873
874
875
876
      return intSingleMesh(vec1, vec2, fct);
    else
      return intDualMesh(vec1, vec2, fct);
  }


  double intSingleMesh(const DOFVector<double> &vec1,
		       const DOFVector<double> &vec2,
		       BinaryAbstractFunction<double, double, double> *fct)
  {
    FUNCNAME("intDualmesh()");
    
    TEST_EXIT(fct)("No function defined!\n");

877
878
    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
879
    Quadrature* quad = 
880
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
881
    FastQuadrature *fastQuad =
882
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
883

884
885
    mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
886
887
888
889

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
890
    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
891
    while (elInfo) {
892
893
894
895
      vec1.getVecAtQPs(elInfo, quad, fastQuad, qp1);
      vec2.getVecAtQPs(elInfo, quad, fastQuad, qp2);

      double tmp = 0.0;
896
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
897
898
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);     
      value += tmp * elInfo->getDet();
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
      
      elInfo = stack.traverseNext(elInfo);
    }

    return value;
  }


  double intDualMesh(const DOFVector<double> &vec1,
		     const DOFVector<double> &vec2,
		     BinaryAbstractFunction<double, double, double> *fct)
  {
    FUNCNAME("intDualmesh()");
    
    TEST_EXIT(fct)("No function defined!\n");

915
916
    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
917
    Quadrature* quad = 
918
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
919
    FastQuadrature *fastQuad =
920
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
921

922
923
    mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
924
925
926
927
928

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    DualTraverse dualTraverse;    
    DualElInfo dualElInfo;
929
930
    bool cont = dualTraverse.traverseFirst(vec1.getFeSpace()->getMesh(),
					   vec2.getFeSpace()->getMesh(),
931
932
933
					   -1, -1, traverseFlag, traverseFlag,
					   dualElInfo);
    while (cont) {
934
935
936
937
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);

      double tmp = 0.0;
938
      for (int iq = 0; iq < quad->getNumPoints(); iq++)
939
940
 	tmp += quad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);      
      value += tmp * dualElInfo.smallElInfo->getDet();
941
942
943
944
945
946
947
      
      cont = dualTraverse.traverseNext(dualElInfo);
    }

    return value;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977