Recovery.cc 29.8 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 "Recovery.h"
23
24
25
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/MeshDistributor.h"
#endif
26

27
28
RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) 
{ 
29
  if (rhs.coords) {
30
31
    if (!coords) 
      coords = new WorldVector<double>;    
Thomas Witkowski's avatar
Thomas Witkowski committed
32
33
    *coords = *rhs.coords;
  } else {
34
    if (coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
35
      delete coords;
36
      coords = nullptr;
37
38
39
40
41
    }
  }

  if (rhs.A) {
    if (!A) 
Thomas Witkowski's avatar
Thomas Witkowski committed
42
      A = new Matrix<double>(rhs.A->getNumRows(), rhs.A->getNumCols());
Thomas Witkowski's avatar
Thomas Witkowski committed
43
44
    *A = *rhs.A;
  } else {
45
    if (A) {
Thomas Witkowski's avatar
Thomas Witkowski committed
46
      delete A;
47
      A = nullptr;
48
49
50
51
52
    }
  }

  if (rhs.rec_uh) {
    if (!rec_uh)
Thomas Witkowski's avatar
Thomas Witkowski committed
53
      rec_uh = new Vector<double>(rhs.rec_uh->getSize());
Thomas Witkowski's avatar
Thomas Witkowski committed
54
55
    *rec_uh = *rhs.rec_uh;
  } else {
56
    if (rec_uh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
57
      delete rec_uh;
58
      rec_uh = nullptr;
59
60
61
62
63
    }
  }

  if (rhs.rec_grdUh) {
    if (!rec_grdUh)
Thomas Witkowski's avatar
Thomas Witkowski committed
64
      rec_grdUh = new Vector<WorldVector<double> >(rhs.rec_grdUh->getSize());
Thomas Witkowski's avatar
Thomas Witkowski committed
65
66
    *rec_grdUh = *rhs.rec_grdUh;
  } else {
67
    if (rec_grdUh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
68
      delete rec_grdUh;
69
      rec_grdUh = nullptr;
70
71
72
73
74
    }
  }

  if (rhs.neighbors) {
    if (!neighbors)
Thomas Witkowski's avatar
Thomas Witkowski committed
75
      neighbors = new std::set<DegreeOfFreedom>;
76
    *neighbors = *rhs.neighbors;
Thomas Witkowski's avatar
Thomas Witkowski committed
77
  } else {
78
    if (neighbors) {
Thomas Witkowski's avatar
Thomas Witkowski committed
79
      delete neighbors;
80
      neighbors = nullptr;
81
82
83
84
    }
  }

  return *this;
85
}
86
87
88
89


void RecoveryStructure::print()
{
90
  FUNCNAME("RecoveryStructure::print()");
91

92
  std::cout << std::endl;
93
94

  MSG("Coordinates of the node: ");
95
  std::cout << *coords << std::endl;
96

97
98
99
100
101
102
103
104
105
106
107
108
109
110
  if (A) {
    MSG("Interior vertex: printing system information.\n\n");
    
    int n = A->getNumRows();
    
    MSG("System matrix:\n");
    for (int i = 0; i < n; i++) {      
      MSG("( ");
      for (int j = 0; j < i; j++)
	std::cout << "* ";
      for (int j = i; j < n; j++)
	std::cout << (*A)[i][j] << " ";
      std::cout << ")" << std::endl;
    }
111

112
113
114
115
116
117
118
119
    MSG("Right hand side:\n");
    for (int i = 0; i < n; i++) {      
      MSG("( ");
      if (rec_grdUh)
	std::cout << (*rec_grdUh)[i];
      else
	std::cout << (*rec_uh)[i];
      std::cout << " )" << std::endl;
120
    }
121
122
123

    if (neighbors) {
      MSG("Printing neighbors vertices\n\n");
124
125

      MSG("Number of neighbors: ");
126
      std::cout << neighbors->size() << std::endl << std::endl;
127
      
128
      MSG("List of neighbors: ");
129
      std::set<DegreeOfFreedom>::const_iterator setIterator;
130
131
      
      for (setIterator = neighbors->begin(); setIterator != neighbors->end();
132
	   ++setIterator)
133
	std::cout << " " << *setIterator;
134
      
135
      std::cout << std::endl << std::endl;
136
    }
137
138
  } else {
    MSG("Boundary vertex or not a vertex node: printing interior neighbors\n\n");
139

140
141
142
143
144
145
146
147
148
    MSG("Number of neighbors: ");
    std::cout << neighbors->size() << std::endl << std::endl;
    
    MSG("List of neighbors: ");
    std::set<DegreeOfFreedom>::const_iterator setIterator;
    
    for (setIterator  = neighbors->begin(); setIterator != neighbors->end();
	 ++setIterator)
      std::cout << " " << *setIterator;
149

150
151
152
153
    std::cout << std::endl << std::endl;
  }
  
  WAIT;
154
155
156
157
158
}


