Recovery.cc 29.3 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 "Recovery.h"
14
15
16
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/MeshDistributor.h"
#endif
17

18
19
RecoveryStructure& RecoveryStructure::operator=(const RecoveryStructure& rhs) 
{ 
20
  if (rhs.coords) {
21
22
    if (!coords) 
      coords = new WorldVector<double>;    
Thomas Witkowski's avatar
Thomas Witkowski committed
23
24
    *coords = *rhs.coords;
  } else {
25
    if (coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
26
      delete coords;
Thomas Witkowski's avatar
Thomas Witkowski committed
27
      coords = NULL;
28
29
30
31
32
    }
  }

  if (rhs.A) {
    if (!A) 
Thomas Witkowski's avatar
Thomas Witkowski committed
33
      A = new Matrix<double>(rhs.A->getNumRows(), rhs.A->getNumCols());
Thomas Witkowski's avatar
Thomas Witkowski committed
34
35
    *A = *rhs.A;
  } else {
36
    if (A) {
Thomas Witkowski's avatar
Thomas Witkowski committed
37
      delete A;
Thomas Witkowski's avatar
Thomas Witkowski committed
38
      A = NULL;
39
40
41
42
43
    }
  }

  if (rhs.rec_uh) {
    if (!rec_uh)
Thomas Witkowski's avatar
Thomas Witkowski committed
44
      rec_uh = new Vector<double>(rhs.rec_uh->getSize());
Thomas Witkowski's avatar
Thomas Witkowski committed
45
46
    *rec_uh = *rhs.rec_uh;
  } else {
47
    if (rec_uh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
48
      delete rec_uh;
49
      rec_uh = NULL;
50
51
52
53
54
    }
  }

  if (rhs.rec_grdUh) {
    if (!rec_grdUh)
Thomas Witkowski's avatar
Thomas Witkowski committed
55
      rec_grdUh = new Vector<WorldVector<double> >(rhs.rec_grdUh->getSize());
Thomas Witkowski's avatar
Thomas Witkowski committed
56
57
    *rec_grdUh = *rhs.rec_grdUh;
  } else {
58
    if (rec_grdUh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
59
      delete rec_grdUh;
Thomas Witkowski's avatar
Thomas Witkowski committed
60
      rec_grdUh = NULL;
61
62
63
64
65
    }
  }

  if (rhs.neighbors) {
    if (!neighbors)
Thomas Witkowski's avatar
Thomas Witkowski committed
66
      neighbors = new std::set<DegreeOfFreedom>;
67
    *neighbors = *rhs.neighbors;
Thomas Witkowski's avatar
Thomas Witkowski committed
68
  } else {
69
    if (neighbors) {
Thomas Witkowski's avatar
Thomas Witkowski committed
70
      delete neighbors;
Thomas Witkowski's avatar
Thomas Witkowski committed
71
      neighbors = NULL;
72
73
74
75
    }
  }

  return *this;
76
}
77
78
79
80


void RecoveryStructure::print()
{
81
  FUNCNAME("RecoveryStructure::print()");
82

83
  std::cout << std::endl;
84
85

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

88
89
90
91
92
93
94
95
96
97
98
99
100
101
  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;
    }
102

103
104
105
106
107
108
109
110
    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;
111
    }
112
113
114

    if (neighbors) {
      MSG("Printing neighbors vertices\n\n");
115
116

      MSG("Number of neighbors: ");
117
      std::cout << neighbors->size() << std::endl << std::endl;
118
      
119
      MSG("List of neighbors: ");
120
      std::set<DegreeOfFreedom>::const_iterator setIterator;
121
122
      
      for (setIterator = neighbors->begin(); setIterator != neighbors->end();
123
	   ++setIterator)
124
	std::cout << " " << *setIterator;
125
      
126
      std::cout << std::endl << std::endl;
127
    }
