DOFVector.cc 32.7 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
  template<>
71
  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* values) const
72
73
74
75
76
77
78
79
80
81
82
83
84
  {
    FUNCNAME("DOFVector<double>::evalAtCoords()");
  
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();
85
    static double value = 0.0;
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
      delete oldElInfo;
    } else
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
      ElementVector uh(lambda.size());
      for (int i = 0; i < lambda.size(); i++)
        uh[i] = operator[](localIndices[i]);
      value = basFcts->evalUh(lambda, uh);
    } else
      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));

    delete [] localIndices;
    if (oldElInfo == NULL)
      delete elInfo;

    if (values != NULL)
      *values = value;
    return value;
  };


  template<>
117
  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
  {
    FUNCNAME("DOFVector<double>::evalAtCoords()");
  
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();

    static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
    WorldVector<double> *val = (NULL != values) ? values : &Values;

    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
      delete oldElInfo;
    } else
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
      mtl::dense_vector<WorldVector<double> > uh(lambda.size());
      for (int i = 0; i < lambda.size(); i++)
        uh[i] = operator[](localIndices[i]);
      basFcts->evalUh(lambda, uh, val);
    } else
      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));

    delete [] localIndices;
    if (oldElInfo == NULL)
      delete elInfo;

    return ((*val));
  };


163
164
165
166
167
168
169
170
171
  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
172
    if (grad) {
173
174
      result = grad;
    } else {
175
      if(result && result->getFeSpace() != feSpace) {
176
177
	delete result;
	result = new DOFVector<WorldVector<double> >(feSpace, "gradient");
178
179
180
181
182
183
184
185
186
      }
    }

    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
187
188
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
189
    std::vector<DimVec<double>*> bary;
190

191
192
    int nNodes = 0;
    int nDofs = 0;
193

Thomas Witkowski's avatar
Thomas Witkowski committed
194
    for (int i = 0; i < dim + 1; i++) {
195
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
196
      int nPositions = mesh->getGeo(geoIndex);
197
      int numPreDofs = admin->getNumberOfPreDofs(i);
198
      for (int j = 0; j < nPositions; j++) {
199
	int dofs = basFcts->getNumberOfDofs(geoIndex);
200
	nNodeDOFs.push_back(dofs);
201
	nDofs += dofs;
202
	nNodePreDofs.push_back(numPreDofs);
203
      }
204
      nNodes += nPositions;
205
206
    }

207
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
208
      ("number of dofs != number of basis functions\n");
209
    
210
    for (int i = 0; i < nDofs; i++)
211
212
      bary.push_back(basFcts->getCoords(i));    

213
    ElementVector localUh(basFcts->getNumber());
214
215

    // traverse mesh
216
    std::vector<bool> visited(getUsedSize(), false);
217
    TraverseStack stack;
218
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
219
220
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

