DOFVector.cc 27.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
#include <boost/numeric/mtl/mtl.hpp>
23
24
25
26
#include "DOFVector.h"
#include "Traverse.h"
#include "DualTraverse.h"
#include "FixVec.h"
27
#include "ElementFunction.h"
28
29
30
31
32
33

namespace AMDiS {

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

36
    switch (coarsenOperation) {
37
38
39
40
41
42
43
    case NO_OPERATION:
      return;
      break;
    case COARSE_RESTRICT:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
      break;
    case COARSE_INTERPOL:
44
45
      TEST_EXIT_DBG(feSpace)("Should not happen!\n");
      TEST_EXIT_DBG(feSpace->getBasisFcts())("Shoud not happen!\n");
46

47
48
49
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
      break;
    default:
50
51
      WARNING("Invalid coarsen operation \"%d\" in vector \"%s\"\n", 
	      coarsenOperation, this->name.c_str());
52
53
54
    }
  }

55

56
57
58
  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList& list, int n)
  {
59
60
61
62
63
64
65
66
67
    switch (refineOperation) {
    case NO_OPERATION:
      return;
      break;
    case REFINE_INTERPOL:
    default:
      (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->refineInter(this, &list, n);
      break;
    }
68
69
  }

70

71
72
73
  template<>
  void DOFVector<WorldVector<double> >::refineInterpol(RCNeighbourList& list, int n)
  {
74
75
76
    if (refineOperation == NO_OPERATION)
      return;
    
77
78
79
    if (n < 1) 
      return;

80
    Element *el = list.getElement(0);
81
    int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
82
83
84
    DegreeOfFreedom dof0 = el->getDof(0, n0);
    DegreeOfFreedom dof1 = el->getDof(1, n0);
    DegreeOfFreedom dof_new = el->getChild(0)->getDof(feSpace->getMesh()->getDim(), n0);
85
86
87
88
89
    vec[dof_new] = vec[dof0];
    vec[dof_new] += vec[dof1];
    vec[dof_new] *= 0.5;
  }

90

91
  template<>
92
93
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p,
					ElInfo *oldElInfo) const
94
  {  
95
96
97
98
99
100
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
101
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
102
103
104
    DimVec<double> lambda(dim, NO_INIT);
  
    ElInfo *elInfo = mesh->createNewElInfo();
105
    double value = 0.0;
106
107
108
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
109
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), nullptr, nullptr);
110
111
      delete oldElInfo;
    } else
112
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
113
114
115
116
117

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Praetorius, Simon's avatar
Praetorius, Simon committed
118
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
119
120
      ElementVector uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
121
122
        uh[i] = operator[](localIndices[i]);
      value = basFcts->evalUh(lambda, uh);
123
124
125
126
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      return 0.0;
#else
127
      ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
128
129
130
#endif
    }

131

132
    if (oldElInfo == nullptr)
133
134
135
      delete elInfo;

    return value;
136
  }
137
138
139


  template<>
140
141
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const
142
  {  
143
144
145
146
147
148
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
  
    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();
  
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
149
150
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
    DimVec<double> lambda(dim, NO_INIT);  
151
    ElInfo *elInfo = mesh->createNewElInfo();
152
    WorldVector<double> value(DEFAULT_VALUE, 0.0);
153
154
155
    bool inside = false;
  
    if (oldElInfo && oldElInfo->getMacroElement()) {
156
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), nullptr, nullptr);
157
158
      delete oldElInfo;
    } else
159
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
160
161
162
163
164

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
Praetorius, Simon's avatar
Praetorius, Simon committed
165
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
166
167
      mtl::dense_vector<WorldVector<double> > uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
168
        uh[i] = operator[](localIndices[i]);
169
      value = basFcts->evalUh(lambda, uh);
170
171
172
    } else {
      ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
    }
173

174
    if (oldElInfo == nullptr)
175
176
      delete elInfo;

177
178
    return value;
  }
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
  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;
230
          this->getVecAtQPs(elInfo, quadSurfaces[face], nullptr, uh_vec);