128
129
  } else {
    MSG("Boundary vertex or not a vertex node: printing interior neighbors\n\n");
130

131
132
133
134
135
136
137
138
139
    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;
140

141
142
143
144
    std::cout << std::endl << std::endl;
  }
  
  WAIT;
145
146
147
148
149
}


void Recovery::set_feSpace(const FiniteElemSpace *fe_space)
{
150
  FUNCNAME("Recovery::set_feSpace()");
151

152
153
154
155
  if (!feSpace || feSpace != fe_space) {
    if (struct_vec) {
      delete struct_vec;
      struct_vec = NULL;
156
157
    }

158
159
160
161
162
    feSpace = fe_space;
    
    // create new structure vector
    struct_vec = new DOFVector<RecoveryStructure>(feSpace, "struct vec");
  }
163
164
}

165

166
167
int Recovery::set_exponents(int degree)
{
168
  FUNCNAME("Recovery::set_exponents()");
169
170
171
172
173

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

  // Computing number of monomials.
174
175
176
177
  if (dow > 1) {
    number_monomials *= (degree + 2);
    number_monomials /= 2;
  }
178

179
180
181
182
  if (dow == 3) {
    number_monomials *= (degree + 3);
    number_monomials /= 3;
  }
183
184

  // Allocating memory.
185
186
187
188
189
190
191
192
193
  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);
  }
194
195

  // Setting vector of exponents.
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
  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;
221

222
223
224
  default:
    ERROR_EXIT("Which dimension have your world???\n");
  }
225

226
227
  TEST_EXIT_DBG(count == n_monomials)("There must be an error!\n");
  
228
229
230
  // Setting matrix of monomials.
  WorldVector<int> sum;

231
232
233
234
235
236
  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);
    }
237
238
239
240

  return n_monomials;
}

241

242
243
244
245
246
247
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)
{
248
  FUNCNAME("Recovery::compute_integrals()");
249

250
  TEST_EXIT_DBG(!(f_vec && f_scal))("Only one diffusion function, please!\n");
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265

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

266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
  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();
    }
281

282
283
    Quadrature *quad = Quadrature::
      provideQuadrature(Global::getGeo(WORLD),
284
			  (*matrix_fcts)[0][i]->getDegree() + deg_f);
285
    int n_points = quad->getNumPoints();
286
    mtl::dense_vector<double> uhAtQP(n_points);
287
288
289
290
291
292
293
294
    
    // Computing contributions to right hand side.
    if (gradient) {    // For gradient recovery.
      double fAtQP = 1.0;
      if (f_scal) {
	if (aux_vec)
	  aux_vec->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
	else
295
	  uh->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
296
297
298
      }

      // Get gradient at quadrature points
299
      mtl::dense_vector<WorldVector<double> > grdAtQP(n_points);
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
      uh->getGrdAtQPs(elInfo, quad, NULL, grdAtQP);
      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
      uh->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
      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();
324
    }
325
  }
326
327
}

328

329
330
331
332
333
334
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)
{
335
  FUNCNAME("Recovery::compute_sums()");
336

337
338
  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");
339
340
341
342

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

346
347
348
349
350
351
352
353
354
355
  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;
    }
356

357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
    // Computing contributions to right hand side.
    double fAtQP = 1.0;
    if (f_scal) {
      if (aux_vec)
	aux_vec->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
      else
	uh->getVecAtQPs(elInfo, quad, NULL, uhAtQP);
    }

    // Get gradient at quadrature points
    uh->getGrdAtQPs(elInfo, quad, NULL, grdAtQP);
    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);
373
      if (f_scal)
374
	fAtQP = (*f_scal)(uhAtQP[k]);
375

376
377
      vec_sum = vec_sum + grdAtQP[k] * fAtQP
	* (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
378
    }
379
380
    (*rec_struct->rec_grdUh)[i] = (*rec_struct->rec_grdUh)[i] + vec_sum;
  }
381
382
}

383