Thomas Witkowski's avatar
Thomas Witkowski committed
221
    while (elInfo) {
222
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
223
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
224
      getLocalVector(elInfo->getElement(), localUh);
225
226

      int localDOFNr = 0;
227
      for (int i = 0; i < nNodes; i++) { // for all nodes
228
229
	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
230
231
232
	  if (!visited[dofIndex]) {	  
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, 
			       localUh, &((*result)[dofIndex]));
233
234
235
236
237
238
239
240
241
242
243
244
245

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }

246

247
248
249
250
251
252
253
254
255
256
257
  template<>
  DOFVector<WorldVector<double> >*
  DOFVector<double>::getRecoveryGradient(DOFVector<WorldVector<double> > *grad) const 
  {
    FUNCNAME("DOFVector<double>::getRecoveryGradient()");

    // define result vector
    static DOFVector<WorldVector<double> > *vec = NULL;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
258
    if (!result) {
259
      if (vec && vec->getFeSpace() != feSpace) {
260
	delete vec;
261
262
	vec = NULL;
      }
263
      if (!vec)
264
	vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
265
      
266
267
268
269
270
271
272
273
274
      result = vec;
    }

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

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

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

277
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    WorldVector<double> grd;
279
280

    // traverse mesh
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    Mesh *mesh = feSpace->getMesh();
282
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
283
284
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET | 
285
					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
286

287
    ElementVector localUh(basFcts->getNumber());
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
    while (elInfo) {
290
      double det = elInfo->getDet();
291
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
292
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
293
      getLocalVector(elInfo->getElement(), localUh);
Thomas Witkowski's avatar
Thomas Witkowski committed
294
      basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
295

Thomas Witkowski's avatar
Thomas Witkowski committed
296
      for (int i = 0; i < dim + 1; i++) {
297
	DegreeOfFreedom dofIndex = dof[i][nPreDofs];
298
299
300
301
302
303
304
305
306
307
	(*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);

308
309
    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
      if (*volIt != 0.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
310
	*grdIt *= 1.0 / (*volIt);
311
312
313
314

    return result;
  }

315

316
  template<>
317
318
  const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo, 
								const Quadrature *quad,
319
								const FastQuadrature *quadFast,
320
								WorldVector<double> *grdAtQPs) const
321
322
323
  {
    FUNCNAME("DOFVector<double>::getGrdAtQPs()");
  
324
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
325

326
327
328
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
329
330
    }

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

    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
335
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
336
337
338
339
340
341
    static WorldVector<double> *grd = NULL;
    WorldVector<double> *result;

    if (grdAtQPs) {
      result = grdAtQPs;
    } else {
342
343
344
      if (grd) 
	delete [] grd;

Thomas Witkowski's avatar
Thomas Witkowski committed
345
      grd = new WorldVector<double>[nPoints];
346
      for (int i = 0; i < nPoints; i++)
347
	grd[i].set(0.0);
348
      
349
350
351
      result = grd;
    }
  
352
353
    ElementVector localVec(nBasFcts);
    getLocalVector(elInfo->getElement(), localVec);
354

355
    mtl::dense_vector<double> grd1(dim + 1);
356
357
358
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

359
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
360
      for (int i = 0; i < nPoints; i++) {
361
	grd1 = 0.0;
362

363
364
	for (int j = 0; j < nBasFcts; j++)
	  for (int k = 0; k < parts; k++)
365
	    grd1[k] += quadFast->getGradient(i, j, k) * localVec[j];
366
	  
367
368
	for (int l=0; l < dow; l++) {
	  result[i][l] = 0.0;
369
	  for (int k = 0; k < parts; k++)
370
371
372
	    result[i][l] += grdLambda[k][l] * grd1[k];
	}
      }
373
374

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

Thomas Witkowski's avatar
Thomas Witkowski committed
378
      for (int i = 0; i < nPoints; i++) {
379
	grd1 = 0.0;
380
381

	for (int j = 0; j < nBasFcts; j++) {
382
	  (*(basFcts->getGrdPhi(j)))(quad->getLambda(i), grdPhi);
383
	  for (int k = 0; k < parts; k++)
384
	    grd1[k] += grdPhi[k] * localVec[j];
385
386
	}

387
388
	for (int l = 0; l < dow; l++) {
	  result[i][l] = 0.0;
389
	  for (int k = 0; k < parts; k++)
390
391
392
393
	    result[i][l] += grdLambda[k][l] * grd1[k];
	}
      }
    }
394

395
396
397
    return const_cast<const WorldVector<double>*>(result);
  }

398

399
  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
402
403
404
405
406
407
408
409
  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");

410
    if (quad && quadFast) {
411
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
412
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428

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

Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
434
      result = grd;
    }
  
435
436
    ElementVector localVec(nBasFcts);
    getLocalVector(largeElInfo->getElement(), localVec);
437
438

    const BasisFunction *basFcts = feSpace->getBasisFcts();    
439
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
442
443

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

444
445
446
    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
447
448
    
    for (int i = 0; i < nPoints; i++) {
449
      grd1 = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
450
451
      
      for (int j = 0; j < nBasFcts; j++) {
452
	grdPhi = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
453
454
455

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

460
	for (int k = 0; k < parts; k++)
461
	  grd1[k] += grdPhi[k] * localVec[j];
Thomas Witkowski's avatar
Thomas Witkowski committed
462
463
464
465
      }
      
      for (int l = 0; l < dow; l++) {
	result[i][l] = 0.0;
466
	for (int k = 0; k < parts; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
467
468
469
470
471
472
473
	  result[i][l] += grdLambda[k][l] * grd1[k];
      }
    }    
  
    return const_cast<const WorldVector<double>*>(result);
  }

474

Thomas Witkowski's avatar
Thomas Witkowski committed
475
476
477
  template<>
  const WorldMatrix<double> *DOFVectorBase<double>::getD2AtQPs(const ElInfo *elInfo, 
							       const Quadrature *quad,
478
479
480
481
482
							       const FastQuadrature *quadFast,
							       WorldMatrix<double>  *d2AtQPs) const
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
483
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
484

485
486
487
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
488
489
    }

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

