Hermes2D  3.0
integrals_h1.h
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 #ifndef __H2D_INTEGRALS_H1_H
17 #define __H2D_INTEGRALS_H1_H
18 
19 #include "../quadrature/limit_order.h"
20 #include "../forms.h"
21 #include "../function/function.h"
22 
23 namespace Hermes
24 {
25  namespace Hermes2D
26  {
27  template<typename Real>
28  Real int_v(int n, double *wt, Func<Real> *v)
29  {
30  Real result = Real(0);
31  for (int i = 0; i < n; i++)
32  result += wt[i] * v->val[i];
33  return result;
34  }
35 
36  template<typename Real, typename Geom>
37  Real int_x_v(int n, double *wt, Func<Real> *v, Geom *e)
38  {
39  Real result = Real(0);
40  for (int i = 0; i < n; i++)
41  result += wt[i] * e->x[i] * v->val[i];
42  return result;
43  }
44 
45  template<typename Real, typename Geom>
46  Real int_y_v(int n, double *wt, Func<Real> *v, Geom *e)
47  {
48  Real result = Real(0);
49  for (int i = 0; i < n; i++)
50  result += wt[i] * e->y[i] * v->val[i];
51  return result;
52  }
53 
54  template<typename Real, typename Scalar>
55  Scalar int_u_v(int n, double *wt, Func<Real> *u, Func<Real> *v)
56  {
57  Scalar result = Scalar(0);
58  for (int i = 0; i < n; i++)
59  result += wt[i] * u->val[i] * v->val[i];
60  return result;
61  }
62 
63  // For residual forms.
64  template<typename Real, typename Scalar>
65  Scalar int_u_ext_v(int n, double *wt, Func<Scalar> *u_ext, Func<Real> *v)
66  {
67  Scalar result = Scalar(0);
68  for (int i = 0; i < n; i++)
69  result += wt[i] * u_ext->val[i] * v->val[i];
70  return result;
71  }
72 
73  template<typename Real, typename Scalar, typename Geom>
74  Scalar int_x_u_v(int n, double *wt, Func<Real> *u, Func<Real> *v, Geom *e)
75  {
76  Scalar result = Scalar(0);
77  for (int i = 0; i < n; i++)
78  result += wt[i] * e->x[i] * u->val[i] * v->val[i];
79  return result;
80  }
81 
82  // For residual forms.
83  template<typename Real, typename Scalar, typename Geom>
84  Scalar int_x_u_ext_v(int n, double *wt, Func<Scalar> *u_ext, Func<Real> *v, Geom *e)
85  {
86  Scalar result = Scalar(0);
87  for (int i = 0; i < n; i++)
88  result += wt[i] * e->x[i] * u_ext->val[i] * v->val[i];
89  return result;
90  }
91 
92  template<typename Real, typename Scalar, typename Geom>
93  Scalar int_y_u_v(int n, double *wt, Func<Real> *u, Func<Real> *v, Geom *e)
94  {
95  Scalar result = Scalar(0);
96  for (int i = 0; i < n; i++)
97  result += wt[i] * e->y[i] * u->val[i] * v->val[i];
98  return result;
99  }
100 
101  // For residual forms.
102  template<typename Real, typename Scalar, typename Geom>
103  Scalar int_y_u_ext_v(int n, double *wt, Func<Scalar> *u_ext, Func<Real> *v, Geom *e)
104  {
105  Scalar result = Scalar(0);
106  for (int i = 0; i < n; i++)
107  result += wt[i] * e->y[i] * u_ext->val[i] * v->val[i];
108  return result;
109  }
110 
111  template<typename Real, typename Scalar, typename Geom>
112  Scalar int_F_v(int n, double *wt, Real(*F)(Real x, Real y), Func<Real> *v, Geom *e)
113  {
114  Scalar result = Scalar(0);
115  for (int i = 0; i < n; i++)
116  result += wt[i] * ((*F)(e->x[i], e->y[i]) * v->val[i]);
117  return result;
118  }
119 
120  template<typename Real, typename Scalar>
121  Scalar int_grad_u_grad_v(int n, double *wt, Func<Real> *u, Func<Real> *v)
122  {
123  Scalar result = Scalar(0);
124  for (int i = 0; i < n; i++)
125  result += wt[i] * (u->dx[i] * v->dx[i] + u->dy[i] * v->dy[i]);
126  return result;
127  }
128 
129  // For residual forms.
130  template<typename Real, typename Scalar>
131  Scalar int_grad_u_ext_grad_v(int n, double *wt, Func<Scalar> *u_ext, Func<Real> *v)
132  {
133  Scalar result = Scalar(0);
134  for (int i = 0; i < n; i++)
135  result += wt[i] * (u_ext->dx[i] * v->dx[i] + u_ext->dy[i] * v->dy[i]);
136  return result;
137  }
138 
139  template<typename Real, typename Scalar, typename Geom>
140  Scalar int_x_grad_u_grad_v(int n, double *wt, Func<Real> *u, Func<Real> *v, Geom *e)
141  {
142  Scalar result = Scalar(0);
143  for (int i = 0; i < n; i++)
144  result += wt[i] * e->x[i] * (u->dx[i] * v->dx[i] + u->dy[i] * v->dy[i]);
145  return result;
146  }
147 
148  // For residual forms.
149  template<typename Real, typename Scalar, typename Geom>
150  Scalar int_x_grad_u_ext_grad_v(int n, double *wt, Func<Scalar> *u_ext, Func<Real> *v, Geom *e)
151  {
152  Scalar result = Scalar(0);
153  for (int i = 0; i < n; i++)
154  result += wt[i] * e->x[i] * (u_ext->dx[i] * v->dx[i] + u_ext->dy[i] * v->dy[i]);
155  return result;
156  }
157 
158  template<typename Real, typename Scalar, typename Geom>
159  Scalar int_y_grad_u_grad_v(int n, double *wt, Func<Real> *u, Func<Real> *v, Geom *e)
160  {
161  Scalar result = Scalar(0);
162  for (int i = 0; i < n; i++)
163  result += wt[i] * e->y[i] * (u->dx[i] * v->dx[i] + u->dy[i] * v->dy[i]);
164  return result;
165  }
166 
167  // For residual forms.
168  template<typename Real, typename Scalar, typename Geom>
169  Scalar int_y_grad_u_ext_grad_v(int n, double *wt, Func<Scalar> *u_ext, Func<Real> *v, Geom *e)
170  {
171  Scalar result = Scalar(0);
172  for (int i = 0; i < n; i++)
173  result += wt[i] * e->y[i] * (u_ext->dx[i] * v->dx[i] + u_ext->dy[i] * v->dy[i]);
174  return result;
175  }
176 
177  template<typename Real, typename Scalar>
178  Scalar int_dudx_v(int n, double *wt, Func<Real> *u, Func<Real> *v)
179  {
180  Scalar result = Scalar(0);
181  for (int i = 0; i < n; i++)
182  result += wt[i] * (u->dx[i] * v->val[i]);
183  return result;
184  }
185 
186  template<typename Real, typename Scalar>
187  Scalar int_dudy_v(int n, double *wt, Func<Real> *u, Func<Real> *v)
188  {
189  Scalar result = Scalar(0);
190  for (int i = 0; i < n; i++)
191  result += wt[i] * (u->dy[i] * v->val[i]);
192  return result;
193  }
194 
195  template<typename Real, typename Scalar>
196  Scalar int_u_dvdx(int n, double *wt, Func<Real> *u, Func<Real> *v)
197  {
198  Scalar result = Scalar(0);
199  for (int i = 0; i < n; i++)
200  result += wt[i] * u->val[i] * v->dx[i];
201  return result;
202  }
203 
204  template<typename Real, typename Scalar, typename Geom>
205  Scalar int_u_dvdx_over_x(int n, double *wt, Func<Real> *u, Func<Real> *v, Geom *e)
206  {
207  Scalar result = Scalar(0);
208  for (int i = 0; i < n; i++)
209  result += wt[i] * ((e->x[i] > 0) ? u->val[i] * v->dx[i] / e->x[i] : 0.0);
210  return result;
211  }
212 
213  // For residual forms.
214  template<typename Real, typename Scalar, typename Geom>
215  Scalar int_u_ext_dvdx_over_x(int n, double *wt, Func<Scalar> *u_ext, Func<Real> *v, Geom *e)
216  {
217  Scalar result = Scalar(0);
218  for (int i = 0; i < n; i++)
219  result += wt[i] * u_ext->val[i] * v->dx[i] / e->x[i];
220  return result;
221  }
222 
223  template<typename Real, typename Scalar>
224  Scalar int_u_dvdy(int n, double *wt, Func<Real> *u, Func<Real> *v)
225  {
226  Scalar result = Scalar(0);
227  for (int i = 0; i < n; i++)
228  result += wt[i] * u->val[i] * v->dy[i];
229  return result;
230  }
231 
232  template<typename Real, typename Scalar, typename Geom>
233  Scalar int_u_dvdy_over_y(int n, double *wt, Func<Real> *u, Func<Real> *v, Geom *e)
234  {
235  Scalar result = Scalar(0);
236  for (int i = 0; i < n; i++)
237  result += wt[i] * ((e->y[i] > 0) ? u->val[i] * v->dy[i] / e->y[i] : 0.0);
238  return result;
239  }
240 
241  // For residual forms.
242  template<typename Real, typename Scalar, typename Geom>
243  Scalar int_u_ext_dvdy_over_y(int n, double *wt, Func<Scalar> *u_ext, Func<Real> *v, Geom *e)
244  {
245  Scalar result = Scalar(0);
246  for (int i = 0; i < n; i++)
247  result += wt[i] * u_ext->val[i] * v->dy[i] / e->y[i];
248  return result;
249  }
250 
251  template<typename Real, typename Scalar>
252  Scalar int_dudx_dvdx(int n, double *wt, Func<Real> *u, Func<Real> *v)
253  {
254  Scalar result = Scalar(0);
255  for (int i = 0; i < n; i++)
256  result += wt[i] * (u->dx[i] * v->dx[i]);
257  return result;
258  }
259 
260  template<typename Real, typename Scalar>
261  Scalar int_dudy_dvdy(int n, double *wt, Func<Real> *u, Func<Real> *v)
262  {
263  Scalar result = Scalar(0);
264  for (int i = 0; i < n; i++)
265  result += wt[i] * (u->dy[i] * v->dy[i]);
266  return result;
267  }
268 
269  template<typename Real, typename Scalar>
270  Scalar int_dudx_dvdy(int n, double *wt, Func<Real> *u, Func<Real> *v)
271  {
272  Scalar result = Scalar(0);
273  for (int i = 0; i < n; i++)
274  result += wt[i] * (u->dx[i] * v->dy[i]);
275  return result;
276  }
277 
278  template<typename Real, typename Scalar>
279  Scalar int_dudy_dvdx(int n, double *wt, Func<Real> *u, Func<Real> *v)
280  {
281  Scalar result = Scalar(0);
282  for (int i = 0; i < n; i++)
283  result += wt[i] * (v->dx[i] * u->dy[i]);
284  return result;
285  }
286 
287  template<typename Real, typename Scalar>
288  Scalar int_w_nabla_u_v(int n, double *wt, Func<Real> *w1, Func<Real> *w2,
289  Func<Real> *u, Func<Real> *v)
290  {
291  Scalar result = Scalar(0);
292  for (int i = 0; i < n; i++)
293  result += wt[i] * (w1->val[i] * u->dx[i] + w2->val[i] * u->dy[i]) * v->val[i];
294  return result;
295  }
296 
298 
299  // the inner integration loops for both constant and non-constant jacobian elements
300  // for expression without partial derivatives - the variables e, quad, o must be already
301  // defined and initialized
302 #define h1_integrate_expression(exp) \
303  {double3* pt = quad->get_points(o, ru->get_active_element()->get_mode()); \
304  unsigned char np = quad->get_num_points(o, ru->get_active_element()->get_mode()); \
305  if (ru->is_jacobian_const()){ \
306  for (int i = 0; i < np; i++) \
307  result += pt[i][2] * (exp); \
308  result *= ru->get_const_jacobian(); \
309  } \
310  else { \
311  double* jac = ru->get_jacobian(o); \
312  for (int i = 0; i < np; i++) \
313  result += pt[i][2] * jac[i] * (exp); \
314  }}
315 
317 
318  template<typename Scalar>
319  inline double int_h1_error(Function<Scalar>* fu, Function<Scalar>* fv, RefMap* ru, RefMap* rv)
320  {
321  Quad2D* quad = fu->get_quad_2d();
322  assert(quad == fv->get_quad_2d());
323 
324  int o = std::max(2 * fu->get_fn_order(), 2 * fv->get_fn_order()) + ru->get_inv_ref_order();
325  limit_order(o, ru->get_active_element()->get_mode());
326  fu->set_quad_order(o);
327  fv->set_quad_order(o);
328 
329  Scalar* fnu = fu->get_fn_values();
330  Scalar* fnv = fv->get_fn_values();
331 
332  Scalar *dudx, *dudy, *dvdx, *dvdy;
333  fu->get_dx_dy_values(dudx, dudy);
334  fv->get_dx_dy_values(dvdx, dvdy);
335 
336  double result = 0.0;
337  h1_integrate_expression(sqr(fnu[i] - fnv[i]) + sqr(dudx[i] - dvdx[i]) + sqr(dudy[i] - dvdy[i]));
338  return result;
339  }
340 
341  template<typename Scalar>
342  inline double int_h1_semi_error(Function<Scalar>* fu, Function<Scalar>* fv, RefMap* ru, RefMap* rv)
343  {
344  Quad2D* quad = fu->get_quad_2d();
345  assert(quad == fv->get_quad_2d());
346 
347  int o = std::max(2 * fu->get_fn_order(), 2 * fv->get_fn_order()) + ru->get_inv_ref_order();
348  limit_order(o, ru->get_active_element()->get_mode());
349  fu->set_quad_order(o);
350  fv->set_quad_order(o);
351 
352  Scalar* fnu = fu->get_fn_values();
353  Scalar* fnv = fv->get_fn_values();
354 
355  Scalar *dudx, *dudy, *dvdx, *dvdy;
356  fu->get_dx_dy_values(dudx, dudy);
357  fv->get_dx_dy_values(dvdx, dvdy);
358 
359  double result = 0.0;
360  h1_integrate_expression(sqr(dudx[i] - dvdx[i]) + sqr(dudy[i] - dvdy[i]));
361  return result;
362  }
363 
364  template<typename Scalar>
365  inline double int_l2_error(Function<Scalar>* fu, Function<Scalar>* fv, RefMap* ru, RefMap* rv)
366  {
367  Quad2D* quad = fu->get_quad_2d();
368  assert(quad == fv->get_quad_2d());
369 
370  int o = std::max(2 * fu->get_fn_order(), 2 * fv->get_fn_order()) + ru->get_inv_ref_order();
371  limit_order(o, ru->get_active_element()->get_mode());
372  fu->set_quad_order(o, H2D_FN_VAL);
373  fv->set_quad_order(o, H2D_FN_VAL);
374 
375  Scalar* fnu = fu->get_fn_values();
376  Scalar* fnv = fv->get_fn_values();
377 
378  double result = 0.0;
379  h1_integrate_expression(sqr(fnu[i] - fnv[i]));
380  return result;
381  }
382 
383  template<typename Scalar>
384  inline double int_dx_error(Function<Scalar>* fu, Function<Scalar>* fv, RefMap* ru, RefMap* rv)
385  {
386  Quad2D* quad = fu->get_quad_2d();
387  assert(quad == fv->get_quad_2d());
388 
389  int o = std::max(2 * fu->get_fn_order(), 2 * fv->get_fn_order()) + ru->get_inv_ref_order();
390  limit_order(o, ru->get_active_element()->get_mode());
391  fu->set_quad_order(o);
392  fv->set_quad_order(o);
393 
394  Scalar *dudx, *dudy, *dvdx, *dvdy;
395  fu->get_dx_dy_values(dudx, dudy);
396  fv->get_dx_dy_values(dvdx, dvdy);
397 
398  double result = 0.0;
399  h1_integrate_expression(sqr(dudx[i] - dvdx[i]));
400  return result;
401  }
402 
403  template<typename Scalar>
404  inline double int_dy_error(Function<Scalar>* fu, Function<Scalar>* fv, RefMap* ru, RefMap* rv)
405  {
406  Quad2D* quad = fu->get_quad_2d();
407  assert(quad == fv->get_quad_2d());
408 
409  int o = std::max(2 * fu->get_fn_order(), 2 * fv->get_fn_order()) + ru->get_inv_ref_order();
410  limit_order(o, ru->get_active_element()->get_mode());
411  fu->set_quad_order(o);
412  fv->set_quad_order(o);
413 
414  Scalar *dudx, *dudy, *dvdx, *dvdy;
415  fu->get_dx_dy_values(dudx, dudy);
416  fv->get_dx_dy_values(dvdx, dvdy);
417 
418  double result = 0.0;
419  h1_integrate_expression(sqr(dudy[i] - dvdy[i]));
420  return result;
421  }
422 
423  template<typename Scalar>
424  inline double int_h1_norm(Function<Scalar>* fu, RefMap* ru)
425  {
426  Quad2D* quad = fu->get_quad_2d();
427 
428  int o = 2 * fu->get_fn_order() + ru->get_inv_ref_order();
429  limit_order(o, ru->get_active_element()->get_mode());
430  fu->set_quad_order(o);
431 
432  Scalar* fnu = fu->get_fn_values();
433  Scalar *dudx, *dudy;
434  fu->get_dx_dy_values(dudx, dudy);
435 
436  double result = 0.0;
437  h1_integrate_expression(sqr(fnu[i]) + sqr(dudx[i]) + sqr(dudy[i]));
438  return result;
439  }
440 
441  template<typename Scalar>
442  inline double int_h1_seminorm(Function<Scalar>* fu, RefMap* ru)
443  {
444  Quad2D* quad = fu->get_quad_2d();
445 
446  int o = 2 * fu->get_fn_order() + ru->get_inv_ref_order();
447  limit_order(o, ru->get_active_element()->get_mode());
448  fu->set_quad_order(o);
449 
450  Scalar* fnu = fu->get_fn_values();
451  Scalar *dudx, *dudy;
452  fu->get_dx_dy_values(dudx, dudy);
453 
454  double result = 0.0;
455  h1_integrate_expression(sqr(dudx[i]) + sqr(dudy[i]));
456  return result;
457  }
458 
459  template<typename Scalar>
460  inline double int_l2_norm(Function<Scalar>* fu, RefMap* ru)
461  {
462  Quad2D* quad = fu->get_quad_2d();
463 
464  int o = 2 * fu->get_fn_order() + ru->get_inv_ref_order();
465  limit_order(o, ru->get_active_element()->get_mode());
466  fu->set_quad_order(o, H2D_FN_VAL);
467 
468  Scalar* fnu = fu->get_fn_values();
469 
470  double result = 0.0;
471  h1_integrate_expression(sqr(fnu[i]));
472  return result;
473  }
474  }
475 }
476 #endif
Hermes
Definition: adapt.h:24
Hermes::Hermes2D::H2D_FN_VAL
const int H2D_FN_VAL
Both components are usually requested together...
Definition: function.h:53