Hermes2D  2.0
global.cpp
1 // This file is part of Hermes2D.
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
15 
16 #include "global.h"
17 #include <algorithm>
18 #include "quadrature/quad_all.h"
19 #include "mesh.h"
20 #include "traverse.h"
21 #include "refmap.h"
22 #include "function/solution.h"
23 #include "quadrature/limit_order.h"
24 #include "integrals/h1.h"
25 #include "discrete_problem.h"
26 
27 namespace Hermes
28 {
29  namespace Hermes2D
30  {
31  class Transformable;
32 
33  template<typename Scalar>
34  double Global<Scalar>::get_l2_norm(Vector<Scalar>* vec)
35  {
36  Scalar val = 0;
37  for (unsigned int i = 0; i < vec->length(); i++)
38  {
39  Scalar inc = vec->get(i);
40  val = val + inc*conj(inc);
41  }
42  return sqrt(std::abs(val));
43  }
44 
45  template<typename Scalar>
46  double Global<Scalar>::calc_abs_error(MeshFunction<Scalar>* sln1, MeshFunction<Scalar>* sln2, int norm_type)
47  {
48  // sanity checks
49  if(sln1 == NULL) throw Hermes::Exceptions::Exception("sln1 is NULL in calc_abs_error().");
50  if(sln2 == NULL) throw Hermes::Exceptions::Exception("sln2 is NULL in calc_abs_error().");
51 
52  Quad2D* quad = &g_quad_2d_std;
53  sln1->set_quad_2d(quad);
54  sln2->set_quad_2d(quad);
55 
56  const Mesh* meshes[2] = { sln1->get_mesh(), sln2->get_mesh() };
57  Transformable* tr[2] = { sln1, sln2 };
58  Traverse trav(true);
59  trav.begin(2, meshes, tr);
60 
61  double error = 0.0;
62  Traverse::State* ee;
63  while ((ee = trav.get_next_state()) != NULL)
64  {
65  update_limit_table(ee->e[0]->get_mode());
66 
67  RefMap* ru = sln1->get_refmap();
68  RefMap* rv = sln2->get_refmap();
69  switch (norm_type)
70  {
71  case HERMES_L2_NORM:
72  error += error_fn_l2(sln1, sln2, ru, rv);
73  break;
74  case HERMES_H1_NORM:
75  error += error_fn_h1(sln1, sln2, ru, rv);
76  break;
77  case HERMES_HCURL_NORM:
78  error += error_fn_hc(sln1, sln2, ru, rv);
79  break;
80  case HERMES_HDIV_NORM:
81  error += error_fn_hdiv(sln1, sln2, ru, rv);
82  break;
83  default: throw Hermes::Exceptions::Exception("Unknown norm in calc_error().");
84  }
85  }
86  trav.finish();
87  return sqrt(error);
88  }
89 
90  template<typename Scalar>
91  double Global<Scalar>::calc_rel_error(MeshFunction<Scalar>* sln, MeshFunction<Scalar>* ref_sln, int norm_type)
92  {
93  double error = calc_abs_error(sln, ref_sln, norm_type);
94  double norm = calc_norm(ref_sln, norm_type);
95 
96  return error/norm;
97  }
98 
99  template<typename Scalar>
100  double Global<Scalar>::calc_norm(MeshFunction<Scalar>* sln, int norm_type)
101  {
102  Quad2D* quad = &g_quad_2d_std;
103  sln->set_quad_2d(quad);
104 
105  double norm = 0.0;
106  Element* e;
107  const Mesh* mesh = sln->get_mesh();
108 
109  for_all_active_elements(e, mesh)
110  {
111  // set maximum integration order for use in integrals, see limit_order()
112  update_limit_table(e->get_mode());
113 
114  sln->set_active_element(e);
115  RefMap* ru = sln->get_refmap();
116 
117  switch (norm_type)
118  {
119  case HERMES_L2_NORM:
120  norm += norm_fn_l2(sln, ru);
121  break;
122  case HERMES_H1_NORM:
123  norm += norm_fn_h1(sln, ru);
124  break;
125  case HERMES_HCURL_NORM:
126  norm += norm_fn_hc(sln, ru);
127  break;
128  case HERMES_HDIV_NORM:
129  norm += norm_fn_hdiv(sln, ru);
130  break;
131  default: throw Hermes::Exceptions::Exception("Unknown norm in calc_norm().");
132  }
133  }
134  return sqrt(norm);
135  }
136 
137  template<typename Scalar>
138  double Global<Scalar>::calc_norms(Hermes::vector<Solution<Scalar>*> slns)
139  {
140  // Calculate norms for all solutions.
141  Hermes::vector<double> norms;
142  int n = slns.size();
143  for (int i = 0; i<n; i++)
144  {
145  switch (slns[i]->get_space_type())
146  {
147  case HERMES_H1_SPACE: norms.push_back(calc_norm(slns[i], HERMES_H1_NORM)); break;
148  case HERMES_HCURL_SPACE: norms.push_back(calc_norm(slns[i], HERMES_HCURL_NORM)); break;
149  case HERMES_HDIV_SPACE: norms.push_back(calc_norm(slns[i], HERMES_HDIV_NORM)); break;
150  case HERMES_L2_SPACE: norms.push_back(calc_norm(slns[i], HERMES_L2_NORM)); break;
151  default: throw Hermes::Exceptions::Exception("Internal in calc_norms(): unknown space type.");
152  }
153  }
154  // Calculate the resulting norm.
155  double result = 0;
156  for (int i = 0; i < n; i++)
157  result += norms[i] * norms[i];
158  return sqrt(result);
159  }
160 
161  template<typename Scalar>
162  double Global<Scalar>::calc_abs_errors(Hermes::vector<Solution<Scalar>*> slns1, Hermes::vector<Solution<Scalar>*> slns2)
163  {
164  // Calculate errors for all solutions.
165  Hermes::vector<double> errors;
166  int n = slns1.size();
167  for (int i = 0; i < n; i++)
168  {
169  switch (slns1[i]->get_space_type())
170  {
171  case HERMES_H1_SPACE: errors.push_back(calc_abs_error(slns1[i], slns2[i], HERMES_H1_NORM)); break;
172  case HERMES_HCURL_SPACE: errors.push_back(calc_abs_error(slns1[i], slns2[i], HERMES_HCURL_NORM)); break;
173  case HERMES_HDIV_SPACE: errors.push_back(calc_abs_error(slns1[i], slns2[i], HERMES_HDIV_NORM)); break;
174  case HERMES_L2_SPACE: errors.push_back(calc_abs_error(slns1[i], slns2[i], HERMES_L2_NORM)); break;
175  default: throw Hermes::Exceptions::Exception("Internal in calc_norms(): unknown space type.");
176  }
177  }
178  // Calculate the resulting error.
179  double result = 0;
180  for (int i = 0; i < n; i++)
181  result += errors[i] * errors[i];
182  return sqrt(result);
183  }
184 
185  template<typename Scalar>
186  double Global<Scalar>::calc_rel_errors(Hermes::vector<Solution<Scalar>*> slns1, Hermes::vector<Solution<Scalar>*> slns2)
187  {
188  return calc_abs_errors(slns1, slns2) / calc_norms(slns2);
189  }
190 
191  template<typename Scalar>
192  double Global<Scalar>::error_fn_h1(MeshFunction<Scalar>* sln1, MeshFunction<Scalar>* sln2, RefMap* ru, RefMap* rv)
193  {
194  Quad2D* quad = sln1->get_quad_2d();
195 
196  int o = 2*std::max(sln1->get_fn_order(), sln2->get_fn_order()) + ru->get_inv_ref_order();
197  limit_order_nowarn(o, ru->get_active_element()->get_mode());
198 
199  sln1->set_quad_order(o);
200  sln2->set_quad_order(o);
201 
202  Scalar *uval, *vval, *dudx, *dudy, *dvdx, *dvdy;
203  uval = sln1->get_fn_values();
204  vval = sln2->get_fn_values();
205  sln1->get_dx_dy_values(dudx, dudy);
206  sln2->get_dx_dy_values(dvdx, dvdy);
207 
208  double result = 0.0;
209  h1_integrate_expression(Hermes::sqr(uval[i] - vval[i]) +
210  sqr(dudx[i] - dvdx[i]) + sqr(dudy[i] - dvdy[i]));
211  return result;
212  }
213 
214  template<typename Scalar>
215  double Global<Scalar>::norm_fn_h1(MeshFunction<Scalar>* sln, RefMap* ru)
216  {
217  Quad2D* quad = sln->get_quad_2d();
218 
219  int o = 2 * sln->get_fn_order() + ru->get_inv_ref_order();
220  limit_order_nowarn(o, ru->get_active_element()->get_mode());
221 
222  sln->set_quad_order(o);
223 
224  Scalar *uval, *dudx, *dudy;
225  uval = sln->get_fn_values();
226  sln->get_dx_dy_values(dudx, dudy);
227 
228  double result = 0.0;
229  h1_integrate_expression(sqr(uval[i]) + sqr(dudx[i]) + sqr(dudy[i]));
230  return result;
231  }
232 
233  template<typename Scalar>
234  double Global<Scalar>::error_fn_l2(MeshFunction<Scalar>* sln1, MeshFunction<Scalar>* sln2, RefMap* ru, RefMap* rv)
235  {
236  Quad2D* quad = sln1->get_quad_2d();
237 
238  int o = 2*std::max(sln1->get_fn_order(), sln2->get_fn_order()) + ru->get_inv_ref_order();
239  limit_order_nowarn(o, ru->get_active_element()->get_mode());
240 
241  sln1->set_quad_order(o, H2D_FN_VAL);
242  sln2->set_quad_order(o, H2D_FN_VAL);
243 
244  Scalar *uval, *vval;
245  uval = sln1->get_fn_values();
246  vval = sln2->get_fn_values();
247 
248  double result = 0.0;
249  h1_integrate_expression(Hermes::sqr(uval[i] - vval[i]));
250  return result;
251  }
252 
253  template<typename Scalar>
254  double Global<Scalar>::norm_fn_l2(MeshFunction<Scalar>* sln, RefMap* ru)
255  {
256  Quad2D* quad = sln->get_quad_2d();
257 
258  int o = 2 *sln->get_fn_order() + ru->get_inv_ref_order();
259  limit_order_nowarn(o, ru->get_active_element()->get_mode());
260 
261  sln->set_quad_order(o, H2D_FN_VAL);
262 
263  Scalar* uval = sln->get_fn_values();
264 
265  double result = 0.0;
266  h1_integrate_expression(Hermes::sqr(uval[i]));
267  return result;
268  }
269 
270  template<typename Scalar>
271  double Global<Scalar>::error_fn_hc(MeshFunction<Scalar>* sln1, MeshFunction<Scalar>* sln2, RefMap* ru, RefMap* rv)
272  {
273  Quad2D* quad = sln1->get_quad_2d();
274 
275  int o = 2 * std::max(sln1->get_fn_order(), sln2->get_fn_order()) + 2 + ru->get_inv_ref_order();
276  limit_order_nowarn(o, ru->get_active_element()->get_mode());
277 
278  sln1->set_quad_order(o);
279  sln2->set_quad_order(o);
280 
281  Scalar *uval0 = sln1->get_fn_values(0), *uval1 = sln1->get_fn_values(1);
282  Scalar *udx1 = sln1->get_dx_values(1), *udy0 = sln1->get_dy_values(0);
283  Scalar *vval0 = sln2->get_fn_values(0), *vval1 = sln2->get_fn_values(1);
284  Scalar *vdx1 = sln2->get_dx_values(1), *vdy0 = sln2->get_dy_values(0);
285 
286  double result = 0.0;
287  h1_integrate_expression(Hermes::sqr(uval0[i] - vval0[i]) + Hermes::sqr(uval1[i] - vval1[i]) +
288  sqr((udx1[i] - udy0[i]) - (vdx1[i] - vdy0[i])));
289  return result;
290  }
291 
292  template<typename Scalar>
293  double Global<Scalar>::norm_fn_hc(MeshFunction<Scalar>* sln, RefMap* ru)
294  {
295  Quad2D* quad = sln->get_quad_2d();
296 
297  int o = 2 * sln->get_fn_order() + 2 + ru->get_inv_ref_order();
298  limit_order_nowarn(o, ru->get_active_element()->get_mode());
299 
300  sln->set_quad_order(o);
301 
302  Scalar *uval0 = sln->get_fn_values(0), *uval1 = sln->get_fn_values(1);
303  Scalar *udx1 = sln->get_dx_values(1), *udy0 = sln->get_dy_values(0);
304 
305  double result = 0.0;
306  h1_integrate_expression(Hermes::sqr(uval0[i]) + Hermes::sqr(uval1[i]) + Hermes::sqr(udx1[i] - udy0[i]));
307  return result;
308  }
309 
310  template<typename Scalar>
311  double Global<Scalar>::error_fn_hcl2(MeshFunction<Scalar>* sln1, MeshFunction<Scalar>* sln2, RefMap* ru, RefMap* rv)
312  {
313  Quad2D* quad = sln1->get_quad_2d();
314 
315  int o = 2 * std::max(sln1->get_fn_order(), sln2->get_fn_order()) + 2 + ru->get_inv_ref_order();
316  limit_order_nowarn(o, ru->get_active_element()->get_mode());
317 
318  sln1->set_quad_order(o);
319  sln2->set_quad_order(o);
320 
321  Scalar *uval0 = sln1->get_fn_values(0), *uval1 = sln1->get_fn_values(1);
322  Scalar *vval0 = sln2->get_fn_values(0), *vval1 = sln2->get_fn_values(1);
323 
324  double result = 0.0;
325  h1_integrate_expression(Hermes::sqr(uval0[i] - vval0[i]) + Hermes::sqr(uval1[i] - vval1[i]));
326  return result;
327  }
328 
329  template<typename Scalar>
330  double Global<Scalar>::norm_fn_hcl2(MeshFunction<Scalar>* sln, RefMap* ru)
331  {
332  Quad2D* quad = sln->get_quad_2d();
333 
334  int o = 2 * sln->get_fn_order() + 2 + ru->get_inv_ref_order();
335  limit_order_nowarn(o, ru->get_active_element()->get_mode());
336 
337  sln->set_quad_order(o);
338 
339  Scalar *uval0 = sln->get_fn_values(0), *uval1 = sln->get_fn_values(1);
340  Scalar *udx1 = sln->get_dx_values(1), *udy0 = sln->get_dy_values(0);
341 
342  double result = 0.0;
343  h1_integrate_expression(Hermes::sqr(uval0[i]) + Hermes::sqr(uval1[i]));
344  return result;
345  }
346 
347  template<typename Scalar>
348  double Global<Scalar>::error_fn_hdiv(MeshFunction<Scalar>* sln1, MeshFunction<Scalar>* sln2, RefMap* ru, RefMap* rv)
349  {
350  throw Hermes::Exceptions::Exception("error_fn_hdiv() not implemented yet.");
351 
352  // Hcurl code
353  Quad2D* quad = sln1->get_quad_2d();
354 
355  int o = 2 * std::max(sln1->get_fn_order(), sln2->get_fn_order()) + 2 + ru->get_inv_ref_order();
356  limit_order_nowarn(o, ru->get_active_element()->get_mode());
357 
358  sln1->set_quad_order(o);
359  sln2->set_quad_order(o);
360 
361  Scalar *uval0 = sln1->get_fn_values(0), *uval1 = sln1->get_fn_values(1);
362  Scalar *udx1 = sln1->get_dx_values(1), *udy0 = sln1->get_dy_values(0);
363  Scalar *vval0 = sln2->get_fn_values(0), *vval1 = sln2->get_fn_values(1);
364  Scalar *vdx1 = sln2->get_dx_values(1), *vdy0 = sln2->get_dy_values(0);
365 
366  double result = 0.0;
367  h1_integrate_expression(Hermes::sqr(uval0[i] - vval0[i]) + Hermes::sqr(uval1[i] - vval1[i]) +
368  Hermes::sqr((udx1[i] - udy0[i]) - (vdx1[i] - vdy0[i])));
369  return result;
370  }
371 
372  template<typename Scalar>
373  double Global<Scalar>::norm_fn_hdiv(MeshFunction<Scalar>* sln, RefMap* ru)
374  {
375  throw Hermes::Exceptions::Exception("norm_fn_hdiv() not implemented yet.");
376 
377  // Hcurl code
378  Quad2D* quad = sln->get_quad_2d();
379 
380  int o = 2 * sln->get_fn_order() + 2 + ru->get_inv_ref_order();
381  limit_order_nowarn(o, ru->get_active_element()->get_mode());
382 
383  sln->set_quad_order(o);
384 
385  Scalar *uval0 = sln->get_fn_values(0), *uval1 = sln->get_fn_values(1);
386  Scalar *udx1 = sln->get_dx_values(1), *udy0 = sln->get_dy_values(0);
387 
388  double result = 0.0;
389  h1_integrate_expression(Hermes::sqr(uval0[i]) + Hermes::sqr(uval1[i]) + Hermes::sqr(udx1[i] - udy0[i]));
390  return result;
391  }
392 
393  template class HERMES_API Global<double>;
394  template class HERMES_API Global<std::complex<double> >;
395  }
396 }