Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add tolerance to box in compute_point_locations_try_all #15102

Merged

Conversation

kronbichler
Copy link
Member

@kronbichler kronbichler commented Apr 14, 2023

In #15035, we observed a test failure of numerics/generalized_interpolation_02 because of the following assertion:

// Ensure that we found all points
AssertThrow(std::get<3>(cell_qpoint_map).empty(),
VectorTools::ExcPointNotAvailableHere());

What is happening there is that a point ends up being outside the triangulation by roundoff (1e-16), which can happen during a typical program run. Looking into the code, I realize the following condition
for (const auto &point_and_id :
p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
bgi::intersects(box)))
performs a check intersects(box) without tolerance. Since I don't know how to teach boost a tolerance (maybe you do @luca-heltai?), I think the right way is to make the bounding box we accept larger by a roundoff margin, say 1e-12 relative tolerance. I added a test that checks this for the actual numbers where it happened (outside by 1e-16), by outside of 1e-8, and outside a lot. I think this way the code makes more sense.

const auto &box = leaf.first;
const double relative_tolerance = 1e-12;
const BoundingBox<spacedim> box = leaf.first.create_extended(
relative_tolerance * leaf.first.side_length(0));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be better if we scale each side length with a relative tolerance. Thanks for tackling this! @kronbichler

Copy link
Member

@peterrum peterrum Apr 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about creating BB::create_extended_relative()? @bergbauer could write this ;)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added @bergbauer's suggestion (via his own commit) to this PR.

Copy link
Member

@peterrum peterrum left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me 👍 I am wondering when we will have eliminated all the tolerance issues!?

const auto &box = leaf.first;
const double relative_tolerance = 1e-12;
const BoundingBox<spacedim> box = leaf.first.create_extended(
relative_tolerance * leaf.first.side_length(0));
Copy link
Member

@peterrum peterrum Apr 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about creating BB::create_extended_relative()? @bergbauer could write this ;)

@kronbichler kronbichler force-pushed the add_tolerance_compute_point_locations branch from e2c4bf3 to cce10ff Compare April 14, 2023 13:58
@bergbauer
Copy link
Contributor

Can we have another review on this? #15035 depends on this to fix the failing tests there.

Copy link
Member

@peterrum peterrum left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If there are no further comments, I will merge this tonight so that we can proceed with #15035.

@peterrum peterrum merged commit 0214143 into dealii:master Apr 17, 2023
14 checks passed
@kronbichler kronbichler deleted the add_tolerance_compute_point_locations branch August 10, 2023 16:39
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants