Hi Andreas,
Unfortunately, I believe that due to roundoff these functions will always
produce questionable results when points are near the cell's boundary. If I
recall correctly, at least for 2D the errors tended to be a couple of
multiples larger than machine precision when doing these sorts of things,
so the library functions won't help much.
One possibility would be, if you have a very easy geometry, to try
MappingCartesian instead. That will have fewer roundoff errors than
MappingQ1/MappingQGeneric.
If you know that the points are close to a cell face then another possible
work around would be to project that point onto the face and then see if
the vector from the point to the projected point faces in or out of the
cell.
I hope that one of these helps.
Thanks,
David Wells
On Thursday, May 19, 2016 at 9:47:36 AM UTC-4, Andreas Krämer wrote:
>
> 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.