231
232
233
234
235
236
237
238
239
          for (int iq = 0; iq < nPoints; iq++)
            normT += quadSurfaces[face]->getWeight(iq) * (uh_vec[iq]);
          result += det * normT;
        }
      }
  
      elInfo = stack.traverseNext(elInfo);
    }
  
240
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
241
    Parallel::mpi::globalAdd(result);
242
#endif
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
    
    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);
300
          this->getVecAtQPs(elInfo, quadSurfaces[face], nullptr, uh_vec);
301
302
303
304
305
306
307
308
309
310
          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);
    }
  
311
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
312
    Parallel::mpi::globalAdd(result);
313
#endif
314
315
316
317
318
    
    return result;
  }


319
  template<>
320
321
322
323
  void DOFVectorBase<double>::getD2AtQPs( const ElInfo *elInfo,
					  const Quadrature *quad,
					  const FastQuadrature *quadFast,
					  mtl::dense_vector<D2Type<double>::type> &d2AtQPs) const
324
325
326
  {
    FUNCNAME("DOFVector<double>::getD2AtQPs()");
  
327
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
328

329
330
331
    if (quad && quadFast) {
      TEST_EXIT_DBG(quad == quadFast->getQuadrature())
      	("quad != quadFast->quadrature\n");
332
333
    }

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

337
338
    Element *el = elInfo->getElement(); 

339
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
340
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
341
  
342
    mtl::dense_vector<double> localVec(nBasFcts);
343
    getLocalVector(el, localVec);
344
345
346
347
348

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

349
    d2AtQPs.change_dim(nPoints);
350
    if (quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
354
355
	    D2Tmp[k][l] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
358
359
	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);
360
361
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
364
	    d2AtQPs[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
367
		d2AtQPs[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
368
369
370
	  }
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
371
      const BasisFunction *basFcts = feSpace->getBasisFcts();
372
      DimMat<double> D2Phi(dim, NO_INIT);
373

Thomas Witkowski's avatar
Thomas Witkowski committed
374
375
376
      for (int iq = 0; iq < nPoints; iq++) {
	for (int k = 0; k < parts; k++)
	  for (int l = 0; l < parts; l++)
377
378
	    D2Tmp[k][l] = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
383
384
	  for (int k = 0; k < parts; k++)
	    for (int l = 0; l < parts; l++)
385
	      D2Tmp[k][l] += localVec[i] * D2Phi[k][l];
386
387
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
388
389
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++) {
390
	    d2AtQPs[iq][i][j] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
	    for (int k = 0; k < parts; k++)
	      for (int l = 0; l < parts; l++)
393
		d2AtQPs[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
394
395
396
397
398
	  }
      }
    }
  }

399

400
401
402
403
  template<>
  void DOFVector<double>::interpol(DOFVector<double> *source, double factor) 
  {

404
    const FiniteElemSpace *sourceFeSpace = source->getFeSpace();
405
406
407
408
409
410
411
412
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    const BasisFunction *sourceBasisFcts = sourceFeSpace->getBasisFcts();

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

    this->set(0.0);

413
    std::vector<DegreeOfFreedom> myLocalIndices(nBasisFcts);
414
    ElementVector sourceLocalCoeffs(nSourceBasisFcts);
415

416
    if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
417
      DimVec<double> *coords = nullptr;
418
419
420
421
422
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS);

