Ensure to look at square of tolerances

[kronbichler]

@bangerth @gassmoeller regarding the tolerance topic I checked the numbers again and believe we should have squares. The reason is that we take the difference between two vertices (each component of which we expect to either satisfy a relative tolerance against the cell diameter or against the absolute tolerance), and then the square in square() or square_distance(). As a result, the numbers should really be on the order of machine_epsilon^2. I looked into prototype numbers of the test that was added in #14092 and can confirm that numbers are really in that order of magnitude: We see squares of around 1e-19, compared to numbers which are 1e13, all of them being squares of either size 1e6 or 1e-9, the roundoff accuracy around 1e6.

[kronbichler]

/rebuild

[gassmoeller]

Please wait before merging this. It seems the test numbers support your interpretation, but I would like to ask your opinion on two questions first:

1. When I consider the problem of how large an tolerance in the difference of two squares is, if the original numbers differ by 1e-16, I would think it works the following way:

x = 1.;
y = 1. + 1e-16;

x^2 = 1;
y^2 = (1^2 + 2e-16 + 1e-32)

y^2 - x^2 = 2e-16 + 1e-32

Did I miss something?

2. Even if my estimate above is not correct, are we not putting too much energy and discussion time in this sanity check? The assert is there to prevent a non-cartesian cell to be used with a MappingCartesian. The worst that can happen if we make the tolerance too lax is that some models with cells that are nearly squares but not quite (with a tolerance between 1e-7 and 1e-14) are still classified as cartesian, even if they are not quite. But this check is there to prevent models with clearly un-cartesian cells (say in a hypershell) to use a MappingCartesian. I.e. even if we do not catch the nearly-cartesian cells, there will be cells that are stronger deformed in an arbitrary MappingQ or whatever is used. So the check still works with a lax tolerance, while it breaks complete models if it is too strict.

I would be in favor of erring on the side of a lax check. Then at least only models that make mistakes crash, and not models that use the code as intended.

[bangerth]

On second thought, I think you are right. I had had in my head that we are comparing

  if (x*x == y*y)

with a tolerance, for which the tolerance should be 1e-16*(x*x). But in reality, we are comparing

  if ( (x-y)*(x-y) == 0)

and you are right that the tolerance in that case should be (1e-16*abs(x))**2 or something proportional.

I'll let @gassmoeller finish his thought, though.

[kronbichler]

Indeed, I agree with what @bangerth wrote. I would be happy with a lax check as well. If we are fine with tolerances around 1e-7/1e-8, I would be happy to just add the comment to the commit and keep the present numbers.

[gassmoeller]

You are right @bangerth, we are computing the square of the difference, not the difference of the square, my mistake. Ok, then go ahead. I still do not think it will matter in the end, but I am at least reasonably certain we are not going to break anything with this PR.