RecoveryEstimator.cc 6.11 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
24
25
26
27
28
29
30
      C(1.0), 
      method(0),
      feSpace(NULL), 
      f_vec(NULL), 
      f_scal(NULL), 
      aux_vec(NULL), 
      rec_struct(NULL)
  {
    FUNCNAME("RecoveryEstimator::constructor()");
31
    
32
33
    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
34
35
36
37
38
    
    GET_PARAMETER(0, name + "->C", "%f", &C);
    C = C > 1e-25 ? sqr(C) : 1.0;
    
    if (norm == H1_NORM) {
39
      feSpace = uh_->getFeSpace();
40
41
42
43
44
45
46
      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 {
47
      degree = uh_->getFeSpace()->getBasisFcts()->getDegree() + 1;    
48
      feSpace = FiniteElemSpace::provideFeSpace(NULL,
49
50
						Lagrange::getLagrange(uh_->getFeSpace()->getMesh()->getDim(), degree),
						uh_->getFeSpace()->getMesh(),
51
52
53
54
55
56
57
58
						name + "->feSpace");

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

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

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

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

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

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

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

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

96
97
98
    recoveryGrdAtQP = mtl::dense_vector<WorldVector<double> >(nPoints);
    uhAtQP = mtl::dense_vector<double>(nPoints);
    recoveryUhAtQP = mtl::dense_vector<double>(nPoints);
99
    
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
    quadFast = FastQuadrature::provideFastQuadrature(basFcts, *quad, INIT_PHI | INIT_GRD_PHI);
    rec_quadFast = FastQuadrature::provideFastQuadrature(rec_basFcts, *quad, INIT_PHI | INIT_GRD_PHI);

    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;
119

120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
    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;
142
    
143
144
145
146
147
148
149
150
151
    if (output) {
      MSG("estimate for component %d = %.8e\n", row, est_sum);
    }
  }

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

152
153
154
      Element *el = elInfo->getElement();
      double det = elInfo->getDet();
      double errEl = 0.0;
155
      double estEl = 0.0;
156
      int dow = Global::getGeo(WORLD);
157
158
159

      if (norm == H1_NORM) {
	// get gradient and recovery gradient at quadrature points
160
	uh->getGrdAtQPs(elInfo, NULL, quadFast, grdAtQP);
161
162
163
164
165
	rec_grd->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryGrdAtQP);
	if (f_scal) {
	  if (aux_vec)
	    aux_vec->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
	  else
166
	    uh->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
167
168
	}

169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
	// 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
186
	uh->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
187
188
189
190
191
192
193
	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]);
      }
      
194
      estEl += C * det * errEl;
195
196
      el->setEstimation(estEl, row);
      est_sum += estEl;
197
      est_max = std::max(est_max, estEl);
198

199
      if (relative) {	
200
201
202
203
204
205
206
207
208
209
210
211
	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]);
212
213
	}

214
215
216
	h1Norm2 += det * normEl;
      }
  }
217
218

}