DOFVector.cc 37.6 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
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
  template<>
  double DOFVector<double>::IntOnBoundary(BoundaryType boundaryType, 
                                          Quadrature* q) const
  {
    FUNCNAME("DOFVector<T>::IntOnBoundary()");
  
    Mesh* mesh = this->feSpace->getMesh();
    int dim = mesh->getDim();
  
    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(dim - 1, deg);
    } else {
      TEST_EXIT(q->getDim() == dim-1)
        ("Provided quadrature has wrong dimension for surfaceQuadrature!\n");
    }
  
    // create barycentric coords for each vertex of each side
    const Element *refElement = Global::getReferenceElement(dim);
    VectorOfFixVecs<DimVec<double> >**coords;
    coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
  
    // for all element sides
    for (int i = 0; i < dim + 1; i++) {
      coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
        DimVec<double>(dim, DEFAULT_VALUE, 0.0));
      // for each vertex of the side
      for (int k = 0; k < dim; k++) {
        int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), 
                                                    i, k);
        (*coords[i])[k][index] = 1.0;
      }
    }
  
    std::vector<SurfaceQuadrature*> quadSurfaces(dim + 1);
    for (int i = 0; i < dim + 1; i++)
      quadSurfaces[i] = new SurfaceQuadrature(q, *coords[i]);
  
    double result = 0.0;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
      Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
    while (elInfo) {
      for (int face = 0; face < dim + 1; face++) {
        if (elInfo->getBoundary(face) == boundaryType) {
          int nPoints = quadSurfaces[face]->getNumPoints();
          mtl::dense_vector<double> uh_vec(nPoints);
          double det = elInfo->calcSurfaceDet(*coords[face]);
          double normT = 0.0;
          this->getVecAtQPs(elInfo, quadSurfaces[face], NULL, uh_vec);
          for (int iq = 0; iq < nPoints; iq++)
            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
          result += det * normT;
        }
      }
  
      elInfo = stack.traverseNext(elInfo);
    }
  
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
    #endif
    
    return result;
  }


  template<>
  double DOFVector<WorldVector<double> >::IntOnBoundaryNormal(
    BoundaryType boundaryType, Quadrature* q) const
  {
    FUNCNAME("DOFVector<T>::IntOnBoundary()");
 
    Mesh* mesh = this->feSpace->getMesh();
    int dim = mesh->getDim();
  
    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(dim - 1, deg);
    } else {
      TEST_EXIT(q->getDim() == dim-1)
        ("Provided quadrature has wrong dimension for surfaceQuadrature!\n");
    }
  
    // create barycentric coords for each vertex of each side
    const Element *refElement = Global::getReferenceElement(dim);
    VectorOfFixVecs<DimVec<double> >**coords;
    coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
  
    // for all element sides
    for (int i = 0; i < dim + 1; i++) {
      coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
        DimVec<double>(dim, DEFAULT_VALUE, 0.0));
      // for each vertex of the side
      for (int k = 0; k < dim; k++) {
        int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), 
                                                    i, k);
        (*coords[i])[k][index] = 1.0;
      }
    }
  
    std::vector<SurfaceQuadrature*> quadSurfaces(dim + 1);
    for (int i = 0; i < dim + 1; i++)
      quadSurfaces[i] = new SurfaceQuadrature(q, *coords[i]);
  
    double result = 0.0;
    WorldVector<double> normal;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
      Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
    while (elInfo) {
      for (int face = 0; face < dim + 1; face++) {
        if (elInfo->getBoundary(face) == boundaryType) {
          elInfo->getNormal(face, normal);
          double det = elInfo->calcSurfaceDet(*coords[face]);
  
          int nPoints = quadSurfaces[face]->getNumPoints();
          mtl::dense_vector<WorldVector<double> > uh_vec(nPoints);
          WorldVector<double> normT; normT.set(0.0);
          this->getVecAtQPs(elInfo, quadSurfaces[face], NULL, uh_vec);
          for (int iq = 0; iq < nPoints; iq++)
            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
          // scalar product between vector-valued solution and normal vector
          result += det * (normT * normal); 
        }
      }
  
      elInfo = stack.traverseNext(elInfo);
    }
  
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
    #endif
    
    return result;
  }


