[drwells]

Rounding issues in the quadratic formula have bitten me several times: in particular, the inverse of a 2D bilinear mapping requires solving a quadratic equation which (if things are skewed) can be very badly conditioned. On the other hand, if one uses a bilinear map on a parallelogram then the quadratic term is gone and one root is zero: also annoying.

In practice I just check the coefficients and then use appropriate versions of the quadratic formula. I guess that I am luckier than Mr. Cook here since I only need one root.

code (some of it was written by me, some improvements are by other people): https://github.com/dealii/dealii/blob/master/source/fe/mappi...

This set of cases (which can probably be further improved) is equivalent to his first answer. I don't think its possible to do better in this test case due to the catastrophic cancelation in subtracting the discriminant from b.

[obl]

inverse bilinear map is a good example where, even if there is an analytical solution, it's often more robust to just do a couple newton iterations until you get down to machine precision.

added benefit is that it carries on in 3D (although with a larger jacobian of course).

[drwells]

You are absolutely right re: robustness (and that code falls back to the iterative method), but IIRC the Newton scheme was at least one order of magnitude slower than the formula. It was also slightly less accurate simply because it did more floating point operations.

I have not yet figured out how to make the 3D formula work so we always do the iterative method in 3D.

A nice trick: we approximate the (bi/tri)linear map as affine (usually quite accurate) to get a good initial guess.