DOFVector.cc 30.9 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<>
Thomas Witkowski's avatar
Thomas Witkowski committed
71
72
  const double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
					       ElInfo *oldElInfo) const
73
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
74
    FUNCNAME("DOFVector<double>::evalAtPoint()");
75
76
77
78
79
80
81
  
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
82
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
83
84
85
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();
86
    static double value = 0.0;
87
88
89
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
90
91
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, 
				       oldElInfo->getMacroElement(), NULL, NULL);
92
93
94
95
96
97
98
99
      delete oldElInfo;
    } else
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Thomas Witkowski's avatar
Thomas Witkowski committed
100
101
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), 
			       localIndices);
102
103
104
105
106
107
108
109
110
111
112
      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."));

    if (oldElInfo == NULL)
      delete elInfo;

    return value;
Thomas Witkowski's avatar
Thomas Witkowski committed
113
  }
114
115
116


  template<>
Thomas Witkowski's avatar
Thomas Witkowski committed
117
118
119
  const WorldVector<double> 
  DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
					       ElInfo *oldElInfo) const
120
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
121
    FUNCNAME("DOFVector<double>::evalAtPoint()");
122
123
124
125
126
127
128
  
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
129
    vector<DegreeOfFreedom> localIndices(nBasFcts);
130
131
132
133
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();

Thomas Witkowski's avatar
Thomas Witkowski committed
134
    WorldVector<double> value(DEFAULT_VALUE, 0.0);
135
136
137
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
138
139
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, 
				       oldElInfo->getMacroElement(), NULL, NULL);
140
141
142
143
144
145
146
147
      delete oldElInfo;
    } else
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Thomas Witkowski's avatar
Thomas Witkowski committed
148
149
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), 
			       localIndices);
150
151
152
      mtl::dense_vector<WorldVector<double> > uh(lambda.size());
      for (int i = 0; i < lambda.size(); i++)
        uh[i] = operator[](localIndices[i]);
Thomas Witkowski's avatar
Thomas Witkowski committed
153
      value = basFcts->evalUh(lambda, uh);
154
155
156
157
158
159
    } else
      throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));

    if (oldElInfo == NULL)
      delete elInfo;

Thomas Witkowski's avatar
Thomas Witkowski committed
160
161
    return value;
  }
162
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
303
  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;
  }


304
  template<>
305
306
307
308
  void DOFVectorBase<double>::getD2AtQPs( const ElInfo *elInfo,
					  const Quadrature *quad,
					  const FastQuadrature *quadFast,
					  mtl::dense_vector<D2Type<double>::type> &d2AtQPs) const
309
310
311
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
312
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
313

314
315
316
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
317
318
    }

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

322
323
    Element *el = elInfo->getElement(); 

324
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
325
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
326
  
327
    mtl::dense_vector<double> localVec(nBasFcts);
328
    getLocalVector(el, localVec);
329
330
331
332
333

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

334
    d2AtQPs.change_dim(nPoints);
335
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
338
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
339
340
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
341
342
343
344
	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);
345
346
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
349
	    d2AtQPs[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
350
351
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
352
		d2AtQPs[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
353
354
355
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
356
      const BasisFunction *basFcts = feSpace->getBasisFcts();
357
      DimMat<double> D2Phi(dim, NO_INIT);
358

Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
361
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
362
363
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
368
369
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
370
	      D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
371
372
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
375
	    d2AtQPs[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
376
377
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
378
		d2AtQPs[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
379
380
381
382
383
	  }
      }
    }
  }

384

385
386
387
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {
388
    FUNCNAME("DOFVector<double>::interpol()");
389

390
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
391
392
393
394
395
396
397
398
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

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

    this->set(0.0);

399
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
400
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
401

402
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
403
404
405
406
407
408
      DimVec<double> *coords = NULL;
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

409
      while (elInfo) {
410
411
	Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
412
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
413
414
	source->getLocalVector(el, sourceLocalCoeffs);

415
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
416
	  if (vec[myLocalIndices[i]] == 0.0) {
417
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
418
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
	  }
	}
	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);
437
      while (nextTraverse) {     
438
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
439
				   myLocalIndices);