384
void Recovery::compute_node_sums(DOFVector<double> *uh, ElInfo *elInfo,
385
				 RecoveryStructure *rec_struct, DimVec<int> preDofs,
386
387
				 int n_vertices, int n_edges, int n_faces)
{
388
  FUNCNAME("Recovery::compute_sums()");
389

390
  TEST_EXIT_DBG(!gradient)
391
    ("SPR of flux or gradient need computing interior sums\n");
392
  TEST_EXIT_DBG(feSpace->getMesh()->getDim() == 1)
393
394
395
    ("At the moment only for linear finite elements.\n");

  WorldVector<double> node;  // For world coordinates at nodes.
396
  const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
397
398
  ElementVector uh_loc(n_vertices);
  uh->getLocalVector(elInfo->getElement(), uh_loc);
399

400
401
  for (int l = 0; l < n_vertices; l++) {
    // Computing contributions of vertex nodes
402
    if (rec_struct->neighbors->insert(dof[l][preDofs[VERTEX]]).second) {
403
404
405
406
407
408
409
410
411
412
413
      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);
      }
414
    }
415
  }
416
417
}

418

419
420
void Recovery::compute_sums_linear(DOFVector<double> *uh, ElInfo *elInfo,
				   RecoveryStructure *rec_struct,
421
				   int vertex, DimVec<int> preDofs,
422
423
				   int n_vertices)
{
424
  FUNCNAME("Recovery::compute_sums_linear()");
425

426
  TEST_EXIT_DBG(!gradient)
427
428
429
    ("SPR of flux or gradient need computing interior sums\n");

  WorldVector<double> node;     // For world coordinates at nodes.
430
  const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
431
432
  ElementVector uh_loc(n_vertices);
  uh->getLocalVector(elInfo->getElement(), uh_loc);
433

434
435
  for (int l = 0;  l < n_vertices; l++) {
    // Computing contributions of vertex nodes
436
    DegreeOfFreedom k = dof[l][preDofs[VERTEX]];
437
438
439
440
441
442
443
444
445
446
447
448
    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);
      }
449
    }
450
  }
451

452
453
  if (vertex > 1 && elInfo->getNeighbour(vertex)) {
    int oppVertex = elInfo->getOppVertex(vertex);
454
    DegreeOfFreedom k = elInfo->getNeighbour(vertex)->getDof(oppVertex)[preDofs[VERTEX]];
455
456
457
458
459
460
461
462
463
464
    
    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);
465
      }
466
467
    }
  }  
468
469
}

470

