RecoveryEstimator.cc 6.44 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
14
15
#include "RecoveryEstimator.h"
#include "Parameters.h"

16
17
namespace AMDiS {

18
  RecoveryEstimator::RecoveryEstimator(std::string name, DOFVector<double> *uh_, int r) 
19
    : Estimator(name, r), 
20
21
      uh(uh_), 
      relative(0), 
22
23
      C(1.0), 
      method(0),
24
      addEstimationToOld(false),
25
26
27
28
29
30
31
      feSpace(NULL), 
      f_vec(NULL), 
      f_scal(NULL), 
      aux_vec(NULL), 
      rec_struct(NULL)
  {
    FUNCNAME("RecoveryEstimator::constructor()");
32
    
33
34
    GET_PARAMETER(0, name + "->rec method", "%d", &method); // 0, 1, or 2 (see Recovery.h)
    GET_PARAMETER(0, name + "->rel error", "%d", &relative); // 0 or 1
35
36
37
38
39
    
    GET_PARAMETER(0, name + "->C", "%f", &C);
    C = C > 1e-25 ? sqr(C) : 1.0;
    
    if (norm == H1_NORM) {
40
      feSpace = uh_->getFeSpace();
41
42
43
44
45
46
47
      degree = feSpace->getBasisFcts()->getDegree();
      
      if (degree <= 2 && C != 1.0) {
	WARNING("Recovery estimators in the H1_NORM usually very good for linear and quadratic finite element; normally you do not need to overwrite the default value of C\n");
	WAIT;
      }
    } else {
48
      degree = uh_->getFeSpace()->getBasisFcts()->getDegree() + 1;    
49
      feSpace = FiniteElemSpace::provideFeSpace(NULL,
50
51
						Lagrange::getLagrange(uh_->getFeSpace()->getMesh()->getDim(), degree),
						uh_->getFeSpace()->getMesh(),
52
53
54
55
56
57
58
59
						name + "->feSpace");

      if (method == 2) {
	ERROR("Simple averaging only for the H1_NORM; using SPR instead\n");
	WAIT;
	method = 0;
      }
    }
60

61
62
    if (method == 2 && degree !=1) {
      ERROR("Simple averaging only for linear elements; using SPR instead\n");
63
64
65
66
      WAIT;
      method = 0;
    }

67
    rec_struct = new Recovery(norm, method);
68
69
  }

70
  void RecoveryEstimator::init(double ts)
71
  {
72
73
74
    FUNCNAME("RecoveryEstimator::init()");

    basFcts = uh->getFeSpace()->getBasisFcts();
75
    int dim = mesh->getDim();
76
77
    h1Norm2 = 0.0;

78
    if (norm == H1_NORM) {    // sets recovery gradient.
79
      if (method == 2)
80
	rec_grd = rec_struct->recovery(uh, f_vec, f_scal, aux_vec);
81
      else
82
	rec_grd = rec_struct->recovery(uh, feSpace, f_vec, f_scal, aux_vec);
83

84
      rec_basFcts = rec_grd->getFeSpace()->getBasisFcts();
85
    } else {                 // sets higher-order recovery solution.
86
      rec_uh = rec_struct->recoveryUh(uh, feSpace);
87
      rec_basFcts = rec_uh->getFeSpace()->getBasisFcts();
88
89
    }

90
    int deg = 2 * std::max(basFcts->getDegree(), rec_basFcts->getDegree());
91
92
    quad = Quadrature::provideQuadrature(dim, deg);
    nPoints = quad->getNumPoints();
93
    
94
95
//  WorldVector<double> quad_pt;
    grdAtQP = new WorldVector<double>[nPoints];
96

97
98
99
    recoveryGrdAtQP = mtl::dense_vector<WorldVector<double> >(nPoints);
    uhAtQP = mtl::dense_vector<double>(nPoints);
    recoveryUhAtQP = mtl::dense_vector<double>(nPoints);
100
    
101
102
103
    quadFast = FastQuadrature::provideFastQuadrature(basFcts, *quad, INIT_PHI | INIT_GRD_PHI);
    rec_quadFast = FastQuadrature::provideFastQuadrature(rec_basFcts, *quad, INIT_PHI | INIT_GRD_PHI);

104
    // clear error indicators
105
    if(!addEstimationToOld) {
106
107
108
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
109
110
111
      elInfo->getElement()->setEstimation(0.0, row);
      elInfo = stack.traverseNext(elInfo);
    }
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
    }