303
304
305
306
307
308
309
310
311
  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
312
    if (grad) {
313
314
      result = grad;
    } else {
315
      if(result && result->getFeSpace() != feSpace) {
316
317
	delete result;
	result = new DOFVector<WorldVector<double> >(feSpace, "gradient");
318
319
320
321
322
323
324
325
326
      }
    }

    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
327
328
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
329
    std::vector<DimVec<double>*> bary;
330

331
332
    int nNodes = 0;
    int nDofs = 0;
333

Thomas Witkowski's avatar
Thomas Witkowski committed
334
    for (int i = 0; i < dim + 1; i++) {
335
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
336
      int nPositions = mesh->getGeo(geoIndex);
337
      int numPreDofs = admin->getNumberOfPreDofs(i);
338
      for (int j = 0; j < nPositions; j++) {
339
	int dofs = basFcts->getNumberOfDofs(geoIndex);
340
	nNodeDOFs.push_back(dofs);
341
	nDofs += dofs;
342
	nNodePreDofs.push_back(numPreDofs);
343
      }
344
      nNodes += nPositions;
345
346
    }

347
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
348
      ("number of dofs != number of basis functions\n");
349
    
350
    for (int i = 0; i < nDofs; i++)
351
352
      bary.push_back(basFcts->getCoords(i));    

353
    ElementVector localUh(basFcts->getNumber());
354
355

    // traverse mesh
356
    std::vector<bool> visited(getUsedSize(), false);
357
    TraverseStack stack;
