DOFVector.cc 29.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include <boost/numeric/mtl/mtl.hpp>
14
15
16
17
18
19
20
21
22
23
#include "DOFVector.h"
#include "Traverse.h"
#include "DualTraverse.h"
#include "FixVec.h"

namespace AMDiS {

  template<>
  void DOFVector<double>::coarseRestrict(RCNeighbourList& list, int n)
  {
24
25
    FUNCNAME("DOFVector<double>::coarseRestrict()");

26
    switch (coarsenOperation) {
27
28
29
30
31
32
33
    case NO_OPERATION:
      return;
      break;
    case COARSE_RESTRICT:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
      break;
    case COARSE_INTERPOL:
34
35
      TEST_EXIT_DBG(feSpace)("Should not happen!\n");
      TEST_EXIT_DBG(feSpace->getBasisFcts())("Shoud not happen!\n");
36

37
38
39
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
      break;
    default:
40
41
      WARNING("Invalid coarsen operation \"%d\" in vector \"%s\"\n", 
	      coarsenOperation, this->name.c_str());
42
43
44
    }
  }

45

46
47
48
49
50
51
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList& list, int n)
  {
    (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->refineInter(this, &list, n);
  }

52

53
54
55
  template<>
  void DOFVector<WorldVector<double> >::refineInterpol(RCNeighbourList& list, int n)
  {
56
57
58
    if (n < 1) 
      return;

59
    Element *el = list.getElement(0);
60
    int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
61
62
63
    DegreeOfFreedom dof0 = el->getDof(0, n0);
    DegreeOfFreedom dof1 = el->getDof(1, n0);
    DegreeOfFreedom dof_new = el->getChild(0)->getDof(feSpace->getMesh()->getDim(), n0);
64
65
66
67
68
    vec[dof_new] = vec[dof0];
    vec[dof_new] += vec[dof1];
    vec[dof_new] *= 0.5;
  }

69

70
71
72
73
74
75
76
77
78
  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
79
    if (grad) {
80
81
      result = grad;
    } else {
82
      if(result && result->getFeSpace() != feSpace) {
83
84
	delete result;
	result = new DOFVector<WorldVector<double> >(feSpace, "gradient");
85
86
87
88
89
90
91
92
93
      }
    }

    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
94
95
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
96
    std::vector<DimVec<double>*> bary;
97

98
99
    int nNodes = 0;
    int nDofs = 0;
100

Thomas Witkowski's avatar
Thomas Witkowski committed
101
    for (int i = 0; i < dim + 1; i++) {
102
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
103
      int nPositions = mesh->getGeo(geoIndex);
104
      int numPreDofs = admin->getNumberOfPreDofs(i);
105
      for (int j = 0; j < nPositions; j++) {
106
	int dofs = basFcts->getNumberOfDofs(geoIndex);
107
	nNodeDOFs.push_back(dofs);
108
	nDofs += dofs;
109
	nNodePreDofs.push_back(numPreDofs);
110
      }
111
      nNodes += nPositions;
112
113
    }

114
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
115
      ("number of dofs != number of basis functions\n");
116
    
117
    for (int i = 0; i < nDofs; i++)
118
119
      bary.push_back(basFcts->getCoords(i));    

120
    ElementVector localUh(basFcts->getNumber());
121
122

    // traverse mesh
123
    std::vector<bool> visited(getUsedSize(), false);
124
    TraverseStack stack;
125
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
126
127
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

Thomas Witkowski's avatar
Thomas Witkowski committed
128
    while (elInfo) {
129
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
130
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
131
      getLocalVector(elInfo->getElement(), localUh);
132
133

      int localDOFNr = 0;
134
      for (int i = 0; i < nNodes; i++) { // for all nodes
135
136
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
Thomas Witkowski's avatar
Thomas Witkowski committed
137
138
139
	  if (!visited[dofIndex]) {	  
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, 
			       localUh, &((*result)[dofIndex]));
140
141
142
143
144
145
146
147
148
149
150
151
152

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }

153

154
155
156
157
158
159
  template<>
  DOFVector<WorldVector<double> >*
  DOFVector<double>::getRecoveryGradient(DOFVector<WorldVector<double> > *grad) const 
  {
    FUNCNAME("DOFVector<double>::getRecoveryGradient()");

160
161
162
163
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif

164
165
166
167
168
    // define result vector
    static DOFVector<WorldVector<double> > *vec = NULL;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
169
    if (!result) {
170
      if (vec && vec->getFeSpace() != feSpace) {
171
	delete vec;
172
173
	vec = NULL;
      }
174
      if (!vec)
175
	vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
176
      
177
178
179
180
181
182
183
184
185
      result = vec;
    }

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

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

    const BasisFunction *basFcts = feSpace->getBasisFcts();
186
    int nPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
187

188
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
Thomas Witkowski's avatar
Thomas Witkowski committed
189
    WorldVector<double> grd;
190
191

    // traverse mesh
Thomas Witkowski's avatar
Thomas Witkowski committed
192
    Mesh *mesh = feSpace->getMesh();
193
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
194
195
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET | 
196
					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
197

198
    ElementVector localUh(basFcts->getNumber());
199

Thomas Witkowski's avatar
Thomas Witkowski committed
200
    while (elInfo) {
201
      double det = elInfo->getDet();
202
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
203
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
204
      getLocalVector(elInfo->getElement(), localUh);
Thomas Witkowski's avatar
Thomas Witkowski committed
205
      basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
206

Thomas Witkowski's avatar
Thomas Witkowski committed
207
      for (int i = 0; i < dim + 1; i++) {
208
	DegreeOfFreedom dofIndex = dof[i][nPreDofs];
209
210
211
212
213
214
215
216
217
218
	(*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);

219
220
    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
      if (*volIt != 0.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
221
	*grdIt *= 1.0 / (*volIt);
222
223
224
225

    return result;
  }

226

227
  template<>
228
229
  const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo, 
								const Quadrature *quad,
230
								const FastQuadrature *quadFast,
231
								WorldVector<double> *grdAtQPs) const
232
233
234
  {
    FUNCNAME("DOFVector<double>::getGrdAtQPs()");
  
235
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
236

237
238
239
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
240
241
    }

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

    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
246
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
247
248
249
250
251
252
    static WorldVector<double> *grd = NULL;
    WorldVector<double> *result;

    if (grdAtQPs) {
      result = grdAtQPs;
    } else {
253
254
255
      if (grd) 
	delete [] grd;

Thomas Witkowski's avatar
Thomas Witkowski committed
256
      grd = new WorldVector<double>[nPoints];
257
      for (int i = 0; i < nPoints; i++)
258
	grd[i].set(0.0);
259
      
260
261
262
      result = grd;
    }
  
263
264
    ElementVector localVec(nBasFcts);
    getLocalVector(elInfo->getElement(), localVec);
265

266
    mtl::dense_vector<double> grd1(dim + 1);
267
268
269
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

270
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
271
      for (int i = 0; i < nPoints; i++) {
272
	grd1 = 0.0;
273

274
275
	for (int j = 0; j < nBasFcts; j++)
	  for (int k = 0; k < parts; k++)
276
	    grd1[k] += quadFast->getGradient(i, j, k) * localVec[j];
277
	  
278
279
	for (int l=0; l < dow; l++) {
	  result[i][l] = 0.0;
280
	  for (int k = 0; k < parts; k++)
281
282
283
	    result[i][l] += grdLambda[k][l] * grd1[k];
	}
      }
284
285

    } else {      
Thomas Witkowski's avatar
Thomas Witkowski committed
286
      const BasisFunction *basFcts = feSpace->getBasisFcts();
287
      mtl::dense_vector<double> grdPhi(dim + 1);
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
      for (int i = 0; i < nPoints; i++) {
290
	grd1 = 0.0;
291
292

	for (int j = 0; j < nBasFcts; j++) {
293
	  (*(basFcts->getGrdPhi(j)))(quad->getLambda(i), grdPhi);
294
	  for (int k = 0; k < parts; k++)
295
	    grd1[k] += grdPhi[k] * localVec[j];
296
297
	}

298
299
	for (int l = 0; l < dow; l++) {
	  result[i][l] = 0.0;
300
	  for (int k = 0; k < parts; k++)
301
302
303
304
	    result[i][l] += grdLambda[k][l] * grd1[k];
	}
      }
    }
305

306
307
308
    return const_cast<const WorldVector<double>*>(result);
  }

309

310
  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
311
312
313
314
315
316
317
318
319
320
  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");

321
    if (quad && quadFast) {
322
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
323
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339

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

    int dow = Global::getGeo(WORLD);
    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];
340
      for (int i = 0; i < nPoints; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
341
	grd[i].set(0.0);
342

Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
345
      result = grd;
    }
  
346
347
    ElementVector localVec(nBasFcts);
    getLocalVector(largeElInfo->getElement(), localVec);
348
349

    const BasisFunction *basFcts = feSpace->getBasisFcts();    
350
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353
354

    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();

355
356
357
    mtl::dense_vector<double> grd1(dim + 1);
    mtl::dense_vector<double> grdPhi(dim + 1);
    mtl::dense_vector<double> tmp(dim + 1);
Thomas Witkowski's avatar
Thomas Witkowski committed
358
359
    
    for (int i = 0; i < nPoints; i++) {
360
      grd1 = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
      
      for (int j = 0; j < nBasFcts; j++) {
363
	grdPhi = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
366

	for (int k = 0; k < nBasFcts; k++) {
	  (*(basFcts->getGrdPhi(k)))(quad->getLambda(i), tmp);
367
	  tmp *= m[j][k];
368
	  grdPhi += tmp;
Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
	}

371
	for (int k = 0; k < parts; k++)
372
	  grd1[k] += grdPhi[k] * localVec[j];
Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
375
376
      }
      
      for (int l = 0; l < dow; l++) {
	result[i][l] = 0.0;
377
	for (int k = 0; k < parts; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
378
379
380
381
382
383
384
	  result[i][l] += grdLambda[k][l] * grd1[k];
      }
    }    
  
    return const_cast<const WorldVector<double>*>(result);
  }

385

Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
  template<>
  const WorldMatrix<double> *DOFVectorBase<double>::getD2AtQPs(const ElInfo *elInfo, 
							       const Quadrature *quad,
389
390
391
392
393
							       const FastQuadrature *quadFast,
							       WorldMatrix<double>  *d2AtQPs) const
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
394
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
395

396
397
398
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
399
400
    }

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

404
405
    Element *el = elInfo->getElement(); 

406
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
407
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
408
409
410
411
412
413
414
415

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

    if (d2AtQPs) {
      result = d2AtQPs;
    } else {
      if(vec) delete [] vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
      vec = new WorldMatrix<double>[nPoints];
      for (int i = 0; i < nPoints; i++) {
418
419
420
421
422
	vec[i].set(0.0);
      }
      result = vec;
    }
  
423
424
    ElementVector localVec(nBasFcts);
    getLocalVector(el, localVec);
425
426
427
428
429

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

430
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
433
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
434
435
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438
439
	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);