440
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
441

442
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
443
	  if (vec[myLocalIndices[i]] == 0.0) {
444
            elInfo1->coordToWorld(*(basisFcts->getCoords(i)), worldVec);
445
446
447
	    elInfo2->worldToCoord(worldVec, &coords2);
	  
	    bool isPositive = true;
448
	    for (int j = 0; j < coords2.size(); j++) {
449
	      if (coords2[j] < -0.00001) {
450
451
452
453
454
		isPositive = false;
		break;
	      }
	    }
	  
455
456
457
	    if (isPositive)
	      vec[myLocalIndices[i]] = 
		sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	   
458
459
460
	  }
	}
      
461
462
	nextTraverse = 
	  dualStack.traverseNext(&elInfo1, &elInfo2, &elInfoSmall, &elInfoLarge);
463
464
465
466
      }
    }
  }

467
468

  template<>
469
470
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
471
472
473
474
475
476
  {
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

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

479
    if (feSpace == vFeSpace)
480
481
482
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
483
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
484

485
    int nBasFcts = basFcts->getNumber();
486
487
    int vNumBasFcts = vBasFcts->getNumber();

488
489
    if (feSpace->getMesh() == vFeSpace->getMesh()) {      
      std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
490
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
491
492
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
493
494
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
495
496
497

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

501
	for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
502
	  if (vec[myLocalIndices[i]] == nul) {
503
	    coords = basFcts->getCoords(i);
504
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs) * factor;
505
506
507
508
509
510
511
512
513
514
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
      ERROR_EXIT("not yet for dual traverse\n");
    }
  }


515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
  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;

531
    if (grad) {
532
533
      result = grad;
    } else {
534
      if (!result) {
535
	result = new WorldVector<DOFVector<double>*>;
536

537
538
	result->set(NULL);
      }
539
      for (int i = 0; i < dow; i++) {
540
	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {
541
542
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
543
544
545
546
547
	}
      }
    }

    // count number of nodes and dofs per node
548
549
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
550
    std::vector<DimVec<double>*> bary;
551

552
553
    int nNodes = 0;
    int nDofs = 0;
554

555
    for (int i = 0; i < dim + 1; i++) {
556
557
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
558
      int numPreDofs = admin->getNumberOfPreDofs(i);
559
      for (int j = 0; j < numPositionNodes; j++) {
560
	int dofs = basFcts->getNumberOfDofs(geoIndex);
561
	nNodeDOFs.push_back(dofs);
562
	nDofs += dofs;
563
	nNodePreDofs.push_back(numPreDofs);
564
      }
565
      nNodes += numPositionNodes;
566
567
    }

568
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
569
      ("number of dofs != number of basis functions\n");
570
    
571
    for (int i = 0; i < nDofs; i++)
572
      bary.push_back(basFcts->getCoords(i));    
573
574

    // traverse mesh
575
    std::vector<bool> visited(getUsedSize(), false);
576
    TraverseStack stack;
577
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
578
579
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
580
    ElementVector localUh(basFcts->getNumber());
581

582
    while (elInfo) {
583
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
584
      getLocalVector(elInfo->getElement(), localUh);
585
586
587
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
588
      for (int i = 0; i < nNodes; i++) { // for all nodes
589
590
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
591

592
	  if (!visited[dofIndex]) {
593
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, grd);
594

595
596
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
597
598
599
600
601
602
603
604
605
606
607
608
609
610

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


611
612
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
613
  {}
614

615

616
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
617
618
  {}

619
620
621
622

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

625
    TEST_EXIT_DBG(vec)("no vector\n");
626

627
    int dow = Global::getGeo(WORLD);
628
629
    static WorldVector<DOFVector<double>*> *result = NULL;

630
    if (!res && !result) {
631
632
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
633
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
634
635
636
637
638
    }

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

    int vecSize = vec->getSize();
639
640
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
641
642
643
644
645
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

646
647


Thomas Witkowski's avatar
Thomas Witkowski committed
648
649
650
651
652
653
654
655
656
657
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    FUNCNAME("DOFVector::assemble()");

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

658
    set_to_zero(this->elementVector);
659
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
660
661
662

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
663
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
666
667
668
669
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	   it != this->operators.end(); 
670
	   ++it, ++factorIt)
671
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
672
	  (*it)->getElementVector(elInfo, this->elementVector, 
673
				  *factorIt ? **factorIt : 1.0);      
674
675
	  addVector = true;
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
    }