358
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
359
360
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    while (elInfo) {
362
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
363
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
364
      getLocalVector(elInfo->getElement(), localUh);
365
366

      int localDOFNr = 0;
367
      for (int i = 0; i < nNodes; i++) { // for all nodes
368
369
	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
370
371
372
	  if (!visited[dofIndex]) {	  
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, 
			       localUh, &((*result)[dofIndex]));
373
374
375
376
377
378
379
380
381
382
383
384
385

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }

386

387
388
389
390
391
392
393
394
395
396
397
  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
398
    if (!result) {
399
      if (vec && vec->getFeSpace() != feSpace) {
400
	delete vec;
401
402
	vec = NULL;
      }
403
      if (!vec)
404
	vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
405
      
406
407
408
409
410
411
412
413
414
      result = vec;
    }

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

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

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

417
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
Thomas Witkowski's avatar
Thomas Witkowski committed
418
    WorldVector<double> grd;
419
420

    // traverse mesh
Thomas Witkowski's avatar
Thomas Witkowski committed
421
    Mesh *mesh = feSpace->getMesh();
422
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
423
424
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET | 
425
					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
426

427
    ElementVector localUh(basFcts->getNumber());
428

Thomas Witkowski's avatar
Thomas Witkowski committed
429
    while (elInfo) {
430
      double det = elInfo->getDet();
431
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
432
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
433
      getLocalVector(elInfo->getElement(), localUh);
Thomas Witkowski's avatar
Thomas Witkowski committed
434
      basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
435

Thomas Witkowski's avatar
Thomas Witkowski committed
436
      for (int i = 0; i < dim + 1; i++) {
437
	DegreeOfFreedom dofIndex = dof[i][nPreDofs];
438
439
440
441
442
443
444
445
446
447
	(*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);

448
449
    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
      if (*volIt != 0.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
450
	*grdIt *= 1.0 / (*volIt);
451
452
453
454

    return result;
  }

455

456
  template<>
457
458
  const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo, 
								const Quadrature *quad,
459
								const FastQuadrature *quadFast,
460
								WorldVector<double> *grdAtQPs) const
461
462
463
  {
    FUNCNAME("DOFVector<double>::getGrdAtQPs()");
  
464
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
465

466
467
468
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
469
470
    }

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

    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
475
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
476
477
478
479
480
481
    static WorldVector<double> *grd = NULL;
    WorldVector<double> *result;

    if (grdAtQPs) {
      result = grdAtQPs;
    } else {
482
483
484
      if (grd) 
	delete [] grd;

Thomas Witkowski's avatar
Thomas Witkowski committed
485
      grd = new WorldVector<double>[nPoints];
486
      for (int i = 0; i < nPoints; i++)
487
	grd[i].set(0.0);
488
      
489
490
491
      result = grd;
    }
  
492
493
    ElementVector localVec(nBasFcts);
    getLocalVector(elInfo->getElement(), localVec);
494

495
    mtl::dense_vector<double> grd1(dim + 1);
496
497
498
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

499
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
500
      for (int i = 0; i < nPoints; i++) {
501
	grd1 = 0.0;
502

503
504
	for (int j = 0; j < nBasFcts; j++)
	  for (int k = 0; k < parts; k++)
505
	    grd1[k] += quadFast->getGradient(i, j, k) * localVec[j];
506
	  
507
508
	for (int l=0; l < dow; l++) {
	  result[i][l] = 0.0;
509
	  for (int k = 0; k < parts; k++)
510
511
512
	    result[i][l] += grdLambda[k][l] * grd1[k];
	}
      }
513
514

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

Thomas Witkowski's avatar
Thomas Witkowski committed
518
      for (int i = 0; i < nPoints; i++) {
519
	grd1 = 0.0;
520
521

	for (int j = 0; j < nBasFcts; j++) {
522
	  (*(basFcts->getGrdPhi(j)))(quad->getLambda(i), grdPhi);
523
	  for (int k = 0; k < parts; k++)
524
	    grd1[k] += grdPhi[k] * localVec[j];
525
526
	}

527
528
	for (int l = 0; l < dow; l++) {
	  result[i][l] = 0.0;
529
	  for (int k = 0; k < parts; k++)
530
531
532
533
	    result[i][l] += grdLambda[k][l] * grd1[k];
	}
      }
    }
534

535
536
537
    return const_cast<const WorldVector<double>*>(result);
  }

538

539
  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
540
541
542
543
544
545
546
547
548
549
  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");

550
    if (quad && quadFast) {
551
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
552
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568

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

Thomas Witkowski's avatar
Thomas Witkowski committed
572
573
574
      result = grd;
    }
  
575
576
    ElementVector localVec(nBasFcts);
    getLocalVector(largeElInfo->getElement(), localVec);
577
578

    const BasisFunction *basFcts = feSpace->getBasisFcts();    
579
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
580
581
582
583

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

584
585
586
    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
587
588
    
    for (int i = 0; i < nPoints; i++) {
589
      grd1 = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
590
591
      
      for (int j = 0; j < nBasFcts; j++) {
592
	grdPhi = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
593
594
595

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

600
	for (int k = 0; k < parts; k++)
601
	  grd1[k] += grdPhi[k] * localVec[j];
Thomas Witkowski's avatar
Thomas Witkowski committed
602
603
604
605
      }
      
      for (int l = 0; l < dow; l++) {
	result[i][l] = 0.0;
606
	for (int k = 0; k < parts; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
607
608
609
610
611
612
613
	  result[i][l] += grdLambda[k][l] * grd1[k];
      }
    }    
  
    return const_cast<const WorldVector<double>*>(result);
  }

614

Thomas Witkowski's avatar
Thomas Witkowski committed
615
616
617
  template<>
  const WorldMatrix<double> *DOFVectorBase<double>::getD2AtQPs(const ElInfo *elInfo, 
							       const Quadrature *quad,
618
619
620
621
622
							       const FastQuadrature *quadFast,
							       WorldMatrix<double>  *d2AtQPs) const
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
623
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
624