void Recovery::set_feSpace(const FiniteElemSpace *fe_space)
{
159
160
161
  if (!feSpace || feSpace != fe_space) {
    if (struct_vec) {
      delete struct_vec;
162
      struct_vec = nullptr;
163
164
    }

165
166
167
168
169
    feSpace = fe_space;
    
    // create new structure vector
    struct_vec = new DOFVector<RecoveryStructure>(feSpace, "struct vec");
  }
170
171
}

172

173
174
int Recovery::set_exponents(int degree)
{
175
  FUNCNAME("Recovery::set_exponents()");
176
177
178
179
180

  int dow = Global::getGeo(WORLD);
  int number_monomials = degree + 1;

  // Computing number of monomials.
181
182
183
184
  if (dow > 1) {
    number_monomials *= (degree + 2);
    number_monomials /= 2;
  }
185

186
187
188
189
  if (dow == 3) {
    number_monomials *= (degree + 3);
    number_monomials /= 3;
  }
190
191

  // Allocating memory.
192
193
194
195
196
197
198
199
200
  if (n_monomials != number_monomials) {
    n_monomials = number_monomials;
    exponents.resize(n_monomials);
    
    if (matrix_fcts)
      matrix_fcts->resize(n_monomials, n_monomials);
    else
      matrix_fcts = new Matrix<Monomial*>(n_monomials, n_monomials);
  }
201
202

  // Setting vector of exponents.
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
  int count = 0;

  switch (dow) {    
  case 1:     // 1D monomials.
    for (int i = 0; i <= degree; i++)
      exponents[count++][0] = i;
    break;

  case 2:     // 2D monomials.
    for (int i = 0; i <= degree; i++)
      for (int j = 0; j <= i; j++) {
	exponents[count][0] = i - j;
	exponents[count++][1] = j;
      }
    break;
    
  case 3:     // 3D monomials.
    for (int i = 0; i <= degree; i++)
      for (int j = 0; j <= i; j++)
	for (int k = 0; k <= j; k++) {
	  exponents[count][0] = i - j;
	  exponents[count][1] = j - k;
	  exponents[count++][2] = k;
	}
    break;
228

229
230
231
  default:
    ERROR_EXIT("Which dimension have your world???\n");
  }
232

233
234
  TEST_EXIT_DBG(count == n_monomials)("There must be an error!\n");
  
235
236
237
  // Setting matrix of monomials.
  WorldVector<int> sum;

238
239
240
241
242
243
  for (int i = 0; i < n_monomials; i++)
    for (int j = i; j < n_monomials; j++) {
      // Computing exponent vector of monomial.
      sum = exponents[i] + exponents[j];
      (*matrix_fcts)[i][j] = new Monomial(sum);
    }
244
245
246
247

  return n_monomials;
}

248