423
      if (basisFcts->isNodal()) {
424
      while (elInfo) {
425
426
	Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
427
	basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
428
429
	source->getLocalVector(el, sourceLocalCoeffs);

430
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
431
	  if (vec[myLocalIndices[i]] == 0.0) {
432
	    coords = basisFcts->getCoords(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
433
	    vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
434
435
436
437
438
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
	// nonnodal base:
	ElementFunctionWorld<double> F(source, factor);
	
	while (elInfo) {
	  F.setElInfo(elInfo);
	  
	  Element *el = elInfo->getElement();
	  
	  basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
	  
	  ElementVector rvec(nBasisFcts);	  
	  basisFcts->interpol(elInfo, nBasisFcts, NULL, &F, rvec);
	  
	  for (int i = 0; i < nBasisFcts; i++) {
	    if (vec[myLocalIndices[i]] == 0.0) {
	      vec[myLocalIndices[i]] = rvec[i];
	    }
	  }
	  elInfo = stack.traverseNext(elInfo);
	}	
      }
    } else {
      if (basisFcts->isNodal()) {
462
463
464
465
466
467
468
469
470
471
472
473
474
      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);
475
      while (nextTraverse) {     
476
	basisFcts->getLocalIndices(elInfo1->getElement(), feSpace->getAdmin(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
477
				   myLocalIndices);
478
	source->getLocalVector(elInfo2->getElement(), sourceLocalCoeffs);
479

480
	for (int i = 0; i < nBasisFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
481
	  if (vec[myLocalIndices[i]] == 0.0) {
482
            elInfo1->coordToWorld(*(basisFcts->getCoords(i)), worldVec);
483
484
485
	    elInfo2->worldToCoord(worldVec, &coords2);
	  
	    bool isPositive = true;
486
	    for (int j = 0; j < coords2.size(); j++) {
487
	      if (coords2[j] < -0.00001) {
488
489
490
491
492
		isPositive = false;
		break;
	      }
	    }
	  
493
494
495
	    if (isPositive)
	      vec[myLocalIndices[i]] = 
		sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);	   
496
497
498
	  }
	}
      
499
500
	nextTraverse = 
	  dualStack.traverseNext(&elInfo1, &elInfo2, &elInfoSmall, &elInfoLarge);
501
      }
502
503
504
505
      } else {
	// nonnodal base
	ERROR_EXIT("not yet implemented\n");
      }
506
507
508
    }
  }

509
510

  template<>
511
512
  void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, 
						 double factor) 