678
679
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
680
681
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
682

Thomas Witkowski's avatar
Thomas Witkowski committed
683
684
685
686
687
688
689
690
691
692
693
694
  template<>
  void DOFVectorBase<double>::assemble2(double factor, 
					ElInfo *mainElInfo, ElInfo *auxElInfo,
					ElInfo *smallElInfo, ElInfo *largeElInfo,
					const BoundaryType *bound, 
					Operator *op)
  {
    FUNCNAME("DOFVector::assemble2()");

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

695
    set_to_zero(this->elementVector);
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
    bool addVector = false;

    TEST_EXIT(!op)("Not yet implemented!\n");

    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator factorIt;
    for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();	
	 it != this->operators.end(); 
	 ++it, ++factorIt)
      if ((*it)->getNeedDualTraverse()) {
	(*it)->getElementVector(mainElInfo, auxElInfo,
				smallElInfo, largeElInfo,
				this->elementVector, 
				*factorIt ? **factorIt : 1.0);
	addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
711
712
      }
  
713
714
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
715
716
717
  }


718
719
720
721
  double integrate(const DOFVector<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct)
  {
722
    if (vec1.getFeSpace()->getMesh() == vec2.getFeSpace()->getMesh())
723
724
725
726
727
728
729
730
731
732
733
734
735
736
      return intSingleMesh(vec1, vec2, fct);
    else
      return intDualMesh(vec1, vec2, fct);
  }


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

737
738
    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
739
    Quadrature* quad = 
740
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
741
    FastQuadrature *fastQuad =
742
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
743

744
745
    mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
746
747
748
749

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
750
    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
751
    while (elInfo) {
752
753
754
755
      vec1.getVecAtQPs(elInfo, quad, fastQuad, qp1);
      vec2.getVecAtQPs(elInfo, quad, fastQuad, qp2);

      double tmp = 0.0;
756
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
757
758
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);     
      value += tmp * elInfo->getDet();
759
760
761
762
      
      elInfo = stack.traverseNext(elInfo);
    }

763
764
765
766
767
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
768
769
770
771
772
773
774
775
776
777
778
779
    return value;
  }


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

780
781
    int deg = std::max(fct->getDegree(), 
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
782
    Quadrature* quad = 
783
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
784
    FastQuadrature *fastQuad =
785
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
786

787
788
    mtl::dense_vector<double> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<double> qp2(fastQuad->getNumPoints());
789
790
791
792
793

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    DualTraverse dualTraverse;    
    DualElInfo dualElInfo;
794
795
    bool cont = dualTraverse.traverseFirst(vec1.getFeSpace()->getMesh(),
					   vec2.getFeSpace()->getMesh(),
796
797
798
					   -1, -1, traverseFlag, traverseFlag,
					   dualElInfo);
    while (cont) {
799
800
801
802
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);

      double tmp = 0.0;
803
      for (int iq = 0; iq < quad->getNumPoints(); iq++)
804
805
 	tmp += quad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);      
      value += tmp * dualElInfo.smallElInfo->getDet();
806
807
808
809
      
      cont = dualTraverse.traverseNext(dualElInfo);
    }

810
811
812
813
814
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif

815
816
817
    return value;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
818
819
820
821
822
823
824
825
  
  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate()");
    
    TEST_EXIT(fct)("No function defined!\n");

826
827
828
829
830
    const FiniteElemSpace *feSpace = vec.getFeSpace();
    Mesh *mesh = feSpace->getMesh();

    int deg = std::max(fct->getDegree(), 2 * feSpace->getBasisFcts()->getDegree());
    Quadrature* quad = Quadrature::provideQuadrature(mesh->getDim(), deg);
Thomas Witkowski's avatar
Thomas Witkowski committed
831
    FastQuadrature *fastQuad =
