Hi,
I have just encountered an issue with projecting points into a cell.
Say, I have a point x which is very close to a cell (e.g. distance =
10^-16) and I want this point to be inside the cell. My first attempt was
to map the point to the unit cell, round the coordinates to 0 or 1, and map
it back to the real cell. However, the mapping introduces an additional
round-off error and cell->point_inside(x) is still false. (I use
MappingQ(1))
As a first workaround, I moved the mapped point into the unit cell by 1e-12.
Is there a better way to do it?
Best,
Andreas
A code snippet is below (the last assertion gives false)
// eliminate round-off errors
dealii::Point<dim> h = pi_unit + increment;
for (size_t i = 0; i < dim; i++) {
if (fabs(h[i]) < 1e-12) {
h[i] = 0;
}
if (fabs(h[i] - 1) < 1e-12) {
h[i] = 1;
}
}
// assert that point is at boundary
if (2 == dim) {
assert(h[0] * h[1] * (1 - h[0]) * (1 - h[1]) == 0);
} else { // 3 == dim
assert(h[0] * h[1] * h[2] * (1 - h[0]) * (1 - h[1]) * (1 - h[2]) ==
0);
}
// map to real cell
p_boundary = m_mapping.transform_unit_to_real_cell(ci, h);
cout << "boundary point: " << p_boundary << endl;
assert(cell->point_inside(p_boundary));
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