440
441
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
442
443
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
444
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
445
446
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
447
448
449
450
		result[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
451
      const BasisFunction *basFcts = feSpace->getBasisFcts();
452
      DimMat<double> D2Phi(dim, NO_INIT);
453

Thomas Witkowski's avatar
Thomas Witkowski committed
454
455
456
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
457
458
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
465
	      D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
466
467
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
470
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
471
472
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
473
		result[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
474
475
476
	  }
      }
    }
477

478
479
480
    return const_cast<const WorldMatrix<double>*>(result);
  }

481

482
483
484
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
485
    FUNCNAME("DOFVector<double>::interpol()");
486

487
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
488
489
490
491
492
493
494
495
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

    int nBasisFcts = basisFcts->getNumber();
    int nSourceBasisFcts = sourceBasisFcts->getNumber();

    this->set(0.0);

496
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
497
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
498

499
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
500
501
502
503
504
505
      DimVec<double> *coords = NULL;
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
509
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
510
511
	source->getLocalVector(el, sourceLocalCoeffs);

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

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

564
565

  template<>
566
567
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
568
569
570
571
572
573
  {
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

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

576
    if (feSpace == vFeSpace)
577
578
579
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
580
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
581

582
    int nBasFcts = basFcts->getNumber();
583
584
    int vNumBasFcts = vBasFcts->getNumber();

585
586
    if (feSpace->getMesh() == vFeSpace->getMesh()) {      
      std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
587
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
588
589
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
590
591
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
592
593
594

      while (elInfo) {
	Element *el = elInfo->getElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
595
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
596
597
	v->getLocalVector(el, vLocalCoeffs);

598
	for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
599
	  if (vec[myLocalIndices[i]] == nul) {
600
	    coords = basFcts->getCoords(i);
601
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs, NULL) * factor;
602
603
604
605
606
607
608
609
610
611
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


612
613
614
615
616
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
    FUNCNAME("DOFVector<double>::getGradient()");

617
618
619
620
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif

621
622
623
624
625
626
627
628
629
630
631
    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;

632
    if (grad) {
633
634
      result = grad;
    } else {
635
      if (!result) {
636
	result = new WorldVector<DOFVector<double>*>;
637

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

    // count number of nodes and dofs per node
649
650
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
651
    std::vector<DimVec<double>*> bary;
652

653
654
    int nNodes = 0;
    int nDofs = 0;
655

656
    for (int i = 0; i < dim + 1; i++) {
657
658
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
659
      int numPreDofs = admin->getNumberOfPreDofs(i);
660
      for (int j = 0; j < numPositionNodes; j++) {
661
	int dofs = basFcts->getNumberOfDofs(geoIndex);
662
	nNodeDOFs.push_back(dofs);
663
	nDofs += dofs;
664
	nNodePreDofs.push_back(numPreDofs);
665
      }
666
      nNodes += numPositionNodes;
667
668
    }

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

    // traverse mesh
676
    std::vector<bool> visited(getUsedSize(), false);
677
    TraverseStack stack;
678
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
679
680
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
681
    ElementVector localUh(basFcts->getNumber());
682

683
    while (elInfo) {
684
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
685
      getLocalVector(elInfo->getElement(), localUh);
686
687
688
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
689
      for (int i = 0; i < nNodes; i++) { // for all nodes
690
691
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
692

693
	  if (!visited[dofIndex]) {
694
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
695

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


712
713
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
714
  {}
715

716

717
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
718
719
  {}

720
721
722
723

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

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

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

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

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

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

    return r;
  }

747

Thomas Witkowski's avatar
Thomas Witkowski committed
748
  template<>
749
750
751
752
753
  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
754
755
  {
    FUNCNAME("DOFVector<double>::getVecAtQPs()");
756
 
Thomas Witkowski's avatar
Thomas Witkowski committed
757
758
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");

759
    if (quad && quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
760
761
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
762
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
763
764
765
766

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

767
    if (smallElInfo->getMesh() == feSpace->getMesh())
768
      return getVecAtQPs(smallElInfo, quad, quadFast, vecAtQPs);    
769

Thomas Witkowski's avatar
Thomas Witkowski committed
770
    const BasisFunction *basFcts = feSpace->getBasisFcts();
771
772
    int nPoints = 
      quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
773

774
775
    vecAtQPs.change_dim(nPoints);
    
776
777
    ElementVector localVec(nBasFcts);
    getLocalVector(largeElInfo->getElement(), localVec);
778
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
779
780

    for (int i = 0; i < nPoints; i++) {
781
      vecAtQPs[i] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
782
783
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
784
	for (int k = 0; k < nBasFcts; k++)
785
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
786
787

	vecAtQPs[i] += localVec[j] * val;
Thomas Witkowski's avatar
Thomas Witkowski committed
788
789
790
791
      }
    }
  }

792

Thomas Witkowski's avatar
Thomas Witkowski committed
793
794
795
796
797
798
799
800
801
802
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    FUNCNAME("DOFVector::assemble()");

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

803
    set_to_zero(this->elementVector);
804
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
805
806
807

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
808
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
809
810
811
812
813
814
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

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

823
824
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
825
826
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
827

Thomas Witkowski's avatar
Thomas Witkowski committed
828
829
830
831
832
833
834
835
836
837
838
839
  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;

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


863
864
865
866
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct)
  {
867
    if (vec1.getFeSpace()->getMesh()