493
494
    Element *el = elInfo->getElement(); 

495
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
496
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
497
498
499
500
501
502
503
504

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

    if (d2AtQPs) {
      result = d2AtQPs;
    } else {
      if(vec) delete [] vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
505
506
      vec = new WorldMatrix<double>[nPoints];
      for (int i = 0; i < nPoints; i++) {
507
508
509
510
511
	vec[i].set(0.0);
      }
      result = vec;
    }
  
512
513
    ElementVector localVec(nBasFcts);
    getLocalVector(el, localVec);
514
515
516
517
518

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

519
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
520
521
522
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
523
524
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
525
526
527
528
	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);
529
530
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
531
532
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
533
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
534
535
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
536
537
538
539
		result[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
540
      const BasisFunction *basFcts = feSpace->getBasisFcts();
541
      DimMat<double> D2Phi(dim, NO_INIT);
542

Thomas Witkowski's avatar
Thomas Witkowski committed
543
544
545
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
546
547
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
552
553
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
554
	      D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
555
556
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
557
558
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
559
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
560
561
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
562
		result[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
563
564
565
	  }
      }
    }
566

567
568
569
    return const_cast<const WorldMatrix<double>*>(result);
  }

570

571
572
573
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
574
    FUNCNAME("DOFVector<double>::interpol()");
575

576
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
577
578
579
580
581
582
583
584
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

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

    this->set(0.0);

585
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
586
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
587

588
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
589
590
591
592
593
594
      DimVec<double> *coords = NULL;
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

595
      while (elInfo) {
596
597
	Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
598
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
599
600
	source->getLocalVector(el, sourceLocalCoeffs);

601
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
602
	  if (vec[myLocalIndices[i]] == 0.0) {
603
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
604
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
	  }
	}
	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);
623
      while (nextTraverse) {     
624
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
625
				   myLocalIndices);
626
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
627

628
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
629
	  if (vec[myLocalIndices[i]] == 0.0) {
630
            elInfo1->coordToWorld(*(basisFcts->getCoords(i)), worldVec);
631
632
633
	    elInfo2->worldToCoord(worldVec, &coords2);
	  
	    bool isPositive = true;
634
	    for (int j = 0; j < coords2.size(); j++) {
635
	      if (coords2[j] < -0.00001) {
636
637
638
639
640
		isPositive = false;
		break;
	      }
	    }
	  
641
642
643
	    if (isPositive)
	      vec[myLocalIndices[i]] = 
		sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	   
644
645
646
	  }
	}
      
647
648
	nextTraverse = 
	  dualStack.traverseNext(&elInfo1, &elInfo2, &elInfoSmall, &elInfoLarge);
649
650
651
652
      }
    }
  }

653
654

  template<>
