24 RefMap::RefMap() : ref_map_shapeset(H1ShapesetJacobi()), ref_map_pss(&ref_map_shapeset)
27 set_quad_2d(&g_quad_2d_std);
28 this->reinit_storage();
38 double* RefMap::get_phys_x(
int order)
40 if (order != this->phys_x_calculated)
41 this->calc_phys_x(order);
45 double* RefMap::get_phys_y(
int order)
47 if (order != this->phys_y_calculated)
48 this->calc_phys_y(order);
52 double3* RefMap::get_tangent(
int edge,
int order)
55 order = quad_2d->get_edge_points(edge, quad_2d->get_max_order(element->get_mode()), element->get_mode());
57 if (order != this->tan_calculated[edge])
58 this->calc_tangent(edge, order);
63 void RefMap::set_quad_2d(
Quad2D* quad_2d)
65 this->quad_2d = quad_2d;
66 ref_map_pss.set_quad_2d(quad_2d);
67 this->reinit_storage();
70 void RefMap::set_active_element(
Element* e)
72 this->reinit_storage();
74 Transformable::set_active_element(e);
75 ref_map_pss.set_active_element(e);
77 this->is_const = element->has_const_ref_map();
80 unsigned short j, k = 0;
81 for (
unsigned char i = 0; i < e->get_nvert(); i++)
87 for (
unsigned char i = 0; i < e->get_nvert(); i++)
89 lin_coeffs[i][0] = e->
vn[i]->
x;
90 lin_coeffs[i][1] = e->
vn[i]->y;
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());
102 if (e->is_quad()) o = H2D_MAKE_QUAD_ORDER(o, o);
106 coeffs = e->
cm->coeffs;
110 this->inv_ref_order = this->element->iro_cache;
114 calc_const_inv_ref_map();
116 const_jacobian = 0.0;
119 void RefMap::set_element_iro_cache(
Element* element)
121 bool is_const = element->has_const_ref_map();
130 #pragma omp critical (element_iro_cache_setting)
133 rm.set_active_element(element);
134 element->
iro_cache = rm.calc_inv_ref_order();
138 void RefMap::reinit_storage()
140 jacobian_calculated = -1;
141 inv_ref_map_calculated = -1;
142 second_ref_map_calculated = -1;
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;
150 void RefMap::push_transform(
int son)
152 Transformable::push_transform(son);
153 const_jacobian *= 0.25;
154 this->reinit_storage();
157 void RefMap::pop_transform()
159 Transformable::pop_transform();
160 const_jacobian *= 4.;
161 this->reinit_storage();
164 void RefMap::calc_inv_ref_map(
int order)
166 int i, j, np = quad_2d->get_num_points(order, element->get_mode());
169 ref_map_pss.force_transform(sub_idx, ctm);
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];
176 for (j = 0; j < np; j++)
177 m_00[j] = m_01[j] = m_10[j] = m_11[j] = 0.;
179 for (i = 0; i < nc; i++)
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++)
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_;
200 double trj = get_transform_jacobian();
201 double2x2* irm = this->inv_ref_map;
202 double* jac = this->jacobian;
204 for (i = 0; i < np; i++)
206 jac[i] = (m_00[i] * m_11[i] - m_01[i] * m_10[i]);
207 double ij = 1.0 / jac[i];
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);
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);
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;
223 this->inv_ref_map_calculated = order;
224 this->jacobian_calculated = order;
227 void RefMap::calc_const_inv_ref_map()
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 } };
233 const_jacobian = 0.25 * (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
235 double ij = 0.5 / const_jacobian;
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;
242 const_jacobian *= get_transform_jacobian();
245 void RefMap::calc_phys_x(
int order)
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++)
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];
260 this->phys_x_calculated = order;
263 void RefMap::calc_phys_y(
int order)
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++)
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];
278 this->phys_y_calculated = order;
281 void RefMap::calc_tangent(
int edge,
int eo)
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);
288 if (!element->is_curved())
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];
297 tan[0][2] *= (edge == 0 || edge == 2) ? ctm->m[0] : ctm->m[1];
299 for (i = 1; i < np; i++)
300 memcpy(tan + i, tan,
sizeof(double3));
307 memset(m, 0, np*
sizeof(double2x2));
308 ref_map_pss.force_transform(sub_idx, ctm);
309 for (i = 0; i < nc; i++)
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++)
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];
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());
328 double ex = (*v2)[0] - (*v1)[0];
329 double ey = (*v2)[1] - (*v1)[1];
330 for (i = 0; i < np; 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];
339 t[2] *= (edge == 0 || edge == 2) ? ctm->m[0] : ctm->m[1];
343 this->tan_calculated[edge] = eo;
346 int RefMap::calc_inv_ref_order()
348 Quad2D* quad = get_quad_2d();
349 int i, o, mo = (quad->get_max_order(element->get_mode()) / 2);
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++)
357 throw Hermes::Exceptions::Exception(
"Element #%d is concave or badly oriented.", element->
id);
363 for (i = 0; i < quad->get_num_points(mo, element->get_mode()); i++, m++)
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];
369 for (o = 0; o < mo; o++)
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++)
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];
381 if ((fabs((exact1 - result1) / exact1) < Hermes::HermesSqrtEpsilon) && (fabs((exact2 - result2) / exact2) < Hermes::HermesSqrtEpsilon))
386 this->warn(
"Element #%d is too distorted (iro ~ %d).", element->
id, o);
391 void RefMap::inv_ref_map_at_point(
double xi1,
double xi2,
double& x,
double& y, double2x2& m)
394 memset(tmp, 0,
sizeof(double2x2));
396 for (
int i = 0; i < nc; i++)
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;
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;
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;
418 #ifdef H2D_USE_SECOND_DERIVATIVES
419 void RefMap::calc_second_ref_map(
int order)
421 int i, j, np = quad_2d->get_num_points(order, element->get_mode());
423 double3x2* k = calloc_with_check<double3x2>(np);
424 ref_map_pss.force_transform(sub_idx, ctm);
425 for (i = 0; i < nc; i++)
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++)
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];
443 double3x2* mm = this->second_ref_map;
444 double2x2* m = get_inv_ref_map(order);
445 for (j = 0; j < np; j++)
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];
452 mm[j][0][0] = -(a * m[j][0][0] + b * m[j][1][0]);
454 mm[j][0][1] = -(a * m[j][0][1] + b * m[j][1][1]);
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];
460 mm[j][1][0] = -(a * m[j][0][0] + b * m[j][1][0]);
462 mm[j][1][1] = -(a * m[j][0][1] + b * m[j][1][1]);
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];
468 mm[j][2][0] = -(a * m[j][0][0] + b * m[j][1][0]);
470 mm[j][2][1] = -(a * m[j][0][1] + b * m[j][1][1]);
475 this->second_ref_map_calculated = order;
478 void RefMap::second_ref_map_at_point(
double xi1,
double xi2,
double& x,
double& y, double3x2& mm)
this-
481 memset(k, 0,
sizeof(double3x2));
483 for (
int i = 0; i < nc; i++)
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;
489 double dxy, dxx, dyy;
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());
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;
504 this->inv_ref_map_at_point(xi1, xi2, x, y, m);
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];
511 mm[0][0] = -(a * m[0][0] + b * m[1][0]);
513 mm[0][1] = -(a * m[0][1] + b * m[1][1]);
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];
519 mm[1][0] = -(a * m[0][0] + b * m[1][0]);
521 mm[1][1] = -(a * m[0][1] + b * m[1][1]);
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];
527 mm[2][0] = -(a * m[0][0] + b * m[1][0]);
529 mm[2][1] = -(a * m[0][1] + b * m[1][1]);
532 double3x2* RefMap::get_second_ref_map(
int order)
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;
542 void RefMap::untransform(
Element* e,
double x,
double y,
double& xi1,
double& xi2)
544 const double TOL = Hermes::HermesSqrtEpsilon;
548 double2* local_coeffs;
551 short local_indices[70];
555 for (
unsigned char i = 0; i < e->get_nvert(); i++)
559 if (e->
cm ==
nullptr)
561 for (
unsigned char i = 0; i < e->get_nvert(); i++)
563 local_lin_coeffs[i][0] = e->
vn[i]->
x;
564 local_lin_coeffs[i][1] = e->
vn[i]->y;
566 local_coeffs = local_lin_coeffs;
567 local_nc = e->get_nvert();
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());
577 o = H2D_MAKE_QUAD_ORDER(o, o);
581 local_coeffs = e->
cm->coeffs;
582 local_nc = e->
cm->nc;
586 if (!e->is_curved() && (e->is_triangle() || e->is_parallelogram()))
588 double dx = e->
vn[0]->
x - x;
589 double dy = e->
vn[0]->y - y;
590 int k = e->is_triangle() ? 2 : 3;
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 }
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);
602 double ij = 0.5 / const_jacobian;
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;
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);
614 double xi1_old = 0.0, xi2_old = 0.0;
622 memset(tmp, 0,
sizeof(double2x2));
624 for (
int i = 0; i < local_nc; i++)
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;
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;
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;
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;
651 Hermes::Mixins::Loggable::Static::warn(
"Could not find reference coordinates - Newton method did not converge.");
661 static bool is_in_ref_domain(
Element* e,
double xi1,
double xi2)
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);
667 return (xi1 - 1.0 <= TOL) && (xi1 + 1.0 >= -TOL) && (xi2 - 1.0 <= TOL) && (xi2 + 1.0 >= -TOL);
670 bool RefMap::is_element_on_physical_coordinates(
Element* e,
double x,
double y,
double* x_reference,
double* y_reference)
675 bool is_triangle = e->is_triangle();
676 bool is_curved = e->is_curved();
680 untransform(e, x, y, xi1, xi2);
682 (*x_reference) = xi1;
684 (*y_reference) = xi2;
686 if (is_in_ref_domain(e, xi1, xi2))
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;
700 vector[2][0] = e->
vn[0]->
x - e->
vn[2]->
x;
701 vector[2][1] = e->
vn[0]->y - e->
vn[2]->y;
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;
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];
719 if ((cross_product_0 * cross_product_1 >= 0) && (cross_product_0 * cross_product_2 >= 0) && (cross_product_1 * cross_product_2 >= 0))
721 if (x_reference || y_reference)
723 untransform(e, x, y, xi1, xi2);
725 (*x_reference) = xi1;
727 (*y_reference) = xi2;
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))
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))
747 if (x_reference || y_reference)
749 untransform(e, x, y, xi1, xi2);
751 (*x_reference) = xi1;
753 (*y_reference) = xi2;
765 Element* RefMap::element_on_physical_coordinates(
bool use_MeshHashGrid, MeshSharedPtr mesh,
double x,
double y,
double* x_reference,
double* y_reference)
774 if (use_MeshHashGrid)
776 if (e = mesh->element_on_physical_coordinates(x, y))
778 if (x_reference || y_reference)
780 untransform(e, x, y, xi1, xi2);
782 (*x_reference) = xi1;
784 (*y_reference) = xi2;
792 std::vector<Element*> improbable_curved_elements;
795 for_all_active_elements(e, mesh)
797 bool is_triangle = e->is_triangle();
798 bool is_curved = e->is_curved();
801 bool is_near_straightened_element =
false;
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;
811 vector[2][0] = e->
vn[0]->
x - e->
vn[2]->
x;
812 vector[2][1] = e->
vn[0]->y - e->
vn[2]->y;
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;
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];
830 if ((cross_product_0 * cross_product_1 >= 0) && (cross_product_0 * cross_product_2 >= 0) && (cross_product_1 * cross_product_2 >= 0))
834 if (x_reference || y_reference)
836 untransform(e, x, y, xi1, xi2);
838 (*x_reference) = xi1;
840 (*y_reference) = xi2;
845 is_near_straightened_element =
true;
847 if (is_curved && !is_near_straightened_element)
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;
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))
867 if (x_reference || y_reference)
869 untransform(e, x, y, xi1, xi2);
871 (*x_reference) = xi1;
873 (*y_reference) = xi2;
878 is_near_straightened_element =
true;
880 if (is_curved && !is_near_straightened_element)
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;
894 if (is_near_straightened_element)
896 untransform(e, x, y, xi1, xi2);
897 if (is_in_ref_domain(e, xi1, xi2))
899 if (x_reference !=
nullptr)
900 (*x_reference) = xi1;
901 if (y_reference !=
nullptr)
902 (*y_reference) = xi2;
907 improbable_curved_elements.push_back(e);
912 for (
unsigned short i = 0; i < improbable_curved_elements.size(); i++)
914 untransform(improbable_curved_elements[i], x, y, xi1, xi2);
915 if (is_in_ref_domain(improbable_curved_elements[i], xi1, xi2))
917 if (x_reference !=
nullptr)
918 (*x_reference) = xi1;
919 if (y_reference !=
nullptr)
920 (*y_reference) = xi2;
921 return improbable_curved_elements[i];
925 Hermes::Mixins::Loggable::Static::warn(
"Point (%g, %g) does not lie in any element.", x, y);
929 void RefMap::force_transform(uint64_t sub_idx,
Trf* ctm)
931 this->sub_idx = sub_idx;
933 this->ctm = stack + top;
935 calc_const_inv_ref_map();
936 this->reinit_storage();
double x
vertex node coordinates
short get_edge_index(unsigned char edge, unsigned short ori, unsigned short order, ElementMode2D mode) const
unsigned char next_vert(unsigned char i) const
Helper functions to obtain the index of the next or previous vertex/edge.
Stores one element of a mesh.
Shape functions based on integrated Jacobi polynomials.
CurvMap * cm
curved mapping, nullptr if not curvilinear
unsigned short order
current polynomial degree of the refmap approximation
Common definitions for Hermes2D.
unsigned short iro_cache
Increase in integration order, see RefMap::calc_inv_ref_order()
short get_vertex_index(int vertex, ElementMode2D mode) const
Returns the index of a vertex shape function associated with the specified vertex.
double2 * get_ref_vertex(int vertex, ElementMode2D mode)
Returns the coordinates of the reference domain vertices.
#define H2D_MAX_NUMBER_VERTICES
A maximum number of vertices of an element.
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.
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...
Node * vn[H2D_MAX_NUMBER_VERTICES]
vertex node pointers