625
626
627
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
628
629
    }

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

633
634
    Element *el = elInfo->getElement(); 

635
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
636
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
637
638
639
640
641
642
643
644

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

    if (d2AtQPs) {
      result = d2AtQPs;
    } else {
      if(vec) delete [] vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
645
646
      vec = new WorldMatrix<double>[nPoints];
      for (int i = 0; i < nPoints; i++) {
647
648
649
650
651
	vec[i].set(0.0);
      }
      result = vec;
    }
  
652
653
    ElementVector localVec(nBasFcts);
    getLocalVector(el, localVec);
654
655
656
657
658

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

659
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
660
661
662
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
663
664
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
665
666
667
668
	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);
669
670
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
671
672
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
673
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
674
675
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
676
677
678
679
		result[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
680
      const BasisFunction *basFcts = feSpace->getBasisFcts();
681
      DimMat<double> D2Phi(dim, NO_INIT);
682

Thomas Witkowski's avatar
Thomas Witkowski committed
683
684
685
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
686
687
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
694
	      D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
695
696
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
697
698
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
699
	    result[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
700
701
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
702
		result[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
703
704
705
	  }
      }
    }
706

707
708
709
    return const_cast<const WorldMatrix<double>*>(result);
  }

710

711
712
713
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
714
    FUNCNAME("DOFVector<double>::interpol()");
715

716
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
717
718
719
720
721
722
723
724
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

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

    this->set(0.0);

725
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
726
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
727

728
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
729
730
731
732
733
734
      DimVec<double> *coords = NULL;
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

735
      while (elInfo) {
736
737
	Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
738
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
739
740
	source->getLocalVector(el, sourceLocalCoeffs);

741
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
742
	  if (vec[myLocalIndices[i]] == 0.0) {
743
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
744
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
	  }
	}
	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);
763
      while (nextTraverse) {     
764
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
765
				   myLocalIndices);
766
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
767

768
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
769
	  if (vec[myLocalIndices[i]] == 0.0) {
770
            elInfo1->coordToWorld(*(basisFcts->getCoords(i)), worldVec);
771
772
773
	    elInfo2->worldToCoord(worldVec, &coords2);
	  
	    bool isPositive = true;
774
	    for (int j = 0; j < coords2.size(); j++) {
775
	      if (coords2[j] < -0.00001) {
776
777
778
779
780
		isPositive = false;
		break;
	      }
	    }
	  
781
782
783
	    if (isPositive)
	      vec[myLocalIndices[i]] = 
		sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	   
784
785
786
	  }
	}
      
787
788
	nextTraverse = 
	  dualStack.traverseNext(&elInfo1, &elInfo2, &elInfoSmall, &elInfoLarge);
789
790
791
792
      }
    }
  }

793
794

  template<>