471
472
473
474
475
void Recovery::fill_struct_vec(DOFVector<double> *uh,
			       AbstractFunction<double, WorldVector<double> > *f_vec,
			       AbstractFunction<double, double> *f_scal,
			       DOFVector<double> *aux_vec)
{
476
  FUNCNAME("Recovery::fill_struct_vec()");
477
478
479

  // Information on the mesh.
  Mesh *mesh = feSpace->getMesh();
480
  int dim = mesh->getDim();
481
482
483

  // Geometric information.
  int n_vertices = Global::getGeo(VERTEX, dim);
484
485
  int n_edges = Global::getGeo(EDGE, dim);
  int n_faces = Global::getGeo(FACE, dim);
486
487
488

  // Information concerning the finite element space.
  const BasisFunction *basis_fcts = feSpace->getBasisFcts();
489
  DimVec<int> *nDOFs = basis_fcts->getNumberOfDofs();
490
491
492

  // Information from DOFAdmin.
  const DOFAdmin *admin = feSpace->getAdmin();
493
  DimVec<int> preDofs = admin->getNumberOfPreDofs();
494
495
496
497
498
499
500
501
502
503
504
505
506

  // 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();
  Quadrature *quad = NULL;
  if (gradient && !method)
    quad = Quadrature::provideQuadrature(Global::getGeo(WORLD), degree);

  DimVec<int> pre_dofs(dim, NO_INIT);
  if (!gradient)
507
    pre_dofs = uh->getFeSpace()->getAdmin()->getNumberOfPreDofs();
508
  
509
  // Variables for traversing the mesh.
510
  Flag fill_flag = 
511
512
513
514
    Mesh::CALL_LEAF_EL |
    Mesh::FILL_COORDS |
    Mesh::FILL_BOUND |
    Mesh::FILL_DET |
515
    Mesh::FILL_GRD_LAMBDA;
516

517
  if (degree == 2 && dim > 1)
518
519
    fill_flag |= Mesh::FILL_NEIGH | Mesh::FILL_OPP_COORDS;

520
521
522
523
524
525
526
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  DofContainer boundaryDofs;
  DofContainerSet boundaryDofSet;
  MeshDistributor::globalMeshDistributor->getAllBoundaryDofs(feSpace, 0, boundaryDofs);
  boundaryDofSet.insert(boundaryDofs.begin(), boundaryDofs.end());
#endif

527
  TraverseStack stack;
528
  ElInfo *el_info = stack.traverseFirst(mesh, -1, fill_flag);
529

530
  while (el_info) {    // traversing the mesh.    
531
    const DegreeOfFreedom **dof = el_info->getElement()->getDof();
532
533
534
    
    int n_neighbors = 0;     // counting interior vertices of element
    for (int i = 0; i < n_vertices; i++) {
535
      DegreeOfFreedom k = dof[i][preDofs[VERTEX]];
536
537
538
539
540
541
542
#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)
543
544
545
546
547
548
549
	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.
550
551
552
553
554
555
#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

556
      DegreeOfFreedom k = dof[i][preDofs[VERTEX]];
557
558
559
560
      
      // Setting world coordinates of node.
      if (!(*struct_vec)[k].coords) {	
	(*struct_vec)[k].coords  = new WorldVector<double>;
561
	*(*struct_vec)[k].coords = el_info->getCoord(i);
562
      }
563

564
      if (isInterior) {	
565
566
567
568
569
570
571
572
573
574
575
576
577
	// 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;
	  }
578
579
	}

580
581
	// Computing the integrals.
	if (method)
582
	  compute_integrals(uh, el_info, &(*struct_vec)[k],
583
584
			    f_vec, f_scal, aux_vec);
	else if (gradient)
585
	  compute_interior_sums(uh, el_info, &(*struct_vec)[k], quad,
586
587
588
589
590
591
				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)
592
	    compute_sums_linear(uh, el_info, &(*struct_vec)[k],
593
594
				i, pre_dofs, n_vertices);
	  else
595
	    compute_node_sums(uh, el_info, &(*struct_vec)[k],
596
			      pre_dofs, n_vertices, n_edges, n_faces);
597
	}
598
599
600
601
602
603
604
605
      } 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]);
      }
    }
606

607
    int n_count = n_vertices;
608

