Mapping points of mesh generated by gmsh

[Giovanni Di Ilio]

Dear Deal.II community,

I am facing a problem with the mapping of a point on the real cell to the unit cell when the mesh is read in from a .msh file.

So what I have is a .geo file (which contains "physical lines" and "physical surfaces", as indicated in Step49). I generate the .msh file and I read in using the GridIn class. The result is a uniform structured mesh.

Then, since I want to compute the shape functions, I need to know the position of the mapped points on the unit cell. What happens is that for some cell (only few cells) the function transform_real_to_unit_cell gives me a wrong result, that is the coordinates of some mapped point are not correct.
If I generate the same identical mesh using subdivided_hyper_rectangle instead, everything works well.

I checked the .msh file and everything looks fine for me. What concerns me the most is that this occurs only for certain "random" cells of the domain.
I don't see any apparent explanation to this. For sure I am missing something but I was not able to figure out the solution to this issue.

I provide a small code that I am testing as well as the .msh and .geo files I am using.
The code computes the mapped point of a generic real point which is placed in the center of a "sick" cell (id=354). I would expect that the coordinates of the unit point are x=0.5, y=0.5. However, what I get is x=0.5, y=0.0. The same apply if you pick an other real point within this cell.

May you please help me? I would really appreciate it.
Thank you very much in advance,

Best,
Giovanni

[Wolfgang Bangerth]

On 07/20/2017 05:53 AM, Giovanni Di Ilio wrote:
>
> I provide a small code that I am testing as well as the .msh and .geo files I
> am using.
> The code computes the mapped point of a generic real point which is placed in
> the center of a "sick" cell (id=354). I would expect that the coordinates of
> the unit point are x=0.5, y=0.5. However, what I get is x=0.5, y=0.0. The same
> apply if you pick an other real point within this cell.

It's often difficult to debug these cases with such a large mesh. But since
you know which cell the problem appears in, take the test.msh file and remove
(by hand) all other cells and all of the vertices not needed. This way you
should end up with a testcase that has only one cell and that should be much
easier to figure out.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Giovanni Di Ilio]

Dear Wolfgang,

thanks for your quick reply and for your suggestion.
I simplified the case to only two cells, the one that was supposed to be "wrong" and one of its neighbours, that was supposed to be ok, for comparison.
By watching at the coordinates of the vertices in the file .msh I found out that the problem was related to some rounding error. I guess there is a bug in the parser. I just converted the coordinate values into double and now the mapping works fine everywhere.

Thank you very much again,

Best,
Giovanni

[Jean-Paul Pelteret]

Dear Giovanni,

Its good to hear that you managed to solve / work around this issue! Would it be possible for you to post the reduced (problematic) mesh so that we can see if there's anything that we can learn from it. I'm guessing that you expected lines with entries such as

6 0.1999999999991087 0 0

to be parsed as

6 0.2 0 0

? Did you ultimately have to output a grid that had the coordinates recorded with a greater precision?

Thanks,

Jean-Paul

[Giovanni Di Ilio]

Dear Jean-Paul,

yes, exactly, if the entries were parsed as you mentioned everything would work fine.
But unfortunately I guess this is not the case. I post here the simplified problem.
Now the mesh described in the .msh file is made of two cells: the one with index 0 which is wrong and one of its neighbours, with index 1, which is ok.
The vertexes of cell 0 are the following:

x1 = 5.399999999999241 , y1 = 0.4999999999988218
x2 = 5.499999999999598 , y2 = 0.499999999998822
x4 = 5.399999999999546 , y4 = 0.59999999999853
x5 = 5.499999999999902 , y5 = 0.59999999999853

For this representation the mapping on cell 0 will fail.
So, let's just consider the coordinate y1. If we approximate the value to 0.499999999998822 the mapping will be ok. However, if we try one of the following, as example:

y1 = 0.499999999998821
y1 = 0.4999999999988221

interesting enough, the mapping will fail again.
Another test: I add say two random digits to each vertex coordinate:

x1 = 5.39999999999924121 , y1 = 0.499999999998821821
x2 = 5.49999999999959821 , y2 = 0.49999999999882221
x4 = 5.39999999999954621 , y4 = 0.5999999999985321
x5 = 5.49999999999990221 , y5 = 0.5999999999985321

Now it works again...this is mind-blowing.
Just something wrong occurs when parsing the values.

Hope this is helpful.
Please let me know if I can provide some more useful information.

Thank you very much again,

Best,
Giovanni

[Jean-Paul Pelteret]

Dear Giovanni,

Hmm... this seems very strange. Thanks a lot for your test case. I'll try to look into it and see if there's some way that we can work around this problem. FYI. Here's the line where the vertex coordinates are read in, but its clear that they are stored as doubles (which matches the max precision of the coordinates in the test case). Curious indeed... I'll report back here if we come up with anything useful.

Best,
Jean-Paul

[Wolfgang Bangerth]

On 07/21/2017 02:39 AM, Giovanni Di Ilio wrote:
>
> So, let's just consider the coordinate y1. If we approximate the value to