795
796
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
797
798
799
800
801
802
  {
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

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

805
    if (feSpace == vFeSpace)
806
807
808
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
809
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
810

811
    int nBasFcts = basFcts->getNumber();
812
813
    int vNumBasFcts = vBasFcts->getNumber();

814
815
    if (feSpace->getMesh() == vFeSpace->getMesh()) {      
      std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
816
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
817
818
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
819
820
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
821
822
823

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

827
	for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
828
	  if (vec[myLocalIndices[i]] == nul) {
829
	    coords = basFcts->getCoords(i);
830
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs, NULL) * factor;
831
832
833
834
835
836
837
838
839
840
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
  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;

857
    if (grad) {
858
859
      result = grad;
    } else {
860
      if (!result) {
861
	result = new WorldVector<DOFVector<double>*>;
862

863
864
	result->set(NULL);
      }
865
      for (int i = 0; i < dow; i++) {
866
	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {
867
868
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
869
870
871
872
873
	}
      }
    }

    // count number of nodes and dofs per node
874
875
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
876
    std::vector<DimVec<double>*> bary;
877

878
879
    int nNodes = 0;
    int nDofs = 0;
880

881
    for (int i = 0; i < dim + 1; i++) {
882
883
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
884
      int numPreDofs = admin->getNumberOfPreDofs(i);
885
      for (int j = 0; j < numPositionNodes; j++) {
886
	int dofs = basFcts->getNumberOfDofs(geoIndex);
887
	nNodeDOFs.push_back(dofs);
888
	nDofs += dofs;
889
	nNodePreDofs.push_back(numPreDofs);
890
      }
891
      nNodes += numPositionNodes;
892
893
    }

894
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
895
      ("number of dofs != number of basis functions\n");
896
    
897
    for (int i = 0; i < nDofs; i++)
898
      bary.push_back(basFcts->getCoords(i));    
899
900

    // traverse mesh
901
    std::vector<bool> visited(getUsedSize(), false);
902
    TraverseStack stack;
903
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
904
905
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
906
    ElementVector localUh(basFcts->getNumber());
907

908
    while (elInfo) {
909
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
910
      getLocalVector(elInfo->getElement(), localUh);
911
912
913
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
914
      for (int i = 0; i < nNodes; i++) { // for all nodes
915
916
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
917

918
	  if (!visited[dofIndex]) {
919
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
920

921
922
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
923
924
925
926
927
928
929
930
931
932
933
934
935
936

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


937
938
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
939
  {}
940

941

942
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
943
944
  {}

945
946
947
948

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

951
    TEST_EXIT_DBG(vec)("no vector\n");
952

953
    int dow = Global::getGeo(WORLD);
954
955
    static WorldVector<DOFVector<double>*> *result = NULL;

956
    if (!res && !result) {
957
958
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
959
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
960
961
962
963
964
    }

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

    int vecSize = vec->getSize();
965
966
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
967
968
969
970
971
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

972

Thomas Witkowski's avatar
Thomas Witkowski committed
973
  template<>
974
975
976
977
978
  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
979
980
  {
    FUNCNAME("DOFVector<double>::getVecAtQPs()");
981
 
Thomas Witkowski's avatar
Thomas Witkowski committed
982
983
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");

984
    if (quad && quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
985
986
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
987
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
988
989
990
991

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

992
    if (smallElInfo->getMesh() == feSpace->getMesh())
993
      return getVecAtQPs(smallElInfo, quad, quadFast, vecAtQPs);    
994

Thomas Witkowski's avatar
Thomas Witkowski committed
995
    const BasisFunction *basFcts = feSpace->getBasisFcts();
996
997
    int nPoints = 
      quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
998

999
1000
    vecAtQPs.change_dim(nPoints);
    
1001
1002
    ElementVector localVec(nBasFcts);
    getLocalVector(largeElInfo->getElement(), localVec);
1003
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
1004
1005

    for (int i = 0; i < nPoints; i++) {
1006
      vecAtQPs[i] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1007
1008
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
1009
	for (int k = 0; k < nBasFcts; k++)
1010
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
1011
1012

	vecAtQPs[i] += localVec[j] * val;
Thomas Witkowski's avatar
Thomas Witkowski committed
1013
1014
1015
1016
      }
    }
  }

1017

Thomas Witkowski's avatar
Thomas Witkowski committed
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    FUNCNAME("DOFVector::assemble()");

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

1028
    set_to_zero(this->elementVector);
1029
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
1030
1031
1032

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
1033
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
1034
1035
1036
1037
1038
1039
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	   it != this->operators.end(); 
1040
	   ++it, ++factorIt)
1041
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1042
	  (*it)->getElementVector(elInfo, this->elementVector, 
1043
				  *factorIt ? **factorIt : 1.0);      
1044
1045
	  addVector = true;
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
1046
1047
    }

1048
1049
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
1050
1051
  }