609
610
611
612
    if (dim > 1) {     // Handling nodes on edges.
      for (int i = 0; i < n_edges; i++) {
	bool isEdgeInterior = el_info->getBoundary(EDGE, i) == INTERIOR;

613
	for (int j = 0; j < (*nDOFs)[EDGE]; j++) {
614
	  DegreeOfFreedom k = dof[n_vertices + i][preDofs[EDGE] + j];
615
616

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

643
644
645
    if (dim == 3)     // Handling nodes on faces.
      for (int i = 0; i < n_faces; i++)
	for (int j = 0; j < (*nDOFs)[FACE]; j++) {
646
	  DegreeOfFreedom k = dof[n_vertices+n_edges + i][preDofs[FACE] + j];
647
648
649
	  
	  if (!(*struct_vec)[k].coords) {
	    // Setting world coordinates of node.
650
	    el_info->coordToWorld(*basis_fcts->getCoords(n_count),
651
652
653
654
655
656
657
				  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>;
	    
658
	    if (el_info->getBoundary(FACE, i) == INTERIOR)
659
660
	      for (int m = 0; m < 3; m++) {
		int l = Global::getReferenceElement(dim)->getVertexOfPosition(FACE, i, m);
661
		if (el_info->getBoundary(VERTEX, l) == INTERIOR)
662
		  (*struct_vec)[k].neighbors->insert(dof[l][preDofs[VERTEX]]);
663
664
665
666
667
668
669
670
671
672
673
	      }
	    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++) {
674
	DegreeOfFreedom k = dof[n_vertices+n_edges+n_faces][preDofs[CENTER] + j];
675
676

	// Setting world coordinates of node.
677
	el_info->coordToWorld(*basis_fcts->getCoords(n_count), coordinates);
678
679
680
681
682
683
684
685
686
687
688
689
	(*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++;
      }
    
690
    el_info = stack.traverseNext(el_info);
691
  }
692
693
}

694

695
696
void Recovery::recoveryUh(DOFVector<double> *uh, DOFVector<double> &rec_vec)
{
697
  FUNCNAME("Recovery::recoveryUh()");
698
699
700
701
702

  clear();

  gradient = false;

703
  const FiniteElemSpace *fe_space = rec_vec.getFeSpace();
704
705
706
707
708
709
710
  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.
711
712
713
714
  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");
715
    }
716
  }
717
718

  // define result vector
719
  DOFVector<double> *result = &rec_vec;
720
721
722
723

  result->set(0.0);

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

726
727
728
729
730
// #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
//   DOFVector<int> neighbourSize(feSpace, "neighbourSize");
//   neighbourSize.set(0);
// #endif

731
732
733
734
735
736
737
738
739
  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++)
740
741
	    *result_it += (*(*struct_vec)[*setIterator].rec_uh)[i] *
	      (*(*matrix_fcts)[0][i])(*(*SV_it).coords, *(*struct_vec)[*setIterator].coords);
742
	}
743
744
745
// #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
// 	neighbourSize[result_it.getDOFIndex()] = (*SV_it).neighbors->size();
// #else
746
	*result_it /= (*SV_it).neighbors->size();
747
748
749
750
	//#endif
      } else {
	*result_it = 0.0;
      }
751
    }
752
  }
753
754
755
756
757
758
759
760
761
762
763

#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()];

  MeshDistributor::globalMeshDistributor->synchVector(rec_vec);
#endif
764
765
}

766

767
768
769
DOFVector<double>*
Recovery::recoveryUh(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
{
770
  FUNCNAME("Recovery::recoveryUh()");
771
772
773
774
775
776
777
778
779
780
781
782

  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.
783
784
785
786
  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");
787
    }
788
789
  }
  
790
791
  // define result vector
  static DOFVector<double> *vec = NULL;
792
  DOFVector<double> *result = NULL;
793
794

  // Allocate memory for result vector
795
  if (vec && vec->getFeSpace() != feSpace) {
796
797
798
799
    delete vec;
    vec = NULL;
  }

800
  if (!vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
801
    vec = new DOFVector<double>(feSpace, "gradient");
802
803
804
805
806
  result = vec;

  result->set(0.0);

  DOFVector<double>::Iterator result_it(result, USED_DOFS);
807
  std::set<DegreeOfFreedom>::const_iterator setIterator;
808
809

  for (SV_it.reset(), result_it.reset(); !result_it.end();
810
811
812
813
814
815
816
817
818
819
820
821
       ++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);
822
	}
823
824
825
	*result_it /= (*SV_it).neighbors->size();
      } else
	*result_it=0.0;
826
    }
827
  }
828
829
830
831

  return result;
}

832