513
  {
514
    FUNCNAME("DOFVector<WorldVector<double> >::interpol()");
515
516
517
518
    WorldVector<double> nul(DEFAULT_VALUE,0.0);

    this->set(nul);

519
    DimVec<double> *coords = nullptr;
520
    const FiniteElemSpace *vFeSpace = v->getFeSpace();
521

522
    if (feSpace == vFeSpace)
523
524
525
      WARNING("same FE-spaces\n");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
526
    const BasisFunction *vBasFcts = vFeSpace->getBasisFcts();
527

528
    int nBasFcts = basFcts->getNumber();
529
530
    int vNumBasFcts = vBasFcts->getNumber();

531
532
    if (feSpace->getMesh() == vFeSpace->getMesh()) {      
      std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
533
      mtl::dense_vector<WorldVector<double> > vLocalCoeffs(vNumBasFcts);
534
535
      Mesh *mesh = feSpace->getMesh();
      TraverseStack stack;
536
537
      ElInfo *elInfo = 
	stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
538
      if (basFcts->isNodal()){
539
540
      while (elInfo) {
	Element *el = elInfo->getElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
541
	basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
542
543
	v->getLocalVector(el, vLocalCoeffs);

544
	for (int i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
545
	  if (vec[myLocalIndices[i]] == nul) {
546
	    coords = basFcts->getCoords(i);
547
	    vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs) * factor;
548
549
550
551
552
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
    } else {
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
	//nonnodal base:
	ElementFunctionWorld<WorldVector<double> > F(v, factor);
	
	while (elInfo) {
	  F.setElInfo(elInfo);
	  
	  Element *el = elInfo->getElement();
	  
	  basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
	  v->getLocalVector(el, vLocalCoeffs);
	  
	  mtl::dense_vector<WorldVector<double> > rvec(nBasFcts);
	  basFcts->interpol(elInfo, nBasFcts, NULL, &F, rvec);
	  
	  for (int i = 0; i < nBasFcts; i++) {
	    if (vec[myLocalIndices[i]] == nul) {
	      vec[myLocalIndices[i]] = rvec[i];
	    }
	  }
	  elInfo = stack.traverseNext(elInfo);
	}	
      }
    } else {
      ERROR_EXIT("not yet implemented for dual traverse\n");
577
578
579
580
    }
  }


581
582
583
  template<>
  WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
  {
584
    FUNCNAME_DBG("DOFVector<double>::getGradient()");
585
586
587
588
589
590
591
592
593
594

    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
595
    static WorldVector<DOFVector<double>*> *result = nullptr; // TODO: REMOVE STATIC
596

597
    if (grad) {
598
599
      result = grad;
    } else {
600
      if (!result) {
601
	result = new WorldVector<DOFVector<double>*>;
602

603
	result->set(nullptr);
604
      }
605
      for (int i = 0; i < dow; i++) {
606
	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {
607
608
	  delete (*result)[i];
	  (*result)[i] = new DOFVector<double>(feSpace, "gradient");
609
610
611
612
613
	}
      }
    }

    // count number of nodes and dofs per node
614
615
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
616
    std::vector<DimVec<double>*> bary;
617

618
619
    int nNodes = 0;
    int nDofs = 0;
620

621
    for (int i = 0; i < dim + 1; i++) {
622
623
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int numPositionNodes = mesh->getGeo(geoIndex);
624
      int numPreDofs = admin->getNumberOfPreDofs(i);
625
      for (int j = 0; j < numPositionNodes; j++) {
626
	int dofs = basFcts->getNumberOfDofs(geoIndex);
627
	nNodeDOFs.push_back(dofs);
628
	nDofs += dofs;
629
	nNodePreDofs.push_back(numPreDofs);
630
      }
631
      nNodes += numPositionNodes;
632
633
    }

634
    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
635
      ("number of dofs != number of basis functions\n");
636
    
637
    for (int i = 0; i < nDofs; i++)
638
      bary.push_back(basFcts->getCoords(i));    
639
640

    // traverse mesh
641
    std::vector<bool> visited(getUsedSize(), false);
642
    TraverseStack stack;
643
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
644
645
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
    WorldVector<double> grd;
646
    ElementVector localUh(basFcts->getNumber());
647

648
    while (elInfo) {
649
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
650
      getLocalVector(elInfo->getElement(), localUh);
651
652
653
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

      int localDOFNr = 0;
654
      for (int i = 0; i < nNodes; i++) { // for all nodes
655
656
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
657

658
	  if (!visited[dofIndex]) {
659
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, grd);
660

661
662
	    for (int k = 0; k < dow; k++)
	      (*(*result)[k])[dofIndex] = grd[k];	    
663
664
665
666
667
668
669
670
671
672
673
674
675
676

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

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


677
678
  DOFVectorDOF::DOFVectorDOF() : 
    DOFVector<DegreeOfFreedom>() 
Thomas Witkowski's avatar
Thomas Witkowski committed
679
  {}
680

681

682
  void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) 
683
684
  {}

685
686
687
688

  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *res)
  {
689
    FUNCNAME_DBG("DOFVector<double>::transform()");
690

691
    TEST_EXIT_DBG(vec)("no vector\n");
692

693
    int dow = Global::getGeo(WORLD);
694
    static WorldVector<DOFVector<double>*> *result = nullptr; // TODO: REMOVE STATIC
695

696
    if (!res && !result) {
697
698
      result = new WorldVector<DOFVector<double>*>;
      for (int i = 0; i < dow; i++)
699
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
700
701
702
703
704
    }

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

    int vecSize = vec->getSize();
705
706
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < dow; j++)
707
708
709
710
711
	(*((*r)[j]))[i] = (*vec)[i][j];

    return r;
  }

712
713


Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
716
717
718
719
720
721
  template<>
  void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
				       const BoundaryType *bound, 
				       Operator *op)
  {
    if (!(op || this->operators.size())) 
      return;

722
    set_to_zero(this->elementVector);
723
    bool addVector = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
724
725
726

    if (op) {
      op->getElementVector(elInfo, this->elementVector);
727
      addVector = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
728
729
730
731
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;

732
      for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
733
	   it != this->operators.end(); 
734
	   ++it, ++factorIt)
735
	if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
736
	  (*it)->getElementVector(elInfo, this->elementVector, 
737
				  *factorIt ? **factorIt : 1.0);      
738
	  addVector = true;	  
739
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
740
741
    }

742
743
    if (addVector)
      addElementVector(factor, this->elementVector, bound, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
746

Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
749
750
751
752
753
754
755
756
757
758
  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;

759
    set_to_zero(this->elementVector);
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
    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
775
776
      }
  
777
778
    if (addVector)
      addElementVector(factor, this->elementVector, bound, mainElInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
779
780
  }

781
782
783
784
785
  
  template<>
  void DOFVectorBase<double>::assembleOperator(Operator &op)
  {
    FUNCNAME("DOFVectorBase::assembleOperator()");
Thomas Witkowski's avatar
Thomas Witkowski committed
786

787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
    TEST_EXIT(op.getRowFeSpace() == feSpace)
      ("Row FE spaces do not fit together!\n");

    Mesh *mesh = feSpace->getMesh();
    mesh->dofCompress();
    const BasisFunction *basisFcts = feSpace->getBasisFcts();

    Flag assembleFlag = getAssembleFlag() | 
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_NEIGH |
      Mesh::FILL_BOUND;

    BoundaryType *bound = new BoundaryType[basisFcts->getNumber()];

    if (getBoundaryManager())
      getBoundaryManager()->initVector(this);


    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
    while (elInfo) {
      basisFcts->getBound(elInfo, bound);

      assemble(1.0, elInfo, bound, &op);

      if (getBoundaryManager())
	getBoundaryManager()->fillBoundaryConditions(elInfo, this);

      elInfo = stack.traverseNext(elInfo);
    }

    finishAssembling();
    getBoundaryManager()->exitVector(this);

    delete [] bound;
  }
  
  
828
829
830
831
832
833
834
835
836
  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(),
Praetorius, Simon's avatar
Praetorius, Simon committed
837
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
838
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
839
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
840
    FastQuadrature *fastQuad =
Praetorius, Simon's avatar
Praetorius, Simon committed
841
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad, INIT_PHI);
842
843
844
845
846
847
848
849
850

    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;
