Use boost multiprecision for GJK

[chrisrichardson]

No description provided.

[jorgensd]

Have you tested this on: #3653

In general, it seems like you have rewritten the algorithm quite a bit. Could you comment on why and how it improves the situation (of course higher precision helps).

[chrisrichardson]

Yes it works, the answer is:
[-1.03648825 -1.01316152 1.08986529]

I tried to make the algorithm clearer - there are some comments throughout. It is possible there is more simplification that can be done.

[jhale]

Looks nice. Multiprecision is header only and we already require boost so that's good.

Do we want the caller to be able to control the internal precision?

[chrisrichardson]

It would be easy to add a template parameter for the internal precision, good for c++, but probably not so much for Python. We could use presets for double and float, I guess.

The functions are already in impl_gjk namespace. Some are duplicated from dolfinx/math.h (e.g. det3) because fma is used there which is not compatible with multiprecision libraries.

[chrisrichardson]

On oneapi seems to be dispatching to float32 incorrectly...

[garth-wells]

It would be easy to add a template parameter for the internal precision, good for c++, but probably not so much for Python. We could use presets for double and float, I guess.

The functions are already in impl_gjk namespace. Some are duplicated from dolfinx/math.h (e.g. det3) because fma is used there which is not compatible with multiprecision libraries.

det3 and friends could be programmed to handled different types to avoid duplication of interfaces.

[jorgensd]

I think we should add a test covering the co-planar case (#3653) so that we don't hit this again, other than that, looks good.

[garth-wells]

This feels a little bit arbitrary. Do we really understand where and why extra precision is needed?

Are we sure that the nanobind interface wasn't casting to lower precision?

What's the performance cost? It could be high for the extended precision.

[chrisrichardson]

To answer those questions:

• The algorithm often enters a limit cycle, and does not converge due to insufficient accuracy. Another alternative would be to reduce the tolerance to be much looser. Note that many other libraries, like CGAL, use high precision for geometry calculations, presumably for the same reasons.
• The nanobind casting issue was only on some platforms. I am 100% certain this is not the main issue.
• Of course, the multiprecision is slower, maybe up to 4 times slower. Another alternative would be to call the interface with double and only go to the extended precision if it fails to converge.

changes made

[jhale]

We should look at unifying the math functions in another PR - with C++ concepts it should be straightforward.

It's not clear to me that we would want to handle increasing precision within the function - shouldn't this be under the control of the caller handling an exception/error code?

[chrisrichardson]

Yes, I agree with the second point. I would like to merge this for now - we can create some actions/issues to address these remaining items. Perhaps it should return a value to indicate (non)convergence, rather than throwing an error.

[jhale]

Probably the latter because I'd consider non-convergence 'normal' rather than exceptional behaviour. Unfortunately we don't have std::expected in C++20 (C++23).

[garth-wells]

Yes, I agree with the second point. I would like to merge this for now - we can create some actions/issues to address these remaining items. Perhaps it should return a value to indicate (non)convergence, rather than throwing an error.

In principle fine, but in practice the details never get fixed later.