Dealii: Incorrect result from Mapping::transform_real_to_unit_cell

[jppelteret]

This relates to the conversation in this forum post. When a regular Cartesian-aligned grid constructed using GridGenerator::subdivided_hyper_rectangle with a non-trivial number of repetitions is refined using refine_global(), for a certain subset of cells Mapping::transform_real_to_unit_cell returns an inaccurately computed point inside the unit cell. Various other tested combinations of grid creation and refinement do not necessary demonstrate a similar pathology, but according to the initial post in the forum thread this can affect read-in grids as well.

An expanded example that demonstrates the issue is attached

[drwells]

These cells are relatively large and not distorted. This is a bit scary.

I'll try to find some time to look at this tomorrow.

[jppelteret]

Thanks @drwells . I also find this is quite alarming!

[drwells]

I lied when I said I would look into it tomorrow. I couldn't help myself and wrote a draft this morning.

The problem is that with a Q1 mapping and parallelograms (not just on a Cartesian grid), in the glorified quadratic equation solver, that both a = eps (purely due to roundoff) and b - sqrt(discriminant) = eps (due to the size of a). However, the current logic only checks the relative size of c, not the relative size of a: hence we are toast if a is small but nonzero and the other terms are O(1).

I think that the best explanation is that, for parallelograms, we don't really need a bilinear mapping so the second-order terms in the mapping go away. I will also change the default to assume that a is small.

I will post a patch in a few minutes.

Edit: somehow I posted another draft above this. That message is deleted now.

[jppelteret]

Fixed by #4659