Hermes2D  3.0
refmap.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 "mesh.h"
18 #include "refmap.h"
19 
20 namespace Hermes
21 {
22  namespace Hermes2D
23  {
24  RefMap::RefMap() : ref_map_shapeset(H1ShapesetJacobi()), ref_map_pss(&ref_map_shapeset)
25  {
26  quad_2d = nullptr;
27  set_quad_2d(&g_quad_2d_std);
28  this->reinit_storage();
29  }
30 
31  RefMap::~RefMap() { }
32 
33  Quad2D* RefMap::get_quad_2d() const
34  {
35  return quad_2d;
36  }
37 
38  double* RefMap::get_phys_x(int order)
39  {
40  if (order != this->phys_x_calculated)
41  this->calc_phys_x(order);
42  return this->phys_x;
43  }
44 
45  double* RefMap::get_phys_y(int order)
46  {
47  if (order != this->phys_y_calculated)
48  this->calc_phys_y(order);
49  return this->phys_y;
50  }
51 
52  double3* RefMap::get_tangent(int edge, int order)
53  {
54  if (order == -1)
55  order = quad_2d->get_edge_points(edge, quad_2d->get_max_order(element->get_mode()), element->get_mode());
56 
57  if (order != this->tan_calculated[edge])
58  this->calc_tangent(edge, order);
59 
60  return tan[edge];
61  }
62 
63  void RefMap::set_quad_2d(Quad2D* quad_2d)
64  {
65  this->quad_2d = quad_2d;
66  ref_map_pss.set_quad_2d(quad_2d);
67  this->reinit_storage();
68  }
69 
70  void RefMap::set_active_element(Element* e)
71  {
72  this->reinit_storage();
73 
74  Transformable::set_active_element(e);
75  ref_map_pss.set_active_element(e);
76 
77  this->is_const = element->has_const_ref_map();
78 
79  // prepare the shapes and coefficients of the reference map
80  unsigned short j, k = 0;
81  for (unsigned char i = 0; i < e->get_nvert(); i++)
82  indices[k++] = ref_map_shapeset.get_vertex_index(i, e->get_mode());
83 
84  // straight-edged element
85  if (e->cm == nullptr)
86  {
87  for (unsigned char i = 0; i < e->get_nvert(); i++)
88  {
89  lin_coeffs[i][0] = e->vn[i]->x;
90  lin_coeffs[i][1] = e->vn[i]->y;
91  }
92  coeffs = lin_coeffs;
93  nc = e->get_nvert();
94  }
95  else // curvilinear element - edge and bubble shapes
96  {
97  unsigned short o = e->cm->order;
98  for (unsigned char i = 0; i < e->get_nvert(); i++)
99  for (j = 2; j <= o; j++)
100  indices[k++] = ref_map_shapeset.get_edge_index(i, 0, j, e->get_mode());
101 
102  if (e->is_quad()) o = H2D_MAKE_QUAD_ORDER(o, o);
103  memcpy(indices + k, ref_map_shapeset.get_bubble_indices(o, e->get_mode()),
104  ref_map_shapeset.get_num_bubbles(o, e->get_mode()) * sizeof(int));
105 
106  coeffs = e->cm->coeffs;
107  nc = e->cm->nc;
108  }
109 
110  this->inv_ref_order = this->element->iro_cache;
111 
112  // constant inverse reference map
113  if (is_const)
114  calc_const_inv_ref_map();
115  else
116  const_jacobian = 0.0;
117  }
118 
119  void RefMap::set_element_iro_cache(Element* element)
120  {
121  bool is_const = element->has_const_ref_map();
122 
123  if (is_const)
124  {
125 #ifdef _DEBUG
126  assert(element->iro_cache == 0);
127 #endif
128  return;
129  }
130 #pragma omp critical (element_iro_cache_setting)
131  {
132  RefMap rm;
133  rm.set_active_element(element);
134  element->iro_cache = rm.calc_inv_ref_order();
135  }
136  }
137 
138  void RefMap::reinit_storage()
139  {
140  jacobian_calculated = -1;
141  inv_ref_map_calculated = -1;
142  second_ref_map_calculated = -1;
143 
144  direct_ref_map_calculated = -1;
145  phys_x_calculated = -1;
146  phys_y_calculated = -1;
147  tan_calculated[0] = tan_calculated[1] = tan_calculated[2] = tan_calculated[3] = -1;
148  }
149 
150  void RefMap::push_transform(int son)
151  {
152  Transformable::push_transform(son);
153  const_jacobian *= 0.25;
154  this->reinit_storage();
155  }
156 
157  void RefMap::pop_transform()
158  {
159  Transformable::pop_transform();
160  const_jacobian *= 4.;
161  this->reinit_storage();
162  }
163 
164  void RefMap::calc_inv_ref_map(int order)
165  {
166  int i, j, np = quad_2d->get_num_points(order, element->get_mode());
167 
168  // construct jacobi matrices of the direct reference map for all integration points
169  ref_map_pss.force_transform(sub_idx, ctm);
170 
171  double* m_00 = this->direct_ref_map[0][0];
172  double* m_01 = this->direct_ref_map[0][1];
173  double* m_10 = this->direct_ref_map[1][0];
174  double* m_11 = this->direct_ref_map[1][1];
175 
176  for (j = 0; j < np; j++)
177  m_00[j] = m_01[j] = m_10[j] = m_11[j] = 0.;
178 
179  for (i = 0; i < nc; i++)
180  {
181  double coeff_0 = coeffs[i][0];
182  double coeff_1 = coeffs[i][1];
183  ref_map_pss.set_active_shape(indices[i]);
184  ref_map_pss.set_quad_order(order, H2D_FN_DX_0 | H2D_FN_DY_0);
185  const double *dx = ref_map_pss.get_dx_values();
186  const double *dy = ref_map_pss.get_dy_values();
187  for (j = 0; j < np; j++)
188  {
189  double dx_ = dx[j];
190  double dy_ = dy[j];
191 
192  m_00[j] += coeff_0 * dx_;
193  m_01[j] += coeff_0 * dy_;
194  m_10[j] += coeff_1 * dx_;
195  m_11[j] += coeff_1 * dy_;
196  }
197  }
198 
199  // calculate the jacobian and inverted matrix
200  double trj = get_transform_jacobian();
201  double2x2* irm = this->inv_ref_map;
202  double* jac = this->jacobian;
203 
204  for (i = 0; i < np; i++)
205  {
206  jac[i] = (m_00[i] * m_11[i] - m_01[i] * m_10[i]);
207  double ij = 1.0 / jac[i];
208 
209  if (!finite(ij))
210  throw Hermes::Exceptions::Exception("1/jac[%d] is infinity when calculating inv. ref. map on element %d for order %d (jac = %g)", i, this->element->id, order, ij);
211  if (ij != ij)
212  throw Hermes::Exceptions::Exception("1/jac[%d] is NaN when calculating inv. ref. map on element %d for order %d (jac = %g)", i, this->element->id, order, ij);
213 
214  // invert and transpose the matrix
215  irm[i][0][0] = m_11[i] * ij;
216  irm[i][0][1] = -m_10[i] * ij;
217  irm[i][1][0] = -m_01[i] * ij;
218  irm[i][1][1] = m_00[i] * ij;
219 
220  jac[i] *= trj;
221  }
222 
223  this->inv_ref_map_calculated = order;
224  this->jacobian_calculated = order;
225  }
226 
227  void RefMap::calc_const_inv_ref_map()
228  {
229  int k = element->is_triangle() ? 2 : 3;
230  double m[2][2] = { { element->vn[1]->x - element->vn[0]->x, element->vn[k]->x - element->vn[0]->x },
231  { element->vn[1]->y - element->vn[0]->y, element->vn[k]->y - element->vn[0]->y } };
232 
233  const_jacobian = 0.25 * (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
234 
235  double ij = 0.5 / const_jacobian;
236 
237  const_inv_ref_map[0][0] = m[1][1] * ij;
238  const_inv_ref_map[1][0] = -m[0][1] * ij;
239  const_inv_ref_map[0][1] = -m[1][0] * ij;
240  const_inv_ref_map[1][1] = m[0][0] * ij;
241 
242  const_jacobian *= get_transform_jacobian();
243  }
244 
245  void RefMap::calc_phys_x(int order)
246  {
247  // transform all x coordinates of the integration points
248  int i, j, np = quad_2d->get_num_points(order, element->get_mode());
249  double* x = this->phys_x;
250  memset(x, 0, np * sizeof(double));
251  ref_map_pss.force_transform(sub_idx, ctm);
252  for (i = 0; i < nc; i++)
253  {
254  ref_map_pss.set_active_shape(indices[i]);
255  ref_map_pss.set_quad_order(order, H2D_FN_VAL_0);
256  const double* fn = ref_map_pss.get_fn_values();
257  for (j = 0; j < np; j++)
258  x[j] += coeffs[i][0] * fn[j];
259  }
260  this->phys_x_calculated = order;
261  }
262 
263  void RefMap::calc_phys_y(int order)
264  {
265  // transform all y coordinates of the integration points
266  int i, j, np = quad_2d->get_num_points(order, element->get_mode());
267  double* y = this->phys_y;
268  memset(y, 0, np * sizeof(double));
269  ref_map_pss.force_transform(sub_idx, ctm);
270  for (i = 0; i < nc; i++)
271  {
272  ref_map_pss.set_active_shape(indices[i]);
273  ref_map_pss.set_quad_order(order, H2D_FN_VAL_0);
274  const double* fn = ref_map_pss.get_fn_values();
275  for (j = 0; j < np; j++)
276  y[j] += coeffs[i][1] * fn[j];
277  }
278  this->phys_y_calculated = order;
279  }
280 
281  void RefMap::calc_tangent(int edge, int eo)
282  {
283  int i, j;
284  unsigned char np = quad_2d->get_num_points(eo, element->get_mode());
285  double3* tan = this->tan[edge];
286  int a = edge, b = element->next_vert(edge);
287 
288  if (!element->is_curved())
289  {
290  // straight edges: the tangent at each point is just the edge length
291  tan[0][0] = element->vn[b]->x - element->vn[a]->x;
292  tan[0][1] = element->vn[b]->y - element->vn[a]->y;
293  tan[0][2] = sqrt(sqr(tan[0][0]) + sqr(tan[0][1]));
294  double inorm = 1.0 / tan[0][2];
295  tan[0][0] *= inorm;
296  tan[0][1] *= inorm;
297  tan[0][2] *= (edge == 0 || edge == 2) ? ctm->m[0] : ctm->m[1];
298 
299  for (i = 1; i < np; i++)
300  memcpy(tan + i, tan, sizeof(double3));
301  }
302  else
303  {
304  // construct jacobi matrices of the direct reference map at integration points along the edge
305  double2x2 m[15];
306  assert(np <= 15);
307  memset(m, 0, np*sizeof(double2x2));
308  ref_map_pss.force_transform(sub_idx, ctm);
309  for (i = 0; i < nc; i++)
310  {
311  ref_map_pss.set_active_shape(indices[i]);
312  ref_map_pss.set_quad_order(eo);
313  const double *dx = ref_map_pss.get_dx_values();
314  const double *dy = ref_map_pss.get_dy_values();
315  for (j = 0; j < np; j++)
316  {
317  m[j][0][0] += coeffs[i][0] * dx[j];
318  m[j][0][1] += coeffs[i][0] * dy[j];
319  m[j][1][0] += coeffs[i][1] * dx[j];
320  m[j][1][1] += coeffs[i][1] * dy[j];
321  }
322  }
323 
324  // multiply them by the vector of the reference edge
325  double2* v1 = ref_map_shapeset.get_ref_vertex(a, element->get_mode());
326  double2* v2 = ref_map_shapeset.get_ref_vertex(b, element->get_mode());
327 
328  double ex = (*v2)[0] - (*v1)[0];
329  double ey = (*v2)[1] - (*v1)[1];
330  for (i = 0; i < np; i++)
331  {
332  double3& t = tan[i];
333  t[0] = m[i][0][0] * ex + m[i][0][1] * ey;
334  t[1] = m[i][1][0] * ex + m[i][1][1] * ey;
335  t[2] = sqrt(sqr(t[0]) + sqr(t[1]));
336  double inorm = 1.0 / t[2];
337  t[0] *= inorm;
338  t[1] *= inorm;
339  t[2] *= (edge == 0 || edge == 2) ? ctm->m[0] : ctm->m[1];
340  }
341  }
342 
343  this->tan_calculated[edge] = eo;
344  }
345 
346  int RefMap::calc_inv_ref_order()
347  {
348  Quad2D* quad = get_quad_2d();
349  int i, o, mo = (quad->get_max_order(element->get_mode()) / 2);
350 
351  // check first the positivity of the jacobian
352  double3* pt = quad->get_points(mo, element->get_mode());
353  double2x2* m = get_inv_ref_map(mo);
354  double* jac = get_jacobian(mo);
355  for (i = 0; i < quad->get_num_points(mo, element->get_mode()); i++)
356  if (jac[i] <= 0.0)
357  throw Hermes::Exceptions::Exception("Element #%d is concave or badly oriented.", element->id);
358 
359  // next, estimate the "exact" value of the typical integral int_grad_u_grad_v
360  // (with grad_u == grad_v == (1, 1)) using the maximum integration rule
361  double exact1 = 0.0;
362  double exact2 = 0.0;
363  for (i = 0; i < quad->get_num_points(mo, element->get_mode()); i++, m++)
364  {
365  exact1 += pt[i][2] * jac[i] * (sqr((*m)[0][0] + (*m)[0][1]) + sqr((*m)[1][0] + (*m)[1][1]));
366  exact2 += pt[i][2] / jac[i];
367  }
368  // find sufficient quadrature degree
369  for (o = 0; o < mo; o++)
370  {
371  pt = quad->get_points(o, element->get_mode());
372  m = get_inv_ref_map(o);
373  jac = get_jacobian(o);
374  double result1 = 0.0;
375  double result2 = 0.0;
376  for (i = 0; i < quad->get_num_points(o, element->get_mode()); i++, m++)
377  {
378  result1 += pt[i][2] * jac[i] * (sqr((*m)[0][0] + (*m)[0][1]) + sqr((*m)[1][0] + (*m)[1][1]));
379  result2 += pt[i][2] / jac[i];
380  }
381  if ((fabs((exact1 - result1) / exact1) < Hermes::HermesSqrtEpsilon) && (fabs((exact2 - result2) / exact2) < Hermes::HermesSqrtEpsilon))
382  break;
383  }
384  if (o >= 10)
385  {
386  this->warn("Element #%d is too distorted (iro ~ %d).", element->id, o);
387  }
388  return o;
389  }
390 
391  void RefMap::inv_ref_map_at_point(double xi1, double xi2, double& x, double& y, double2x2& m)
392  {
393  double2x2 tmp;
394  memset(tmp, 0, sizeof(double2x2));
395  x = y = 0;
396  for (int i = 0; i < nc; i++)
397  {
398  double val = ref_map_shapeset.get_fn_value(indices[i], xi1, xi2, 0, element->get_mode());
399  x += coeffs[i][0] * val;
400  y += coeffs[i][1] * val;
401 
402  double dx = ref_map_shapeset.get_dx_value(indices[i], xi1, xi2, 0, element->get_mode());
403  double dy = ref_map_shapeset.get_dy_value(indices[i], xi1, xi2, 0, element->get_mode());
404  tmp[0][0] += coeffs[i][0] * dx;
405  tmp[0][1] += coeffs[i][0] * dy;
406  tmp[1][0] += coeffs[i][1] * dx;
407  tmp[1][1] += coeffs[i][1] * dy;
408  }
409 
410  // inverse matrix
411  double jac = tmp[0][0] * tmp[1][1] - tmp[0][1] * tmp[1][0];
412  m[0][0] = tmp[1][1] / jac;
413  m[0][1] = -tmp[1][0] / jac;
414  m[1][0] = -tmp[0][1] / jac;
415  m[1][1] = tmp[0][0] / jac;
416  }
417 
418 #ifdef H2D_USE_SECOND_DERIVATIVES
419  void RefMap::calc_second_ref_map(int order)
420  {
421  int i, j, np = quad_2d->get_num_points(order, element->get_mode());
422 
423  double3x2* k = calloc_with_check<double3x2>(np);
424  ref_map_pss.force_transform(sub_idx, ctm);
425  for (i = 0; i < nc; i++)
426  {
427  ref_map_pss.set_active_shape(indices[i]);
428  ref_map_pss.set_quad_order(order, H2D_FN_ALL);
429  const double *dxx = ref_map_pss.get_dxx_values();
430  const double *dyy = ref_map_pss.get_dyy_values();
431  const double *dxy = ref_map_pss.get_dxy_values();
432  for (j = 0; j < np; j++)
433  {
434  k[j][0][0] += coeffs[i][0] * dxx[j];
435  k[j][0][1] += coeffs[i][1] * dxx[j];
436  k[j][1][0] += coeffs[i][0] * dxy[j];
437  k[j][1][1] += coeffs[i][1] * dxy[j];
438  k[j][2][0] += coeffs[i][0] * dyy[j];
439  k[j][2][1] += coeffs[i][1] * dyy[j];
440  }
441  }
442 
443  double3x2* mm = this->second_ref_map;
444  double2x2* m = get_inv_ref_map(order);
445  for (j = 0; j < np; j++)
446  {
447  double a, b;
448  // coefficients in second derivative with respect to xx
449  a = sqr(m[j][0][0])*k[j][0][0] + 2 * m[j][0][0] * m[j][0][1] * k[j][1][0] + sqr(m[j][0][1])*k[j][2][0];
450  b = sqr(m[j][0][0])*k[j][0][1] + 2 * m[j][0][0] * m[j][0][1] * k[j][1][1] + sqr(m[j][0][1])*k[j][2][1];
451  // du/dx
452  mm[j][0][0] = -(a * m[j][0][0] + b * m[j][1][0]);
453  // du/dy
454  mm[j][0][1] = -(a * m[j][0][1] + b * m[j][1][1]);
455 
456  // coefficients in second derivative with respect to xy
457  a = m[j][0][0] * m[j][1][0] * k[j][0][0] + (m[j][0][1] * m[j][1][0] + m[j][0][0] * m[j][1][1])*k[j][1][0] + m[j][0][1] * m[j][1][1] * k[j][2][0];
458  b = m[j][0][0] * m[j][1][0] * k[j][0][1] + (m[j][0][1] * m[j][1][0] + m[j][0][0] * m[j][1][1])*k[j][1][1] + m[j][0][1] * m[j][1][1] * k[j][2][1];
459  // du/dx
460  mm[j][1][0] = -(a * m[j][0][0] + b * m[j][1][0]);
461  // du/dy
462  mm[j][1][1] = -(a * m[j][0][1] + b * m[j][1][1]);
463 
464  // coefficients in second derivative with respect to yy
465  a = sqr(m[j][1][0])*k[j][0][0] + 2 * m[j][1][0] * m[j][1][1] * k[j][1][0] + sqr(m[j][1][1])*k[j][2][0];
466  b = sqr(m[j][1][0])*k[j][0][1] + 2 * m[j][1][0] * m[j][1][1] * k[j][1][1] + sqr(m[j][1][1])*k[j][2][1];
467  // du/dx
468  mm[j][2][0] = -(a * m[j][0][0] + b * m[j][1][0]);
469  // du/dy
470  mm[j][2][1] = -(a * m[j][0][1] + b * m[j][1][1]);
471  }
472 
473  free_with_check(k);
474 
475  this->second_ref_map_calculated = order;
476  }
477 
478  void RefMap::second_ref_map_at_point(double xi1, double xi2, double& x, double& y, double3x2& mm)this-
479  {
480  double3x2 k;
481  memset(k, 0, sizeof(double3x2));
482  x = y = 0;
483  for (int i = 0; i < nc; i++)
484  {
485  double val = ref_map_shapeset.get_fn_value(indices[i], xi1, xi2, 0, element->get_mode());
486  x += coeffs[i][0] * val;
487  y += coeffs[i][1] * val;
488 
489  double dxy, dxx, dyy;
490 
491  dxx = ref_map_shapeset.get_dxx_value(indices[i], xi1, xi2, 0, element->get_mode());
492  dxy = ref_map_shapeset.get_dxy_value(indices[i], xi1, xi2, 0, element->get_mode());
493  dyy = ref_map_shapeset.get_dxy_value(indices[i], xi1, xi2, 0, element->get_mode());
494 
495  k[0][0] += coeffs[i][0] * dxx;
496  k[0][1] += coeffs[i][1] * dxx;
497  k[1][0] += coeffs[i][0] * dxy;
498  k[1][1] += coeffs[i][1] * dxy;
499  k[2][0] += coeffs[i][0] * dyy;
500  k[2][1] += coeffs[i][1] * dyy;
501  }
502 
503  double2x2 m;
504  this->inv_ref_map_at_point(xi1, xi2, x, y, m);
505  double a, b;
506 
507  // coefficients in second derivative with respect to xx
508  a = sqr(m[0][0])*k[0][0] + 2 * m[0][0] * m[0][1] * k[1][0] + sqr(m[0][1])*k[2][0];
509  b = sqr(m[0][0])*k[0][1] + 2 * m[0][0] * m[0][1] * k[1][1] + sqr(m[0][1])*k[2][1];
510  // du/dx
511  mm[0][0] = -(a * m[0][0] + b * m[1][0]);
512  // du/dy
513  mm[0][1] = -(a * m[0][1] + b * m[1][1]);
514 
515  // coefficients in second derivative with respect to xy
516  a = m[0][0] * m[1][0] * k[0][0] + (m[0][1] * m[1][0] + m[0][0] * m[1][1])*k[1][0] + m[0][1] * m[1][1] * k[2][0];
517  b = m[0][0] * m[1][0] * k[0][1] + (m[0][1] * m[1][0] + m[0][0] * m[1][1])*k[1][1] + m[0][1] * m[1][1] * k[2][1];
518  // du/dx
519  mm[1][0] = -(a * m[0][0] + b * m[1][0]);
520  // du/dy
521  mm[1][1] = -(a * m[0][1] + b * m[1][1]);
522 
523  // coefficients in second derivative with respect to yy
524  a = sqr(m[1][0])*k[0][0] + 2 * m[1][0] * m[1][1] * k[1][0] + sqr(m[1][1])*k[2][0];
525  b = sqr(m[1][0])*k[0][1] + 2 * m[1][0] * m[1][1] * k[1][1] + sqr(m[1][1])*k[2][1];
526  // du/dx
527  mm[2][0] = -(a * m[0][0] + b * m[1][0]);
528  // du/dy
529  mm[2][1] = -(a * m[0][1] + b * m[1][1]);
530  }
531 
532  double3x2* RefMap::get_second_ref_map(int order)
533  {
534  if (this->is_const)
535  throw Hermes::Exceptions::Exception("RefMap::get_second_ref_map() called with a const jacobian.");
536  if (order != this->second_ref_map_calculated)
537  this->calc_second_ref_map(order);
538  return this->second_ref_map;
539  }
540 #endif
541 
542  void RefMap::untransform(Element* e, double x, double y, double& xi1, double& xi2)
543  {
544  const double TOL = Hermes::HermesSqrtEpsilon;
545 
546  // Newton Method
547  int local_nc;
548  double2* local_coeffs;
549  double2 local_lin_coeffs[H2D_MAX_NUMBER_VERTICES];
550  H1ShapesetJacobi shapeset;
551  short local_indices[70];
552 
553  // prepare the shapes and coefficients of the reference map
554  int j, k = 0;
555  for (unsigned char i = 0; i < e->get_nvert(); i++)
556  local_indices[k++] = shapeset.get_vertex_index(i, e->get_mode());
557 
558  // straight-edged element
559  if (e->cm == nullptr)
560  {
561  for (unsigned char i = 0; i < e->get_nvert(); i++)
562  {
563  local_lin_coeffs[i][0] = e->vn[i]->x;
564  local_lin_coeffs[i][1] = e->vn[i]->y;
565  }
566  local_coeffs = local_lin_coeffs;
567  local_nc = e->get_nvert();
568  }
569  else // curvilinear element - edge and bubble shapes
570  {
571  int o = e->cm->order;
572  for (unsigned char i = 0; i < e->get_nvert(); i++)
573  for (j = 2; j <= o; j++)
574  local_indices[k++] = shapeset.get_edge_index(i, 0, j, e->get_mode());
575 
576  if (e->is_quad())
577  o = H2D_MAKE_QUAD_ORDER(o, o);
578  memcpy(local_indices + k, shapeset.get_bubble_indices(o, e->get_mode()),
579  shapeset.get_num_bubbles(o, e->get_mode()) * sizeof(short));
580 
581  local_coeffs = e->cm->coeffs;
582  local_nc = e->cm->nc;
583  }
584 
585  // Constant reference mapping.
586  if (!e->is_curved() && (e->is_triangle() || e->is_parallelogram()))
587  {
588  double dx = e->vn[0]->x - x;
589  double dy = e->vn[0]->y - y;
590  int k = e->is_triangle() ? 2 : 3;
591  double m[2][2] =
592  {
593  { e->vn[1]->x - e->vn[0]->x, e->vn[k]->x - e->vn[0]->x },
594  { e->vn[1]->y - e->vn[0]->y, e->vn[k]->y - e->vn[0]->y }
595  };
596 
597  double const_jacobian = 0.25 * (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
598  double2x2 const_inv_ref_map;
599  if (const_jacobian <= 0.0)
600  throw Hermes::Exceptions::Exception("Element #%d is concave or badly oriented in RefMap::untransform().", e->id);
601 
602  double ij = 0.5 / const_jacobian;
603 
604  const_inv_ref_map[0][0] = m[1][1] * ij;
605  const_inv_ref_map[1][0] = -m[0][1] * ij;
606  const_inv_ref_map[0][1] = -m[1][0] * ij;
607  const_inv_ref_map[1][1] = m[0][0] * ij;
608 
609  xi1 = -1.0 - (const_inv_ref_map[0][0] * dx + const_inv_ref_map[1][0] * dy);
610  xi2 = -1.0 - (const_inv_ref_map[0][1] * dx + const_inv_ref_map[1][1] * dy);
611  }
612  else
613  {
614  double xi1_old = 0.0, xi2_old = 0.0;
615  double vx, vy;
616  double2x2 m;
617  // number of Newton iterations
618  int it = 0;
619  while (1)
620  {
621  double2x2 tmp;
622  memset(tmp, 0, sizeof(double2x2));
623  vx = vy = 0;
624  for (int i = 0; i < local_nc; i++)
625  {
626  double val = shapeset.get_fn_value(local_indices[i], xi1_old, xi2_old, 0, e->get_mode());
627  vx += local_coeffs[i][0] * val;
628  vy += local_coeffs[i][1] * val;
629 
630  double dx = shapeset.get_dx_value(local_indices[i], xi1_old, xi2_old, 0, e->get_mode());
631  double dy = shapeset.get_dy_value(local_indices[i], xi1_old, xi2_old, 0, e->get_mode());
632  tmp[0][0] += local_coeffs[i][0] * dx;
633  tmp[0][1] += local_coeffs[i][0] * dy;
634  tmp[1][0] += local_coeffs[i][1] * dx;
635  tmp[1][1] += local_coeffs[i][1] * dy;
636  }
637 
638  // inverse matrix
639  double jac = tmp[0][0] * tmp[1][1] - tmp[0][1] * tmp[1][0];
640  m[0][0] = tmp[1][1] / jac;
641  m[0][1] = -tmp[1][0] / jac;
642  m[1][0] = -tmp[0][1] / jac;
643  m[1][1] = tmp[0][0] / jac;
644 
645  xi1 = xi1_old - (m[0][0] * (vx - x) + m[1][0] * (vy - y));
646  xi2 = xi2_old - (m[0][1] * (vx - x) + m[1][1] * (vy - y));
647  if (fabs(xi1 - xi1_old) < TOL && fabs(xi2 - xi2_old) < TOL) return;
648  if (it > 1 && (xi1 > 1.5 || xi2 > 1.5 || xi1 < -1.5 || xi2 < -1.5)) return;
649  if (it > 100)
650  {
651  Hermes::Mixins::Loggable::Static::warn("Could not find reference coordinates - Newton method did not converge.");
652  return;
653  }
654  xi1_old = xi1;
655  xi2_old = xi2;
656  it++;
657  }
658  }
659  }
660 
661  static bool is_in_ref_domain(Element* e, double xi1, double xi2)
662  {
663  const double TOL = 1e-11;
664  if (e->get_nvert() == 3)
665  return (xi1 + xi2 <= TOL) && (xi1 + 1.0 >= -TOL) && (xi2 + 1.0 >= -TOL);
666  else
667  return (xi1 - 1.0 <= TOL) && (xi1 + 1.0 >= -TOL) && (xi2 - 1.0 <= TOL) && (xi2 + 1.0 >= -TOL);
668  }
669 
670  bool RefMap::is_element_on_physical_coordinates(Element* e, double x, double y, double* x_reference, double* y_reference)
671  {
672  // utility reference points that serve for the case when x_reference, y_reference are not passed.
673  double xi1, xi2;
674 
675  bool is_triangle = e->is_triangle();
676  bool is_curved = e->is_curved();
677 
678  if (is_curved)
679  {
680  untransform(e, x, y, xi1, xi2);
681  if (x_reference)
682  (*x_reference) = xi1;
683  if (y_reference)
684  (*y_reference) = xi2;
685 
686  if (is_in_ref_domain(e, xi1, xi2))
687  return true;
688  else
689  return false;
690  }
691 
692  // edge vectors.
693  double2 vector[4];
694  vector[0][0] = e->vn[1]->x - e->vn[0]->x;
695  vector[0][1] = e->vn[1]->y - e->vn[0]->y;
696  vector[1][0] = e->vn[2]->x - e->vn[1]->x;
697  vector[1][1] = e->vn[2]->y - e->vn[1]->y;
698  if (is_triangle)
699  {
700  vector[2][0] = e->vn[0]->x - e->vn[2]->x;
701  vector[2][1] = e->vn[0]->y - e->vn[2]->y;
702  }
703  else
704  {
705  vector[2][0] = e->vn[3]->x - e->vn[2]->x;
706  vector[2][1] = e->vn[3]->y - e->vn[2]->y;
707  vector[3][0] = e->vn[0]->x - e->vn[3]->x;
708  vector[3][1] = e->vn[0]->y - e->vn[3]->y;
709  }
710 
711  // calculate cross products
712  // -> if all cross products of edge vectors (vector[*]) x vector (thePoint - aVertex) are positive (negative),
713  // the point is inside of the element.
714  double cross_product_0 = (x - e->vn[0]->x) * vector[0][1] - (y - e->vn[0]->y) * vector[0][0];
715  double cross_product_1 = (x - e->vn[1]->x) * vector[1][1] - (y - e->vn[1]->y) * vector[1][0];
716  double cross_product_2 = (x - e->vn[2]->x) * vector[2][1] - (y - e->vn[2]->y) * vector[2][0];
717  if (is_triangle)
718  {
719  if ((cross_product_0 * cross_product_1 >= 0) && (cross_product_0 * cross_product_2 >= 0) && (cross_product_1 * cross_product_2 >= 0))
720  {
721  if (x_reference || y_reference)
722  {
723  untransform(e, x, y, xi1, xi2);
724  if (x_reference)
725  (*x_reference) = xi1;
726  if (y_reference)
727  (*y_reference) = xi2;
728  }
729  return true;
730  }
731  else
732  {
733  return false;
734  }
735  }
736  else
737  {
738  double cross_product_3 = (x - e->vn[3]->x) * vector[3][1] - (y - e->vn[3]->y) * vector[3][0];
739  if ((cross_product_0 * cross_product_1 >= 0) && (cross_product_0 * cross_product_2 >= 0) &&
740  (cross_product_1 * cross_product_2 >= 0) && (cross_product_1 * cross_product_3 >= 0) &&
741  (cross_product_2 * cross_product_3 >= 0) && (cross_product_0 * cross_product_3 >= 0))
742  {
743  if ((cross_product_0 * cross_product_1 >= 0) && (cross_product_0 * cross_product_2 >= 0) &&
744  (cross_product_1 * cross_product_2 >= 0) && (cross_product_1 * cross_product_3 >= 0) &&
745  (cross_product_2 * cross_product_3 >= 0) && (cross_product_0 * cross_product_3 >= 0))
746  {
747  if (x_reference || y_reference)
748  {
749  untransform(e, x, y, xi1, xi2);
750  if (x_reference)
751  (*x_reference) = xi1;
752  if (y_reference)
753  (*y_reference) = xi2;
754  }
755  return true;
756  }
757  else
758  return false;
759  }
760  }
761 
762  return false;
763  }
764 
765  Element* RefMap::element_on_physical_coordinates(bool use_MeshHashGrid, MeshSharedPtr mesh, double x, double y, double* x_reference, double* y_reference)
766  {
767  // utility element pointer.
768  Element *e;
769 
770  // utility reference points that serve for the case when x_reference, y_reference are not passed.
771  double xi1, xi2;
772 
773  // Optionally try the fastest approach for a multitude of successive calls - using the MeshHashGrid grid.
774  if (use_MeshHashGrid)
775  {
776  if (e = mesh->element_on_physical_coordinates(x, y))
777  {
778  if (x_reference || y_reference)
779  {
780  untransform(e, x, y, xi1, xi2);
781  if (x_reference)
782  (*x_reference) = xi1;
783  if (y_reference)
784  (*y_reference) = xi2;
785  }
786  return e;
787  }
788  }
789 
790  // vector for curved elements that do not contain the point when considering straightened edges.
791  // these are then checked at the end of this method, as they are really slow to look for the point in.
792  std::vector<Element*> improbable_curved_elements;
793 
794  // main loop over all active elements.
795  for_all_active_elements(e, mesh)
796  {
797  bool is_triangle = e->is_triangle();
798  bool is_curved = e->is_curved();
799 
800  // For curved elements.
801  bool is_near_straightened_element = false;
802 
803  // edge vectors.
804  double2 vector[4];
805  vector[0][0] = e->vn[1]->x - e->vn[0]->x;
806  vector[0][1] = e->vn[1]->y - e->vn[0]->y;
807  vector[1][0] = e->vn[2]->x - e->vn[1]->x;
808  vector[1][1] = e->vn[2]->y - e->vn[1]->y;
809  if (is_triangle)
810  {
811  vector[2][0] = e->vn[0]->x - e->vn[2]->x;
812  vector[2][1] = e->vn[0]->y - e->vn[2]->y;
813  }
814  else
815  {
816  vector[2][0] = e->vn[3]->x - e->vn[2]->x;
817  vector[2][1] = e->vn[3]->y - e->vn[2]->y;
818  vector[3][0] = e->vn[0]->x - e->vn[3]->x;
819  vector[3][1] = e->vn[0]->y - e->vn[3]->y;
820  }
821 
822  // calculate cross products
823  // -> if all cross products of edge vectors (vector[*]) x vector (thePoint - aVertex) are positive (negative),
824  // the point is inside of the element.
825  double cross_product_0 = (x - e->vn[0]->x) * vector[0][1] - (y - e->vn[0]->y) * vector[0][0];
826  double cross_product_1 = (x - e->vn[1]->x) * vector[1][1] - (y - e->vn[1]->y) * vector[1][0];
827  double cross_product_2 = (x - e->vn[2]->x) * vector[2][1] - (y - e->vn[2]->y) * vector[2][0];
828  if (is_triangle)
829  {
830  if ((cross_product_0 * cross_product_1 >= 0) && (cross_product_0 * cross_product_2 >= 0) && (cross_product_1 * cross_product_2 >= 0))
831  {
832  if (!is_curved)
833  {
834  if (x_reference || y_reference)
835  {
836  untransform(e, x, y, xi1, xi2);
837  if (x_reference)
838  (*x_reference) = xi1;
839  if (y_reference)
840  (*y_reference) = xi2;
841  }
842  return e;
843  }
844  else
845  is_near_straightened_element = true;
846  }
847  if (is_curved && !is_near_straightened_element)
848  {
849  double gravity_center_x = (e->vn[0]->x + e->vn[1]->x + e->vn[2]->x) / 3.;
850  double gravity_center_y = (e->vn[0]->y + e->vn[1]->y + e->vn[2]->y) / 3.;
851  double distance_x = std::abs(x - gravity_center_x);
852  double distance_y = std::abs(y - gravity_center_y);
853  if (distance_x < std::abs(e->vn[0]->x - gravity_center_x) && distance_x < std::abs(e->vn[1]->x - gravity_center_x) && distance_x < std::abs(e->vn[2]->x - gravity_center_x) &&
854  distance_y < std::abs(e->vn[0]->y - gravity_center_y) && distance_y < std::abs(e->vn[1]->y - gravity_center_y) && distance_y < std::abs(e->vn[2]->y - gravity_center_y))
855  is_near_straightened_element = true;
856  }
857  }
858  else
859  {
860  double cross_product_3 = (x - e->vn[3]->x) * vector[3][1] - (y - e->vn[3]->y) * vector[3][0];
861  if ((cross_product_0 * cross_product_1 >= 0) && (cross_product_0 * cross_product_2 >= 0) &&
862  (cross_product_1 * cross_product_2 >= 0) && (cross_product_1 * cross_product_3 >= 0) &&
863  (cross_product_2 * cross_product_3 >= 0) && (cross_product_0 * cross_product_3 >= 0))
864  {
865  if (!is_curved)
866  {
867  if (x_reference || y_reference)
868  {
869  untransform(e, x, y, xi1, xi2);
870  if (x_reference)
871  (*x_reference) = xi1;
872  if (y_reference)
873  (*y_reference) = xi2;
874  }
875  return e;
876  }
877  else
878  is_near_straightened_element = true;
879  }
880  if (is_curved && !is_near_straightened_element)
881  {
882  double gravity_center_x = (e->vn[0]->x + e->vn[1]->x + e->vn[2]->x + e->vn[3]->x) / 4.;
883  double gravity_center_y = (e->vn[0]->y + e->vn[1]->y + e->vn[2]->y + e->vn[3]->y) / 4.;
884  double distance_x = std::abs(x - gravity_center_x);
885  double distance_y = std::abs(y - gravity_center_y);
886  if (distance_x <= std::abs(e->vn[0]->x - gravity_center_x) && distance_x <= std::abs(e->vn[1]->x - gravity_center_x) && distance_x <= std::abs(e->vn[2]->x - gravity_center_x) && distance_x <= std::abs(e->vn[3]->x - gravity_center_x) &&
887  distance_y <= std::abs(e->vn[0]->y - gravity_center_y) && distance_y <= std::abs(e->vn[1]->y - gravity_center_y) && distance_y <= std::abs(e->vn[2]->y - gravity_center_y) && distance_y <= std::abs(e->vn[3]->y - gravity_center_y))
888  is_near_straightened_element = true;
889  }
890  }
891 
892  if (is_curved)
893  {
894  if (is_near_straightened_element)
895  {
896  untransform(e, x, y, xi1, xi2);
897  if (is_in_ref_domain(e, xi1, xi2))
898  {
899  if (x_reference != nullptr)
900  (*x_reference) = xi1;
901  if (y_reference != nullptr)
902  (*y_reference) = xi2;
903  return e;
904  }
905  }
906  else
907  improbable_curved_elements.push_back(e);
908  }
909  }
910 
911  // loop through the improbable curved elements.
912  for (unsigned short i = 0; i < improbable_curved_elements.size(); i++)
913  {
914  untransform(improbable_curved_elements[i], x, y, xi1, xi2);
915  if (is_in_ref_domain(improbable_curved_elements[i], xi1, xi2))
916  {
917  if (x_reference != nullptr)
918  (*x_reference) = xi1;
919  if (y_reference != nullptr)
920  (*y_reference) = xi2;
921  return improbable_curved_elements[i];
922  }
923  }
924 
925  Hermes::Mixins::Loggable::Static::warn("Point (%g, %g) does not lie in any element.", x, y);
926  return nullptr;
927  }
928 
929  void RefMap::force_transform(uint64_t sub_idx, Trf* ctm)
930  {
931  this->sub_idx = sub_idx;
932  stack[top] = *ctm;
933  this->ctm = stack + top;
934  if (this->is_const)
935  calc_const_inv_ref_map();
936  this->reinit_storage();
937  }
938  }
939 }
Hermes::Hermes2D::Node::x
double x
vertex node coordinates
Definition: element.h:63
Hermes::Hermes2D::Shapeset::get_edge_index
short get_edge_index(unsigned char edge, unsigned short ori, unsigned short order, ElementMode2D mode) const
Definition: shapeset.cpp:208
Hermes
Definition: adapt.h:24
Hermes::Hermes2D::Element::next_vert
unsigned char next_vert(unsigned char i) const
Helper functions to obtain the index of the next or previous vertex/edge.
Definition: element.h:202
Hermes::Hermes2D::Element::id
int id
element id number
Definition: element.h:112
Hermes::Hermes2D::Element
Stores one element of a mesh.
Definition: element.h:107
Hermes::Hermes2D::H1ShapesetJacobi
Shape functions based on integrated Jacobi polynomials.
Definition: shapeset_h1_all.h:27
Hermes::Hermes2D::Element::cm
CurvMap * cm
curved mapping, nullptr if not curvilinear
Definition: element.h:188
Hermes::Hermes2D::CurvMap::order
unsigned short order
current polynomial degree of the refmap approximation
Definition: curved.h:122
global.h
Common definitions for Hermes2D.
Hermes::Hermes2D::Element::iro_cache
unsigned short iro_cache
Increase in integration order, see RefMap::calc_inv_ref_order()
Definition: element.h:199
Hermes::Hermes2D::Shapeset::get_vertex_index
short get_vertex_index(int vertex, ElementMode2D mode) const
Returns the index of a vertex shape function associated with the specified vertex.
Definition: shapeset.cpp:197
Hermes::Hermes2D::Shapeset::get_ref_vertex
double2 * get_ref_vertex(int vertex, ElementMode2D mode)
Returns the coordinates of the reference domain vertices.
Definition: shapeset.cpp:360
Hermes::Hermes2D::Trf
2D transformation.
Definition: transformable.h:29
H2D_MAX_NUMBER_VERTICES
#define H2D_MAX_NUMBER_VERTICES
A maximum number of vertices of an element.
Definition: global.h:32
Hermes::Hermes2D::Shapeset::get_num_bubbles
virtual unsigned short get_num_bubbles(unsigned short order, ElementMode2D mode) const
Returns the number of bubble functions for an element of the given order.
Definition: shapeset.cpp:233
Hermes::Hermes2D::Shapeset::get_bubble_indices
virtual short * get_bubble_indices(unsigned short order, ElementMode2D mode) const
Returns a complete set of indices of bubble functions for an element of the given order...
Definition: shapeset.cpp:221
Hermes::Hermes2D::Element::vn
Node * vn[H2D_MAX_NUMBER_VERTICES]
vertex node pointers
Definition: element.h:134