    est_sum = 0.0;
    est_max = 0.0;
    est_t_sum = 0.0;
    est_t_max = 0.0;

    traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;

  }

  void RecoveryEstimator::exit(bool output)
  {
    FUNCNAME("RecoveryEstimator::exit()");

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double send_est_sum = est_sum;
    double send_est_max = est_max;
130

131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
    MPI::COMM_WORLD.Allreduce(&send_est_sum, &est_sum, 1, MPI_DOUBLE, MPI_SUM);
    MPI::COMM_WORLD.Allreduce(&send_est_max, &est_max, 1, MPI_DOUBLE, MPI_MAX);
#endif

    // Computing relative errors
    if (relative) {
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      while (elInfo) {
	double estEl = elInfo->getElement()->getEstimation(row);
	estEl /= h1Norm2;
	elInfo->getElement()->setEstimation(estEl, row);
	elInfo = stack.traverseNext(elInfo);
      }
      
      est_max /= h1Norm2;
      est_sum /= h1Norm2;
    }

    est_sum = sqrt(est_sum);

    delete [] grdAtQP;
153
    
154
155
156
157
158
159
160
161
162
    if (output) {
      MSG("estimate for component %d = %.8e\n", row, est_sum);
    }
  }

  void RecoveryEstimator::estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo)
  {    
    FUNCNAME("RecoveryEstimator::estimateElement()");

163
164
165
      Element *el = elInfo->getElement();
      double det = elInfo->getDet();
      double errEl = 0.0;
166
167
      double estEl = el->getEstimation(row);
      int dow = Global::getGeo(WORLD);
168
169
170

      if (norm == H1_NORM) {
	// get gradient and recovery gradient at quadrature points
171
	uh->getGrdAtQPs(elInfo, NULL, quadFast, grdAtQP);
172
173
174
175
176
	rec_grd->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryGrdAtQP);
	if (f_scal) {
	  if (aux_vec)
	    aux_vec->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
	  else
177
	    uh->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
178
179
	}

180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
	// calc h1 error
	for (int i = 0; i < nPoints; i++) {
	  double  err2 = 0.0;
	  double fAtQP = 1.0;
	  if (f_scal)
	    fAtQP = (*f_scal)(uhAtQP[i]);
	  if (f_vec) {
	    elInfo->coordToWorld(quad->getLambda(i), quad_pt);
	    fAtQP = (*f_vec)(quad_pt);
	  }
	  
	  for (int j = 0; j < dow; j++)
	    err2 += sqr(recoveryGrdAtQP[i][j] - fAtQP * grdAtQP[i][j]);
	  errEl += quad->getWeight(i) * err2;
	}
      } else {
	// get vector and recovery vector at quadrature points
197
	uh->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
198
199
200
201
202
203
204
	rec_uh->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryUhAtQP);
	
	// calc l2 error
	for (int i = 0; i < nPoints; i++)
	  errEl += quad->getWeight(i) * sqr(recoveryUhAtQP[i] - uhAtQP[i]);
      }
      
205
      estEl += C * det * errEl;
206
207
      el->setEstimation(estEl, row);
      est_sum += estEl;
208
      est_max = std::max(est_max, estEl);
209

210
      if (relative) {	
211
212
213
214
215
216
217
218
219
220
221
222
	double normEl = 0.0;
	
	if (norm == H1_NORM) {
	  for (int i = 0; i < nPoints; i++) {
	    double norm2 = 0.0;
	    for (int j = 0; j < dow; j++)
	      norm2 += sqr(recoveryGrdAtQP[i][j]);
	    normEl += quad->getWeight(i) * norm2;
	  }
	} else {
	  for (int i = 0; i < nPoints; i++)
	    normEl += quad->getWeight(i) * sqr(recoveryUhAtQP[i]);
223
224
	}

225
226
227
	h1Norm2 += det * normEl;
      }
  }
228
229

}