> 0.49999999999882*2* the mapping will be ok. However, if we try one of the
> following, as example:
>
> y1 = 0.49999999999882*1*
> y1 = 0.49999999999882*21
>
> *interesting enough, the mapping will fail again.

> Another test: I add say two random digits to each vertex coordinate:
>

> x1 = 5.399999999999241*21* , y1 = 0.4999999999988218*21*
> x2 = 5.499999999999598*21* , y2 = 0.499999999998822*21*
> x4 = 5.399999999999546*21* , y4 = 0.59999999999853*21*
> x5 = 5.499999999999902*21* , y5 = 0.59999999999853*21*

>
> Now it works again...this is mind-blowing.

No, this is round-off error in action :-) Floating point numbers have only
approximately 16 digits of accuracy. The differences between the numbers you
have are at approximately at the level of round-off, and any computation you
do with these numbers will add to the effect of round-off.

If you have numbers where you depend on differences at the level of round-off,
arithmetic works very differently than for regular real numbers. They behave
more like integers because every operation involves rounding to the next
representable floating point number, and this rounding operation changes the
number significantly (relative to the differences between numbers you are
considering).

There is nothing you can do about this other than not rely on differences that
are this small relative to the absolute size of the numbers. One way you can
avoid this is if you have a small domain somewhere in R^2, move it to the
origin so that the relative distance between points is no longer small
compared to the distance from the origin.

[Giovanni Di Ilio]

Thanks Wolfgang and Jean-Paul for explanations.

Yes, I agree, the issue is definitely related to round-off error. But I still don't understand why the mapping fails, if the vertex coordinates read in from my .msh file are apparently stored as doubles...
Anyway, yes, the problem is easily overcome by "manipulating" the values before mapping.

Best,

Giovanni

[Wolfgang Bangerth]

On 07/22/2017 04:42 AM, Giovanni Di Ilio wrote:
>
> Yes, I agree, the issue is definitely related to round-off error. But I still
> don't understand why the mapping fails, if the vertex coordinates read in from
> my .msh file are apparently stored as doubles...
> Anyway, yes, the problem is easily overcome by "manipulating" the values
> before mapping.

The numbers you consider differ in the last two or three bits when stored as
double precision. Any computations you do with them, for example in the
mapping, will then necessarily be completely off.

[Giovanni Di Ilio]

Sorry, I thought I worked around the problem but actually I did not. The mapping fails again when I work with a refined mesh, for both the cases of a mesh read in from gmsh and a mesh generated with Deal.

I try to explain better what I mean:
1) I generate a grid, using: GridGenerator::subdivided_hyper_rectangle(triangulation,repetitions,p1,p2) ;
2) I refine one time the grid with: triangulation.refine_global(1) ;
3) I run over the cells and I try to get the mapped position of a given point, which I set to be in the centre of the considered cell: mapping.transform_real_to_unit_cell.

Since the point I am testing is exactly in the centre of the cell I expect the mapped position, in the unit cell, to be x = 0.5; y = 0.5.
However, for some "random" cell, again, the mapping fails.

Here I am not dealing anymore with the digits of accuracy.
The numbers I am dealing with are the coordinates of the vertexes provided by the GridGenerator and the coordinates of a testing point, which have just a couple of digits after comma. Moreover, the spacing between vertexes for the cells I am considering is of the order of 0.1.

I attach a very simple example, where I consider a uniform grid of 5 x 5 elements (repetitions[0]=repetitions[1]=5), refined one time, so 100 elements in total.
For cell with index 96 of my example, which vertexes are:
vertex(0): 0.8 0.8
vertex(1): 0.9 0.8
vertex(2): 0.8 0.9
vertex(3): 0.9 0.9
I test the point (0.85, 0.85) and the mapped position in the unit cell results to be: (0.5, 0.46875) instead of (0.5, 0.5).

You may notice that, by running the example with refineglobal(false), the mapping on an equivalent grid made of 10 x 10 elements, with no refinements (repetitions[0]=repetitions[1]=10), works just fine.
I don't see any difference between the two cases.
What is missing here? Can you help me please?

Thank you very much,
Best,

Giovanni

[Jean-Paul Pelteret]

Dear Giovanni,

Thank you very much for the clear and simple test case! I've played around with it and I don't see any fault in the way that you're trying to use the mapping functionality, so I think that its quite clear that there's something wrong with Mapping::transform_real_to_unit_cell. I've tested this with both deal.II 8.5 and the development version and the problem exists in both cases. I've opened a issue on GitHub that uses your example to demonstrate the problem. I'll report back here when we make some headway towards solving the issue!

Best regards,

Jean-Paul

[Giovanni Di Ilio]

Dear Jean-Paul,

(I tried also version 8.4.1 and same problem).
Ok, thanks a lot!

Best,
Giovanni

[Jean-Paul Pelteret]

Dear Giovanni,

This bug was confirmed and fixed, so hopefully if you try this now with the development version you should no longer have any issues. Thanks again for reporting the problem!

Best,

Jean-Paul

[Giovanni Di Ilio]

Great!
Thank you very much Jean-Paul for the very quick fix!

Best,
Giovanni

[Jean-Paul Pelteret]

David Wells fixed the issue, so all of the credit goes to him :-)