Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

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 95
	int dofs = basFcts->getNumberOfDOFs(geoIndex);
	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
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
Thomas Witkowski's avatar
Thomas Witkowski committed
176
    WorldVector<double> grd;
177 178

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
194
      for (int i = 0; i < dim + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
195
	DegreeOfFreedom dofIndex = dof[i][nPreDOFs];
196 197 198 199 200 201 202 203 204 205
	(*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);

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

    return result;
  }

213

214
  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;
    }
  
251 252
    const ElementVector& localVec = localVectors[myRank];
    getLocalVector(elInfo->getElement(), const_cast<ElementVector&>(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
    return const_cast<const WorldVector<double>*>(result);
  }

299

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

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

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

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

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

    DimVec<double> &grd1 = *grdTmp[myRank];
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();

    DimVec<double>* grdPhi = grdPhis[myRank];
    DimVec<double> tmp(dim, DEFAULT_VALUE, 0.0);
    
    for (int i = 0; i < nPoints; i++) {
351
      for (int j = 0; j < dim + 1; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
352 353 354 355 356 357 358
	grd1[j] = 0.0;
      
      for (int j = 0; j < nBasFcts; j++) {

	grdPhi->set(0.0);
	for (int k = 0; k < nBasFcts; k++) {
	  (*(basFcts->getGrdPhi(k)))(quad->getLambda(i), tmp);
359
	  tmp *= m[j][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
360 361 362
	  *grdPhi += tmp;
	}

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

377

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

474

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

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

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

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

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

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

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

557 558

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

    this->set(nul);

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

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

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

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

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

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

	for (int i = 0; i < numBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
592
	  if (vec[myLocalIndices[i]] == nul) {
593
	    coords = basFcts->getCoords(i);
594
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs, NULL) * factor;
595 596 597 598 599 600 601 602 603 604
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


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

709

710
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
711 712
  {}

713 714 715 716

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

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

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

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

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

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

    return r;
  }

740

Thomas Witkowski's avatar
Thomas Witkowski committed
741
  template<>
742 743 744 745 746
  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
747 748
  {
    FUNCNAME("DOFVector<double>::getVecAtQPs()");
749
 
Thomas Witkowski's avatar
Thomas Witkowski committed
750 751
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");

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

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

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

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

767 768 769 770
    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
771

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

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

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

786

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

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

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

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

      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	   it != this->operators.end(); 
809
	   ++it, ++factorIt)
810
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
811
	  (*it)->getElementVector(elInfo, this->elementVector,