833
834
835
836
837
838
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)
{
Thomas Witkowski's avatar
Thomas Witkowski committed
839
  FUNCNAME("Recovery::recovery()");
840
841
842
843
844
845
846
847
848
849
850
851

  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
852
853
  for (SV_it.reset(); !SV_it.end(); ++SV_it) {
    if ((*SV_it).A) {
854
855
856
857
      DBG_VAR(int error =)
	Cholesky::solve((*SV_it).A, 
			(*SV_it).rec_grdUh,
			(*SV_it).rec_grdUh);
Thomas Witkowski's avatar
Thomas Witkowski committed
858
859
      TEST_EXIT_DBG(error)
	("There must be some error, matrix is not positive definite.\n");
860
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
861
  }
862
863
864

  // define result vector
  static DOFVector<WorldVector<double> > *vec = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
865
  DOFVector<WorldVector<double> > *result = NULL;
866
867

  // Allocate memory for result vector
868
  if (vec && vec->getFeSpace() != feSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
869
    delete vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
870
871
872
    vec = NULL;
  }
  
873
  if (!vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
874
    vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
Thomas Witkowski's avatar
Thomas Witkowski committed
875

876
877
878
879
880
  result = vec;

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
883
884
885
886
887
888
889
890
891
892
893
894
895
  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());
896
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
897
  }
898
899
900
901

  return result;
}

902

903
904
905
906
907
908
DOFVector<WorldVector<double> >*
Recovery::recovery(DOFVector<double> *uh,
		   AbstractFunction<double, WorldVector<double> > *f_vec,
		   AbstractFunction<double, double> *f_scal,
		   DOFVector<double> *aux_vec)
{
Thomas Witkowski's avatar
Thomas Witkowski committed
909
  FUNCNAME("Recovery::simpleAveraging()");
910

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

913
  const FiniteElemSpace *fe_space = uh->getFeSpace();
914
915
916

  // define result vector
  static DOFVector<WorldVector<double> > *vec = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
917
  DOFVector<WorldVector<double> > *result = NULL;
918
919

  // Allocate memory for result vector
920
  if (vec && vec->getFeSpace() != fe_space) {
Thomas Witkowski's avatar
Thomas Witkowski committed
921
    delete vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
922
923
    vec = NULL;
  }
924

925
926
  if (!vec) 
    vec = new DOFVector<WorldVector<double> >(fe_space, "gradient");  
927

928
  result = vec;
929
930
931
932
933
934
  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
935
  int dim = mesh->getDim();
936
  const BasisFunction *basFcts = fe_space->getBasisFcts();
Thomas Witkowski's avatar
Thomas Witkowski committed
937
  DOFAdmin *admin = fe_space->getAdmin();
938
  int numPreDofs = admin->getNumberOfPreDofs(0);
939
940
941
942
943
  DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
  WorldVector<double> barycenter;     // For world coordinates at barycenter

  // traverse mesh
  TraverseStack stack;
944
945
  Flag fillFlag =
    Mesh::CALL_LEAF_EL | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
946
947
  ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

948
949
  ElementVector localUh(basFcts->getNumber());

Thomas Witkowski's avatar
Thomas Witkowski committed
950
951
  while (elInfo) {
    double det = elInfo->getDet();
952
    const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
953
954
955
    uh->getLocalVector(elInfo->getElement(), localUh);


Thomas Witkowski's avatar
Thomas Witkowski committed
956
957
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    WorldVector<double> grd;
958
    basFcts->evalGrdUh(bary, grdLambda, localUh, grd);
Thomas Witkowski's avatar
Thomas Witkowski committed
959
960
961
962
    
    double fAtBary = 1.0;
    
    if (f_vec) {
963
      elInfo->coordToWorld(bary, barycenter);
Thomas Witkowski's avatar
Thomas Witkowski committed
964
965
966
967
968
      fAtBary = (*f_vec)(barycenter);
    }
    
    if (f_scal) {
      if (aux_vec)
969
	aux_vec->getLocalVector(elInfo->getElement(), localUh);
Thomas Witkowski's avatar
Thomas Witkowski committed
970
971
972
      
      fAtBary = basFcts->evalUh(bary, localUh);
      fAtBary = (*f_scal)(fAtBary);
973
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
974
    
Thomas Witkowski's avatar
Thomas Witkowski committed
975
    for (int i = 0; i < dim + 1; i++) {