832
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);
Thomas Witkowski's avatar
Thomas Witkowski committed
833
834
835
836
837
838
839

    mtl::dense_vector<double> qp(fastQuad->getNumPoints());
    WorldVector<double> coords;

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
840
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
Thomas Witkowski's avatar
Thomas Witkowski committed
841
842
843
844
845
846
847
848
849
850
851
852
853
854
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);

      double tmp = 0.0;
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
 	tmp += fastQuad->getWeight(iq) * (qp[iq] * (*fct)(coords));
      }

      value += tmp * elInfo->getDet();
      
      elInfo = stack.traverseNext(elInfo);
    }

855
856
857
858
859
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
860
861
    return value;
  }
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894


  double integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<double, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate()");
    
    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(), feSpace->getBasisFcts()->getDegree());
    Quadrature* quad = 
      Quadrature::provideQuadrature(feSpace->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);

    WorldVector<double> coords;

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;      
    TraverseStack stack;    
    ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, traverseFlag);
    while (elInfo) {
      double tmp = 0.0;
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
 	tmp += fastQuad->getWeight(iq) * (*fct)(coords);
      }

      value += tmp * elInfo->getDet();
      
      elInfo = stack.traverseNext(elInfo);
    }

895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return value;
  }



  double integrateGeneral(const std::vector<DOFVector<double>*> &vecs,
		       AbstractFunction<double, std::vector<double> > *fct)
  {
    FUNCNAME("integrateGeneral()");

    TEST_EXIT(fct)("No function defined!\n");
    TEST_EXIT(vecs.size()>0)("No DOFVectors provided!\n");

    int deg = std::max(fct->getDegree(),
		       2 * vecs[0].getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad =
      Quadrature::provideQuadrature(vecs[0].getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vecs[0].getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

    std::vector<mtl::dense_vector<double> > qp(vecs.size());
    std::vector<double> qp_local(vecs.size());
    for (size_t i = 0; i < vecs.size(); i++)
      qp[i].change_dim(fastQuad->getNumPoints());

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
	vecs[i].getVecAtQPs(elInfo, quad, fastQuad, qp[i]);

      double tmp = 0.0;
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	for (size_t i = 0; i < vecs.size(); i++)
	  qp_local[i] = qp[i][iq];
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp_local);
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return value;
  }

  double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
				  const std::vector<DOFVector<double>*> &grds,
				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > *fct)
  {
    FUNCNAME("integrateGeneral()");

    TEST_EXIT(fct)("No function defined!\n");
    TEST_EXIT(vecs.size()>0)("No DOFVectors provided!\n");
    TEST_EXIT(grds.size()>0)("No DOFVectors for gradients provided!\n");

    int deg = std::max(fct->getDegree(),
		       2 * vecs[0].getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad =
      Quadrature::provideQuadrature(vecs[0].getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vecs[0].getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

    std::vector<mtl::dense_vector<double> > qp(vecs.size());
    std::vector<mtl::dense_vector<WorldVector<double> > > qpGrd(vecs.size());
    std::vector<double> qp_local(vecs.size());
    std::vector<WorldVector<double> > grd_local(grds.size());
    for (size_t i = 0; i < vecs.size(); i++)
      qp[i].change_dim(fastQuad->getNumPoints());
    for (size_t i = 0; i < grds.size(); i++)
      qpGrd[i].change_dim(fastQuad->getNumPoints());

    double value = 0.0;
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
	vecs[i].getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
      for (size_t i = 0; i < grds.size(); i++)
	grds[i].getGradientAtQPs(elInfo, quad, fastQuad, qpGrd[i]);

      double tmp = 0.0;
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	for (size_t i = 0; i < vecs.size(); i++)
	  qp_local[i] = qp[i][iq];
	for (size_t i = 0; i < grds.size(); i++)
	  grd_local[i] = qpGrd[i][iq];
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp_local, grd_local);
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

1001
1002
1003
1004
1005
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localValue = value;
    MPI::COMM_WORLD.Allreduce(&localValue, &value, 1, MPI_DOUBLE, MPI_SUM);
#endif

1006
1007
    return value;
  }
1008
1009
}