655
656
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
657
658
659
660
661
662
  {
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

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

665
    if (feSpace == vFeSpace)
666
667
668
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
669
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
670

671
    int nBasFcts = basFcts->getNumber();
672
673
    int vNumBasFcts = vBasFcts->getNumber();

674
675
    if (feSpace->getMesh() == vFeSpace->getMesh()) {      
      std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
676
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
677
678
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
679
680
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
681
682
683

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

687
	for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
688
	  if (vec[myLocalIndices[i]] == nul) {
689
	    coords = basFcts->getCoords(i);
690
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs, NULL) * factor;
691
692
693
694
695
696
697
698
699
700
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
    FUNCNAME("DOFVector<double>::getGradient()");

    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;

717
    if (grad) {
718
719
      result = grad;
    } else {
720
      if (!result) {
721
	result = new WorldVector<DOFVector<double>*>;
722

723
724
	result->set(NULL);
      }
725
      for (int i = 0; i < dow; i++) {
726
	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {
727
728
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
729
730
731
732
733
	}
      }
    }

    // count number of nodes and dofs per node
734
735
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
736
    std::vector<DimVec<double>*> bary;
737

738
739
    int nNodes = 0;
    int nDofs = 0;
740

741
    for (int i = 0; i < dim + 1; i++) {
742
743
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
744
      int numPreDofs = admin->getNumberOfPreDofs(i);
745
      for (int j = 0; j < numPositionNodes; j++) {
746
	int dofs = basFcts->getNumberOfDofs(geoIndex);
747
	nNodeDOFs.push_back(dofs);
748
	nDofs += dofs;
749
	nNodePreDofs.push_back(numPreDofs);
750
      }
751
      nNodes += numPositionNodes;
752
753
    }

754
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
755
      ("number of dofs != number of basis functions\n");
756
    
757
    for (int i = 0; i < nDofs; i++)
758
      bary.push_back(basFcts->getCoords(i));    
759
760

    // traverse mesh
761
    std::vector<bool> visited(getUsedSize(), false);
762
    TraverseStack stack;
763
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
764
765
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
766
    ElementVector localUh(basFcts->getNumber());
767

768
    while (elInfo) {
769
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
770
      getLocalVector(elInfo->getElement(), localUh);
771
772
773
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
774
      for (int i = 0; i < nNodes; i++) { // for all nodes
775
776
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
777

778
	  if (!visited[dofIndex]) {
779
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
780

781
782
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
783
784
785
786
787
788
789
790
791
792
793
794
795
796

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


797
798
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
799
  {}
800

801

802
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
803
804
  {}

805
806
807
808

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

811
    TEST_EXIT_DBG(vec)("no vector\n");
812

813
    int dow = Global::getGeo(WORLD);
814
815
    static WorldVector<DOFVector<double>*> *result = NULL;

816
    if (!res && !result) {
817
818
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
819
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
820
821
822
823
824
    }

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

    int vecSize = vec->getSize();
825
826
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
827
828
829
830
831
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

832

Thomas Witkowski's avatar
Thomas Witkowski committed
833
  template<>
834
835
836
837
838
  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
839
840
  {
    FUNCNAME("DOFVector<double>::getVecAtQPs()");
841
 
Thomas Witkowski's avatar
Thomas Witkowski committed
842
843
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");

844
    if (quad && quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
845
846
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
847
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
848
849
850
851

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

852
    if (smallElInfo->getMesh() == feSpace->getMesh())
853
      return getVecAtQPs(smallElInfo, quad, quadFast, vecAtQPs);    
854

Thomas Witkowski's avatar
Thomas Witkowski committed
855
    const BasisFunction *basFcts = feSpace->getBasisFcts();
856
857
    int nPoints = 
      quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
858

859
860
    vecAtQPs.change_dim(nPoints);
    
861
862
    ElementVector localVec(nBasFcts);
    getLocalVector(largeElInfo->getElement(), localVec);
863
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
864
865

    for (int i = 0; i < nPoints; i++) {
866
      vecAtQPs[i] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
867
868
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
869
	for (int k = 0; k < nBasFcts; k++)
870
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
871
872

	vecAtQPs[i] += localVec[j] * val;
Thomas Witkowski's avatar
Thomas Witkowski committed
873
874
875
876
      }
    }
  }

877

Thomas Witkowski's avatar
Thomas Witkowski committed
878
879
880
881
882
883
884
885
886
887
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    FUNCNAME("DOFVector::assemble()");

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

888
    set_to_zero(this->elementVector);
889
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
890
891
892

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
893
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
894
895
896
897
898
899
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	   it != this->operators.end(); 
900
	   ++it, ++factorIt)
901
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
902
	  (*it)->getElementVector(elInfo, this->elementVector, 
903
				  *factorIt ? **factorIt : 1.0);