249
250
251
252
253
254
void Recovery::compute_integrals(DOFVector<double> *uh, ElInfo *elInfo,
				 RecoveryStructure *rec_struct,
				 AbstractFunction<double, WorldVector<double> > *f_vec,
				 AbstractFunction<double, double> *f_scal,
				 DOFVector<double> *aux_vec)
{
255
  FUNCNAME_DBG("Recovery::compute_integrals()");
256

257
  TEST_EXIT_DBG(!(f_vec && f_scal))("Only one diffusion function, please!\n");
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272

  WorldVector<double> vec_sum;
  WorldVector<double> quad_pts;  // For world coordinates of quadrature points.

  int deg_f = 0;
  if (f_vec)
    deg_f = f_vec->getDegree();
  if (f_scal)
    deg_f = f_scal->getDegree();

  if (gradient)
    deg_f += feSpace->getBasisFcts()->getDegree() - 1;
  else
    deg_f += feSpace->getBasisFcts()->getDegree();

273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
  for (int i = 0; i < n_monomials; i++) {
    // Computing contributions to system matrix.
    for (int j = i; j < n_monomials; j++) {
      double sum  = 0.0;
      Quadrature *quad = Quadrature::provideQuadrature(Global::getGeo(WORLD),
						       (*matrix_fcts)[i][j]->getDegree());
      int n_points = quad->getNumPoints();

      for (int k = 0; k < n_points; k++) {
	elInfo->coordToWorld(quad->getLambda(k), quad_pts);
	sum += quad->getWeight(k) *
	  (*(*matrix_fcts)[i][j])(quad_pts, *rec_struct->coords);
      }
      (*(rec_struct->A))[i][j] += sum * elInfo->getDet();
    }
288

289
290
    Quadrature *quad = Quadrature::
      provideQuadrature(Global::getGeo(WORLD),
291
			  (*matrix_fcts)[0][i]->getDegree() + deg_f);
292
    int n_points = quad->getNumPoints();
293
    mtl::dense_vector<double> uhAtQP(n_points);
294
295
296
297
298
299
    
    // Computing contributions to right hand side.
    if (gradient) {    // For gradient recovery.
      double fAtQP = 1.0;
      if (f_scal) {
	if (aux_vec)
300
	  aux_vec->getVecAtQPs(elInfo, quad, nullptr, uhAtQP);
301
	else
302
	  uh->getVecAtQPs(elInfo, quad, nullptr, uhAtQP);
303
304
305
      }

      // Get gradient at quadrature points
306
      mtl::dense_vector<WorldVector<double> > grdAtQP(n_points);
307
      uh->getGrdAtQPs(elInfo, quad, nullptr, grdAtQP);
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
      vec_sum = 0.0;
      for (int k = 0; k < n_points; k++) {
	elInfo->coordToWorld(quad->getLambda(k), quad_pts);
	if (f_vec)
	  fAtQP = (*f_vec)(quad_pts);
	if (f_scal)
	  fAtQP = (*f_scal)(uhAtQP[k]);
	
	vec_sum = vec_sum + grdAtQP[k] * fAtQP * quad->getWeight(k)
	  * (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
      }
      (*rec_struct->rec_grdUh)[i] = (*rec_struct->rec_grdUh)[i]
	+ vec_sum * elInfo->getDet();
    } else {         // For recovery of DOFVector.
      // Get uh at quadrature points
323
      uh->getVecAtQPs(elInfo, quad, nullptr, uhAtQP);
324
325
326
327
328
329
330
      double sum = 0.0;
      for (int k = 0; k < n_points; k++) {
	elInfo->coordToWorld(quad->getLambda(k), quad_pts);
	sum += uhAtQP[k] * quad->getWeight(k)
	  * (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
      }
      (*rec_struct->rec_uh)[i] += sum * elInfo->getDet();
331
    }
332
  }
333
334
}

335

336
337
338
339
340
341
void Recovery::compute_interior_sums(DOFVector<double> *uh, ElInfo *elInfo,
				     RecoveryStructure *rec_struct, Quadrature *quad,
				     AbstractFunction<double, WorldVector<double> > *f_vec,
				     AbstractFunction<double, double> *f_scal,
				     DOFVector<double> *aux_vec)
{
342
  FUNCNAME_DBG("Recovery::compute_sums()");
343

344
345
  TEST_EXIT_DBG(gradient)("SPR of solution need computing node sums.\n");
  TEST_EXIT_DBG(!(f_vec && f_scal))("Only one diffusion function, please!\n");
346
347
348
349

  WorldVector<double> vec_sum;
  int n_points = quad->getNumPoints();
  WorldVector<double> quad_pts;  // For world coordinates of quadrature points.
350
  mtl::dense_vector<double> uhAtQP(n_points);
351
  mtl::dense_vector<WorldVector<double> > grdAtQP(n_points);
352

353
354
355
356
357
358
359
360
361
362
  for (int i = 0; i < n_monomials; i++) {
    // Computing contributions to system matrix.
    for (int j = i; j < n_monomials; j++) {
      double sum = 0.0;
      for (int k = 0; k < n_points; k++) {
	elInfo->coordToWorld(quad->getLambda(k), quad_pts);
	sum += (*(*matrix_fcts)[i][j])(quad_pts, *rec_struct->coords);
      }
      (*(rec_struct->A))[i][j] += sum;
    }
363

364
365
366
367
    // Computing contributions to right hand side.
    double fAtQP = 1.0;
    if (f_scal) {
      if (aux_vec)
368
	aux_vec->getVecAtQPs(elInfo, quad, nullptr, uhAtQP);
369
      else
370
	uh->getVecAtQPs(elInfo, quad, nullptr, uhAtQP);
371
372
373
    }

    // Get gradient at quadrature points
374
    uh->getGrdAtQPs(elInfo, quad, nullptr, grdAtQP);
375
376
377
378
379
    vec_sum = 0.0;
    for (int k = 0; k < n_points; k++) {
      elInfo->coordToWorld(quad->getLambda(k), quad_pts);
      if (f_vec)
	fAtQP = (*f_vec)(quad_pts);
380
      if (f_scal)
381
	fAtQP = (*f_scal)(uhAtQP[k]);
382

383
384
      vec_sum = vec_sum + grdAtQP[k] * fAtQP
	* (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
385
    }
386
387
    (*rec_struct->rec_grdUh)[i] = (*rec_struct->rec_grdUh)[i] + vec_sum;
  }
388
389
}

390

391
void Recovery::compute_node_sums(DOFVector<double> *uh, ElInfo *elInfo,
392
				 RecoveryStructure *rec_struct, DimVec<int> preDofs,
393
394
				 int n_vertices, int n_edges, int n_faces)
{
395
  FUNCNAME_DBG("Recovery::compute_sums()");
396

397
  TEST_EXIT_DBG(!gradient)
398
    ("SPR of flux or gradient need computing interior sums\n");
399
  TEST_EXIT_DBG(feSpace->getMesh()->getDim() == 1)
400
401
402
    ("At the moment only for linear finite elements.\n");

  WorldVector<double> node;  // For world coordinates at nodes.
403
  const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
404
405
  ElementVector uh_loc(n_vertices);
  uh->getLocalVector(elInfo->getElement(), uh_loc);
406

407
408
  for (int l = 0; l < n_vertices; l++) {
    // Computing contributions of vertex nodes
409
    if (rec_struct->neighbors->insert(dof[l][preDofs[VERTEX]]).second) {
410
411
412
413
414
415
416
417
418
419
420
      node = elInfo->getCoord(l);
      for (int i = 0; i < n_monomials; i++) {
	// Computing contributions to system matrix.
	for (int j = i; j < n_monomials; j++)
	  (*(rec_struct->A))[i][j] += (*(*matrix_fcts)[i][j])(node,
							      *rec_struct->coords);
	
	// Computing contributions to right hand side.
	(*rec_struct->rec_uh)[i] += uh_loc[l]
	  * (*(*matrix_fcts)[0][i])(node, *rec_struct->coords);
      }
421
    }
422
  }
423
424
}

425

426
427
void Recovery::compute_sums_linear(DOFVector<double> *uh, ElInfo *elInfo,
				   RecoveryStructure *rec_struct,
428
				   int vertex, DimVec<int> preDofs,
429
430
				   int n_vertices)
{
431
  FUNCNAME_DBG("Recovery::compute_sums_linear()");
432

433
  TEST_EXIT_DBG(!gradient)
434
435
436
    ("SPR of flux or gradient need computing interior sums\n");

  WorldVector<double> node;     // For world coordinates at nodes.
437
  const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
438
439
  ElementVector uh_loc(n_vertices);
  uh->getLocalVector(elInfo->getElement(), uh_loc);
440

441
442
  for (int l = 0;  l < n_vertices; l++) {
    // Computing contributions of vertex nodes
443
    DegreeOfFreedom k = dof[l][preDofs[VERTEX]];
444
445
446
447
448
449
450
451
452
453
454
455
    if (rec_struct->neighbors->insert(k).second) {
      node = elInfo->getCoord(l);
      for (int i = 0; i < n_monomials; i++) {
	// Computing contributions to system matrix.
	for (int j = i; j < n_monomials; j++)
	  (*(rec_struct->A))[i][j] += (*(*matrix_fcts)[i][j])(node,
							      *rec_struct->coords);
	
	// Computing contributions to right hand side.
	(*rec_struct->rec_uh)[i] += uh_loc[l]
	  * (*(*matrix_fcts)[0][i])(node, *rec_struct->coords);
      }
456
    }
457
  }
458

459
460
  if (vertex > 1 && elInfo->getNeighbour(vertex)) {
    int oppVertex = elInfo->getOppVertex(vertex);
461
    DegreeOfFreedom k = elInfo->getNeighbour(vertex)->getDof(oppVertex)[preDofs[VERTEX]];
462
463
464
465
466
467
468
469
470
471
    
    if (rec_struct->neighbors->insert(k).second) {
      node = elInfo->getOppCoord(vertex);
      for (int i = 0; i < n_monomials; i++) {
	// Computing contributions to system matrix.
	for (int j = i; j < n_monomials; j++)
	  (*(rec_struct->A))[i][j] += (*(*matrix_fcts)[i][j])(node, *rec_struct->coords);
	
	// Computing contributions to right hand side.
	(*rec_struct->rec_uh)[i] += (*uh)[k] * (*(*matrix_fcts)[0][i])(node, *rec_struct->coords);
472
      }
473
474
    }
  }  
475
476
}

477

478
479
480
481
482
void Recovery::fill_struct_vec(DOFVector<double> *uh,
			       AbstractFunction<double, WorldVector<double> > *f_vec,
			       AbstractFunction<double, double> *f_scal,
			       DOFVector<double> *aux_vec)
{
483
  FUNCNAME_DBG("Recovery::fill_struct_vec()");
484
485
486

  // Information on the mesh.
  Mesh *mesh = feSpace->getMesh();
487
  int dim = mesh->getDim();
488
489
490

  // Geometric information.
  int n_vertices = Global::getGeo(VERTEX, dim);
491
492
  int n_edges = Global::getGeo(EDGE, dim);
  int n_faces = Global::getGeo(FACE, dim);
493
494
495

  // Information concerning the finite element space.
  const BasisFunction *basis_fcts = feSpace->getBasisFcts();
496
  DimVec<int> *nDOFs = basis_fcts->getNumberOfDofs();
497
498
499

  // Information from DOFAdmin.
  const DOFAdmin *admin = feSpace->getAdmin();
500
  DimVec<int> preDofs = admin->getNumberOfPreDofs();
501
502
503
504
505
506
507

  // Variables for storing temporary information.
  DimVec<DegreeOfFreedom> interior_vertices(dim);
  WorldVector<double> coordinates;

  // Variables for passing information to integration routines.
  int degree = basis_fcts->getDegree();
508
  Quadrature *quad = nullptr;
509
510
511
512
513
  if (gradient && !method)
    quad = Quadrature::provideQuadrature(Global::getGeo(WORLD), degree);

  DimVec<int> pre_dofs(dim, NO_INIT);
  if (!gradient)
514
    pre_dofs = uh->getFeSpace()->getAdmin()->getNumberOfPreDofs();
515
  
516
  // Variables for traversing the mesh.
517
  Flag fill_flag = 
518
519
520
521
    Mesh::CALL_LEAF_EL |
    Mesh::FILL_COORDS |
    Mesh::FILL_BOUND |
    Mesh::FILL_DET |
522
    Mesh::FILL_GRD_LAMBDA;
523

524
  if (degree == 2 && dim > 1)
525
526
    fill_flag |= Mesh::FILL_NEIGH | Mesh::FILL_OPP_COORDS;

527
528
529
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  DofContainer boundaryDofs;
  DofContainerSet boundaryDofSet;
530
  Parallel::MeshDistributor::globalMeshDistributor->getAllBoundaryDofs(feSpace, 0, boundaryDofs);
531
532
533
  boundaryDofSet.insert(boundaryDofs.begin(), boundaryDofs.end());
#endif

534
  TraverseStack stack;
535
  ElInfo *el_info = stack.traverseFirst(mesh, -1, fill_flag);
536

537
  while (el_info) {    // traversing the mesh.    
538
    const DegreeOfFreedom **dof = el_info->getElement()->getDof();
539
540
541
    
    int n_neighbors = 0;     // counting interior vertices of element
    for (int i = 0; i < n_vertices; i++) {
542
      DegreeOfFreedom k = dof[i][preDofs[VERTEX]];
543
544
545
546
547
548
549
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      bool isInterior = (el_info->getBoundary(VERTEX, i) == INTERIOR) && (boundaryDofSet.count(dof[i]) == 0);
#else
      bool isInterior = el_info->getBoundary(VERTEX, i) == INTERIOR;
#endif

      if (isInterior)
550
551
552
553
554
555
556
	interior_vertices[n_neighbors++] = k;
    }

    TEST_EXIT_DBG(n_neighbors)
      ("Each element should have a least one interior vertex!\n");

    for (int i = 0; i < n_vertices; i++) {     // Handling nodes on vertices.
557
558
559
560
561
562
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      bool isInterior = (el_info->getBoundary(VERTEX, i) == INTERIOR) && (boundaryDofSet.count(dof[i]) == 0);
#else
      bool isInterior = el_info->getBoundary(VERTEX, i) == INTERIOR;
#endif

563
      DegreeOfFreedom k = dof[i][preDofs[VERTEX]];
564
565
566
567
      
      // Setting world coordinates of node.
      if (!(*struct_vec)[k].coords) {	
	(*struct_vec)[k].coords  = new WorldVector<double>;
568
	*(*struct_vec)[k].coords = el_info->getCoord(i);
569
      }
570

571
      if (isInterior) {	
572
573
574
575
576
577
578
579
580
581
582
583
584
	// Allocating memory for matrix and right hand side.
	if (!(*struct_vec)[k].A) {
	  (*struct_vec)[k].A = new Matrix<double>(n_monomials, n_monomials);
	  *(*struct_vec)[k].A = 0.0;
	  
	  if (gradient) {	    
	    (*struct_vec)[k].rec_grdUh = new Vector<WorldVector<double> >(n_monomials);
	    for (int j = 0; j < n_monomials; j++)
	      (*(*struct_vec)[k].rec_grdUh)[j] = 0.0;
	  } else {
	    (*struct_vec)[k].rec_uh = new Vector<double>(n_monomials);
	    *(*struct_vec)[k].rec_uh = 0.0;
	  }
585
586
	}

587
588
	// Computing the integrals.
	if (method)
589
	  compute_integrals(uh, el_info, &(*struct_vec)[k],
590
591
			    f_vec, f_scal, aux_vec);
	else if (gradient)
592
	  compute_interior_sums(uh, el_info, &(*struct_vec)[k], quad,
593
594
595
596
597
598
				f_vec, f_scal, aux_vec);
	else {
	  if (!(*struct_vec)[k].neighbors)
	    (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;

	  if (degree == 2 && dim > 1)
599
	    compute_sums_linear(uh, el_info, &(*struct_vec)[k],
600
601
				i, pre_dofs, n_vertices);
	  else
602
	    compute_node_sums(uh, el_info, &(*struct_vec)[k],
603
			      pre_dofs, n_vertices, n_edges, n_faces);
604
	}
605
606
607
608
609
610
611
612
      } else {     // Setting list of adjacent interior vertices.
	if (!(*struct_vec)[k].neighbors)
	  (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
	
	for (int j = 0; j < n_neighbors; j++)
	  (*struct_vec)[k].neighbors->insert(interior_vertices[j]);
      }
    }
613

614
    int n_count = n_vertices;
615

616
617
618
619
    if (dim > 1) {     // Handling nodes on edges.
      for (int i = 0; i < n_edges; i++) {
	bool isEdgeInterior = el_info->getBoundary(EDGE, i) == INTERIOR;

620
	for (int j = 0; j < (*nDOFs)[EDGE]; j++) {
621
	  DegreeOfFreedom k = dof[n_vertices + i][preDofs[EDGE] + j];
622
623

	  if (!(*struct_vec)[k].coords) {
624
	    // Setting world coordinates of node.
625
	    el_info->coordToWorld(*basis_fcts->getCoords(n_count),
626
627
				  coordinates);
	    (*struct_vec)[k].coords = new WorldVector<double>;
628
	    *(*struct_vec)[k].coords = coordinates;
629
	    
630
	    // Setting list of adjacent interior vertices.
Thomas Witkowski's avatar
Thomas Witkowski committed
631
	    (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
632
	    
633
	    if (isEdgeInterior) {
634
635
	      for (int m = 0; m < 2; m++) {
		int l = Global::getReferenceElement(dim)->getVertexOfEdge(i, m);
636
		if (el_info->getBoundary(VERTEX, l) == INTERIOR)
637
		  (*struct_vec)[k].neighbors->insert(dof[l][preDofs[VERTEX]]);
638
639
	      } 
	    } else {
640
641
	      for (int m = 0; m < n_neighbors; m++)
		(*struct_vec)[k].neighbors->insert(interior_vertices[m]);
642
	    }
643
	  }
644
645
646
	  
	  n_count++;
	}
647
648
      }
    }
649

650
651
652
    if (dim == 3)     // Handling nodes on faces.
      for (int i = 0; i < n_faces; i++)
	for (int j = 0; j < (*nDOFs)[FACE]; j++) {
653
	  DegreeOfFreedom k = dof[n_vertices+n_edges + i][preDofs[FACE] + j];
654
655
656
	  
	  if (!(*struct_vec)[k].coords) {
	    // Setting world coordinates of node.
657
	    el_info->coordToWorld(*basis_fcts->getCoords(n_count),
658
659
660
661
662
663
664
				  coordinates);
	    (*struct_vec)[k].coords  = new WorldVector<double>;
	    *(*struct_vec)[k].coords = coordinates;
	    
	    // Setting list of adjacent interior vertices.
	    (*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
	    
665
	    if (el_info->getBoundary(FACE, i) == INTERIOR)
666
667
	      for (int m = 0; m < 3; m++) {
		int l = Global::getReferenceElement(dim)->getVertexOfPosition(FACE, i, m);
668
		if (el_info->getBoundary(VERTEX, l) == INTERIOR)
669
		  (*struct_vec)[k].neighbors->insert(dof[l][preDofs[VERTEX]]);
670
671
672
673
674
675
676
677
678
679
680
	      }
	    else
	      for (int m = 0; m < n_neighbors; m++)
		(*struct_vec)[k].neighbors->insert(interior_vertices[m]);
	  }
	  
	  n_count++;
	}
    
    if ((*nDOFs)[CENTER])    // Handling nodes on center of element.
      for (int j = 0; j < (*nDOFs)[CENTER]; j++) {
681
	DegreeOfFreedom k = dof[n_vertices+n_edges+n_faces][preDofs[CENTER] + j];
682
683

	// Setting world coordinates of node.
684
	el_info->coordToWorld(*basis_fcts->getCoords(n_count), coordinates);
685
686
687
688
689
690
691
692
693
694
695
696
	(*struct_vec)[k].coords = new WorldVector<double>;
	*(*struct_vec)[k].coords = coordinates;

	// Setting list of adjacent interior vertices.
	(*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
	
	for (int m = 0; m < n_neighbors; m++)
	  (*struct_vec)[k].neighbors->insert(interior_vertices[m]);
	
	n_count++;
      }
    
697
    el_info = stack.traverseNext(el_info);
698
  }
699
700
}

701

702
703
void Recovery::recoveryUh(DOFVector<double> *uh, DOFVector<double> &rec_vec)
{
704
  FUNCNAME("Recovery::recoveryUh()");
705
706
707
708
709

  clear();

  gradient = false;

710
  const FiniteElemSpace *fe_space = rec_vec.getFeSpace();
711
712
713
714
715
716
717
  set_feSpace(fe_space);                                  // Setting feSpace.
  set_exponents(feSpace->getBasisFcts()->getDegree());    // Setting exponents.
  fill_struct_vec(uh);               // Filling vector of recovery structures.

  DOFVector<RecoveryStructure>::Iterator SV_it(struct_vec, USED_DOFS);

  // Solving local systems.
718
719
720
721
  for (SV_it.reset(); !SV_it.end(); ++SV_it) {
    if ((*SV_it).A) {
      TEST(Cholesky::solve((*SV_it).A, (*SV_it).rec_uh, (*SV_it).rec_uh))
	("There must be an error, matrix is not positive definite.\n");
722
    }
723
  }
724
725

  // define result vector
726
  DOFVector<double> *result = &rec_vec;
727
728
729
730

  result->set(0.0);

  DOFVector<double>::Iterator result_it(result, USED_DOFS);
731
  std::set<DegreeOfFreedom>::const_iterator setIterator;
732

733
734
735
736
737
// #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
//   DOFVector<int> neighbourSize(feSpace, "neighbourSize");
//   neighbourSize.set(0);
// #endif

738
739
740
741
742
743
744
745
746
  for (SV_it.reset(), result_it.reset(); !result_it.end(); ++SV_it, ++result_it) {
    if ((*SV_it).rec_uh) {
      *result_it = (*(*SV_it).rec_uh)[0];
    } else {
      if ((*SV_it).neighbors) {
	for (setIterator  = (*SV_it).neighbors->begin();
	     setIterator != (*SV_it).neighbors->end();
	     ++setIterator) {
	  for (int i = 0; i < n_monomials; i++)
747
748
	    *result_it += (*(*struct_vec)[*setIterator].rec_uh)[i] *
	      (*(*matrix_fcts)[0][i])(*(*SV_it).coords, *(*struct_vec)[*setIterator].coords);
749
	}
750
751
752
// #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
// 	neighbourSize[result_it.getDOFIndex()] = (*SV_it).neighbors->size();
// #else
753
	*result_it /= (*SV_it).neighbors->size();
754
755
756
757
	//#endif
      } else {
	*result_it = 0.0;
      }
758
    }
759
  }
760
761
762
763
764
765
766
767
768

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
//   MeshDistributor::globalMeshDistributor->synchAddVector(rec_vec);
//   MeshDistributor::globalMeshDistributor->synchAddVector(neighbourSize);

//   for (result_it.reset(); !result_it.end(); ++result_it)
//     if (neighbourSize[result_it.getDOFIndex()] > 0) 
//       *result_it /= neighbourSize[result_it.getDOFIndex()];

769
  Parallel::MeshDistributor::globalMeshDistributor->synchVector(rec_vec);
770
#endif
771
772
}

773

774
775
776
DOFVector<double>*
Recovery::recoveryUh(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
{
777
  FUNCNAME("Recovery::recoveryUh()");
778
779
780
781
782
783
784
785
786
787
788
789

  clear();

  gradient = false;

  set_feSpace(fe_space);                                  // Setting feSpace.
  set_exponents(feSpace->getBasisFcts()->getDegree());    // Setting exponents.
  fill_struct_vec(uh);               // Filling vector of recovery structures.

  DOFVector<RecoveryStructure>::Iterator SV_it(struct_vec, USED_DOFS);

  // Solving local systems.
790
791
792
793
  for (SV_it.reset(); !SV_it.end(); ++SV_it) {
    if ((*SV_it).A) {
      TEST(Cholesky::solve((*SV_it).A, (*SV_it).rec_uh, (*SV_it).rec_uh))
	("There must be an error, matrix is not positive definite.\n");
794
    }
795
796
  }
  
797
  // define result vector
798
  static DOFVector<double> *vec = nullptr;// TODO: REMOVE STATIC
799
  DOFVector<double> *result = nullptr;
800
801

  // Allocate memory for result vector
802
  if (vec && vec->getFeSpace() != feSpace) {
803
    delete vec;
804
    vec = nullptr;
805
806
  }

807
  if (!vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
808
    vec = new DOFVector<double>(feSpace, "gradient");
809
810
811
812
813
  result = vec;

  result->set(0.0);

  DOFVector<double>::Iterator result_it(result, USED_DOFS);
814
  std::set<DegreeOfFreedom>::const_iterator setIterator;
815
816

  for (SV_it.reset(), result_it.reset(); !result_it.end();
817
818
819
820
821
822
823
824
825
826
827
828
       ++SV_it, ++result_it) {
    if ((*SV_it).rec_uh) {
      *result_it = (*(*SV_it).rec_uh)[0];
    } else {
      if ((*SV_it).neighbors) {
	for (setIterator  = (*SV_it).neighbors->begin();
	     setIterator != (*SV_it).neighbors->end();
	     ++setIterator) {
	  for (int i = 0; i < n_monomials; i++)
	    *result_it = *result_it + (*(*struct_vec)[*setIterator].rec_uh)[i] *
	      (*(*matrix_fcts)[0][i])(*(*SV_it).coords,
				      *(*struct_vec)[*setIterator].coords);
829
	}
830
831
832
	*result_it /= (*SV_it).neighbors->size();
      } else
	*result_it=0.0;
833
    }
834
  }
835
836
837
838

  return result;
}

839

840
841
842
843
844
845
DOFVector<WorldVector<double> >*
Recovery::recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space,
		   AbstractFunction<double, WorldVector<double> > *f_vec,
		   AbstractFunction<double, double> *f_scal,
		   DOFVector<double> *aux_vec)
{
846
  FUNCNAME_DBG("Recovery::recovery()");
847
848
849
850
851
852
853
854
855
856
857
858

  clear();

  gradient = true;

  set_feSpace(fe_space);                                  // Setting feSpace.
  set_exponents(feSpace->getBasisFcts()->getDegree());    // Setting exponents.
  fill_struct_vec(uh, f_vec, f_scal, aux_vec);  // Filling vec. of rec. struct.

  DOFVector<RecoveryStructure>::Iterator SV_it(struct_vec, USED_DOFS);

  // Solving local systems.
Thomas Witkowski's avatar
Thomas Witkowski committed
859
860
  for (SV_it.reset(); !SV_it.end(); ++SV_it) {
    if ((*SV_it).A) {
861
862
863
864
      DBG_VAR(int error =)
	Cholesky::solve((*SV_it).A, 
			(*SV_it).rec_grdUh,
			(*SV_it).rec_grdUh);
Thomas Witkowski's avatar
Thomas Witkowski committed
865
866
      TEST_EXIT_DBG(error)
	("There must be some error, matrix is not positive definite.\n");
867
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
868
  }
869
870

  // define result vector
871
  static DOFVector<WorldVector<double> > *vec = nullptr;// TODO: REMOVE STATIC
872
  DOFVector<WorldVector<double> > *result = nullptr;
873
874

  // Allocate memory for result vector
875
  if (vec && vec->getFeSpace() != feSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
876
    delete vec;
877
    vec = nullptr;
Thomas Witkowski's avatar
Thomas Witkowski committed
878
879
  }
  
880
  if (!vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
881
    vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
Thomas Witkowski's avatar
Thomas Witkowski committed
882

883
884
885
886
887
  result = vec;

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

  DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS);
888
  std::set<DegreeOfFreedom>::const_iterator setIterator;
889

Thomas Witkowski's avatar
Thomas Witkowski committed
890
891
892
893
894
895
896
897
898
899
900
901
902
  for (SV_it.reset(), grdIt.reset(); !grdIt.end(); ++SV_it, ++grdIt) {
    if ((*SV_it).rec_grdUh) {
      *grdIt = (*(*SV_it).rec_grdUh)[0];
    } else {
      for (setIterator  = (*SV_it).neighbors->begin();
	   setIterator != (*SV_it).neighbors->end();
	   ++setIterator) {
	for (int i = 0; i < n_monomials; i++)
	  *grdIt = *grdIt + (*(*struct_vec)[*setIterator].rec_grdUh)[i] *
	    (*(*matrix_fcts)[0][i])(*(*SV_it).coords,
				    *(*struct_vec)[*setIterator].coords);
      }
      *grdIt = *grdIt * (1.0 / (*SV_it).neighbors->size());
903
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
904
  }
905
906
907
908

  return result;
}

909

910
911
912
913
914
915
DOFVector<WorldVector<double> >*
Recovery::recovery(DOFVector<double> *uh,
		   AbstractFunction<double, WorldVector<double> > *f_vec,
		   AbstractFunction<double, double> *f_scal,
		   DOFVector<double> *aux_vec)
{
916
  FUNCNAME_DBG("Recovery::simpleAveraging()");
917

918
  TEST_EXIT_DBG(!(f_vec && f_scal))("Only one diffusion function, please!\n");
919

920
  const FiniteElemSpace *fe_space = uh->getFeSpace();
921
922

  // define result vector
923
  static DOFVector<WorldVector<double> > *vec = nullptr;// TODO: REMOVE STATIC
924
  DOFVector<WorldVector<double> > *result = nullptr;
925
926

  // Allocate memory for result vector
927
  if (vec && vec->getFeSpace() != fe_space) {
Thomas Witkowski's avatar
Thomas Witkowski committed
928
    delete vec;
929
    vec = nullptr;
Thomas Witkowski's avatar
Thomas Witkowski committed
930
  }
931

932
933
  if (!vec) 
    vec = new DOFVector<WorldVector<double> >(fe_space, "gradient");  
934

935
  result = vec;
936
937
938
939
940
941
  result->set(WorldVector<double>(DEFAULT_VALUE, 0.0));

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

  Mesh *mesh = fe_space->getMesh();
Thomas Witkowski's avatar
Thomas Witkowski committed
942
  int dim = mesh->getDim();
943
  const BasisFunction *basFcts = fe_space->getBasisFcts();
Thomas Witkowski's avatar
Thomas Witkowski committed
944
  DOFAdmin *admin = fe_space->getAdmin();
945
  int numPreDofs = admin->getNumberOfPreDofs(0);
946
947
948
949
950
  DimVec<double> bary(