Praetorius, Simon's avatar
Praetorius, Simon committed
851
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
852
853
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
854
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
855
856
857
858
859
860
861
862
863
864
865
866
867

      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
868
    Parallel::mpi::globalAdd(value);
869
870
871
872
873
874
875
#endif

    return value;
  }

  double integrateGeneralGradient(const std::vector<DOFVector<double>*> &vecs,
				  const std::vector<DOFVector<double>*> &grds,
Praetorius, Simon's avatar
Praetorius, Simon committed
876
				  BinaryAbstractFunction<double, std::vector<double>, std::vector<WorldVector<double> > > *fct)
877
878
879
880
881
882
883
884
  {
    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(),
Praetorius, Simon's avatar
Praetorius, Simon committed
885
		       2 * vecs[0]->getFeSpace()->getBasisFcts()->getDegree());
886
    Quadrature* quad =
Praetorius, Simon's avatar
Praetorius, Simon committed
887
      Quadrature::provideQuadrature(vecs[0]->getFeSpace()->getMesh()->getDim(), deg);
888
    FastQuadrature *fastQuad =
889
      FastQuadrature::provideFastQuadrature(vecs[0]->getFeSpace()->getBasisFcts(), *quad,  INIT_PHI | INIT_GRD_PHI);
890
891
892
893
894
895
896
897
898
899
900

    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;
901
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;
902
    TraverseStack stack;
Praetorius, Simon's avatar
Praetorius, Simon committed
903
    ElInfo *elInfo = stack.traverseFirst(vecs[0]->getFeSpace()->getMesh(), -1, traverseFlag);
904
905
    while (elInfo) {
      for (size_t i = 0; i < vecs.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
906
	vecs[i]->getVecAtQPs(elInfo, quad, fastQuad, qp[i]);
907
      for (size_t i = 0; i < grds.size(); i++)
Praetorius, Simon's avatar
Praetorius, Simon committed
908
	grds[i]->getGrdAtQPs(elInfo, quad, fastQuad, qpGrd[i]);
909
910
911
912
913
914
915
916
917
918
919
920
921
922

      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);
    }

923
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
924
    Parallel::mpi::globalAdd(value);
925
926
#endif

927
928
    return value;
  }
929
930
}