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