Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-27 Thread Jean-Paul Pelteret 
David Wells fixed the issue, so all of the credit goes to him :-)

On Thursday, July 27, 2017 at 10:12:16 AM UTC+2, Giovanni Di Ilio wrote:
>
> Great!
> Thank you very much Jean-Paul for the very quick fix!
>
> Best,
> Giovanni
>
>
>
> On Thursday, July 27, 2017 at 9:47:48 AM UTC+2, Jean-Paul Pelteret wrote:
>>
>> 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
>>
>> On Tuesday, July 25, 2017 at 10:13:51 AM UTC+2, Giovanni Di Ilio wrote:
>>>
>>> Dear Jean-Paul,
>>>
>>> (I tried also version 8.4.1 and same problem).
>>> Ok, thanks a lot!
>>>
>>> Best,
>>> Giovanni
>>>
>>> On Tuesday, July 25, 2017 at 9:08:53 AM UTC+2, Jean-Paul Pelteret wrote:

 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

 On Tuesday, July 25, 2017 at 12:29:51 AM UTC+2, Giovanni Di Ilio wrote:
>
> 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
>
>
>
> On Saturday, July 22, 2017 at 8:57:57 PM UTC+2, Wolfgang Bangerth 
> wrote:
>>
>> 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. 
>>
>> Best 
>>   W. 
>>
>> -- 
>>  
>>
>> Wolfgang Bangerth  email: 
>> bang...@colostate.edu 
>> www: 
>> http://www.math.colostate.edu/~bangerth/ 
>>

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-27 Thread Giovanni Di Ilio 
Great!
Thank you very much Jean-Paul for the very quick fix!

Best,
Giovanni



On Thursday, July 27, 2017 at 9:47:48 AM UTC+2, Jean-Paul Pelteret wrote:
>
> 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
>
> On Tuesday, July 25, 2017 at 10:13:51 AM UTC+2, Giovanni Di Ilio wrote:
>>
>> Dear Jean-Paul,
>>
>> (I tried also version 8.4.1 and same problem).
>> Ok, thanks a lot!
>>
>> Best,
>> Giovanni
>>
>> On Tuesday, July 25, 2017 at 9:08:53 AM UTC+2, Jean-Paul Pelteret wrote:
>>>
>>> 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
>>>
>>> On Tuesday, July 25, 2017 at 12:29:51 AM UTC+2, Giovanni Di Ilio wrote:

 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



 On Saturday, July 22, 2017 at 8:57:57 PM UTC+2, Wolfgang Bangerth wrote:
>
> 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. 
>
> Best 
>   W. 
>
> -- 
>  
>
> Wolfgang Bangerth  email: 
> bang...@colostate.edu 
> www: 
> http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-25 Thread Giovanni Di Ilio 
Dear Jean-Paul,

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

Best,
Giovanni

On Tuesday, July 25, 2017 at 9:08:53 AM UTC+2, Jean-Paul Pelteret wrote:
>
> 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
>
> On Tuesday, July 25, 2017 at 12:29:51 AM UTC+2, Giovanni Di Ilio wrote:
>>
>> 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
>>
>>
>>
>> On Saturday, July 22, 2017 at 8:57:57 PM UTC+2, Wolfgang Bangerth wrote:
>>>
>>> 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. 
>>>
>>> Best 
>>>   W. 
>>>
>>> -- 
>>>  
>>> Wolfgang Bangerth  email: bang...@colostate.edu 
>>> www: 
>>> http://www.math.colostate.edu/~bangerth/ 
>>>
>>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-25 Thread 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

On Tuesday, July 25, 2017 at 12:29:51 AM UTC+2, Giovanni Di Ilio wrote:
>
> 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
>
>
>
> On Saturday, July 22, 2017 at 8:57:57 PM UTC+2, Wolfgang Bangerth wrote:
>>
>> 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. 
>>
>> Best 
>>   W. 
>>
>> -- 
>>  
>> Wolfgang Bangerth  email: bang...@colostate.edu 
>> www: http://www.math.colostate.edu/~bangerth/ 
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-24 Thread 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



On Saturday, July 22, 2017 at 8:57:57 PM UTC+2, Wolfgang Bangerth wrote:
>
> 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. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

using namespace dealii;

class RefineTest
{

public:
static constexpr int dim  = 2 ;

public:

RefineTest(bool _refine)
:
triangulation(),
dof_handler(triangulation),
refine(_refine)
{
make_grid(); 
test_mapping();
}

private:
void make_grid()
{
if(refine)
{
std::vector< unsigned  int > repetitions(dim);
repetitions[0] = repetitions[1] = 5;
std::array L = {1.0, 1.0};
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point(0.0,0.0), Point(L[0],L[1]), true);
 	triangulation.refine_global(1);
}
else
{
std::vector< unsigned  int > repetitions(dim);
repetitions[0] = repetitions[1] = 10;
std::array L = {1.0, 1.0};
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point(0.0,0.0), Point(L[0],L[1]), true);
 	triangulation.refine_global(0); 
}
}

void test_mapping()
{
auto cell = dof_handler.begin_active();
auto endc = dof_handler.end();
for (; cell!=endc; ++cell)
{

auto cellN = cell->index();
auto P1 = cell->vertex(0);
	auto P2 = cell->vertex(1);
	auto P3 = cell->vertex(2);

	Point test;
test(0) = (P1[0]+P2[0])/2.0; 
test(1) = (P2[1]+P3[1])/2.0; 

MappingQ mapping(1);

Point dp = mapping.transform_real_to_unit_cell(cell,test);
Point dp_real = mapping.transform_unit_to_real_cell(cell,dp);

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-22 Thread 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



On Friday, July 21, 2017 at 5:29:29 PM UTC+2, Wolfgang Bangerth wrote:
>
> 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.499882*2* the mapping will be ok. However, if we try one of 
> the 
> > following, as example: 
> > 
> > y1 = 0.499882*1* 
> > y1 = 0.499882*21 
> > 
> > *interesting enough, the mapping will fail again. 
> > Another test: I add say two random digits to each vertex coordinate: 
> > 
> > x1 = 5.3999241*21* , y1 = 0.49988218*21* 
> > x2 = 5.4999598*21* , y2 = 0.4998822*21* 
> > x4 = 5.3999546*21* , y4 = 0.599853*21* 
> > x5 = 5.402*21* , y5 = 0.599853*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. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-21 Thread 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.499882*2* the mapping will be ok. However, if we try one of the 
following, as example:


y1 = 0.499882*1*
y1 = 0.499882*21

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

x1 = 5.3999241*21* , y1 = 0.49988218*21*
x2 = 5.4999598*21* , y2 = 0.4998822*21*
x4 = 5.3999546*21* , y4 = 0.599853*21*
x5 = 5.402*21* , y5 = 0.599853*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.


Best
 W.

--

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

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-21 Thread 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

On Friday, July 21, 2017 at 10:39:55 AM UTC+2, Giovanni Di Ilio wrote:
>
> 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.3999241 , y1 = 0.49988218
> x2 = 5.4999598 , y2 = 0.4998822
> x4 = 5.3999546 , y4 = 0.599853
> x5 = 5.402 , y5 = 0.599853
>
> 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.499882*2* the mapping will be ok. However, if we try one of the 
> following, as example:
>
> y1 = 0.499882*1*
> y1 = 0.499882
>
> *21*interesting enough, the mapping will fail again.
> Another test: I add say two random digits to each vertex coordinate:
>
> x1 = 5.3999241*21* , y1 = 0.49988218*21*
> x2 = 5.4999598*21* , y2 = 0.4998822*21*
> x4 = 5.3999546*21* , y4 = 0.599853*21*
> x5 = 5.402*21* , y5 = 0.599853*21*
>
> 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
>
>
>
> On Thursday, July 20, 2017 at 9:08:34 PM UTC+2, Jean-Paul Pelteret wrote:
>>
>> 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.19991087 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
>>
>> On Thursday, July 20, 2017 at 7:47:04 PM UTC+2, Giovanni Di Ilio wrote:
>>>
>>> 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
>>>
>>> On Thursday, July 20, 2017 at 4:32:43 PM UTC+2, Wolfgang Bangerth wrote:

 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/ 



-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-21 Thread 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.3999241 , y1 = 0.49988218
x2 = 5.4999598 , y2 = 0.4998822
x4 = 5.3999546 , y4 = 0.599853
x5 = 5.402 , y5 = 0.599853

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.499882*2* the mapping will be ok. However, if we try one of the 
following, as example:

y1 = 0.499882*1*
y1 = 0.499882

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

x1 = 5.3999241*21* , y1 = 0.49988218*21*
x2 = 5.4999598*21* , y2 = 0.4998822*21*
x4 = 5.3999546*21* , y4 = 0.599853*21*
x5 = 5.402*21* , y5 = 0.599853*21*

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



On Thursday, July 20, 2017 at 9:08:34 PM UTC+2, Jean-Paul Pelteret wrote:
>
> 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.19991087 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
>
> On Thursday, July 20, 2017 at 7:47:04 PM UTC+2, Giovanni Di Ilio wrote:
>>
>> 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
>>
>> On Thursday, July 20, 2017 at 4:32:43 PM UTC+2, Wolfgang Bangerth wrote:
>>>
>>> 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/ 
>>>
>>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

using namespace dealii;

class GmeshTest
{

public:
static constexpr int dim  = 2 ;

public:

GmeshTest(bool _use_hypercube)
:
triangulation(),
dof_handler(triangulation),
	use_hypercube(_use_hypercube)
{
make_grid(); 
test_mapping();
}

private:
void make_grid()
{

if(use_hypercube)
{
std::vector< unsigned  int > repetitions(dim);
repetitions[0] = 2;
repetitions[1] = 1;
GridGenerator::subdivided_hyper_rectangle(triangulation, repetitions, Point(5.4,0.5),

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-20 Thread 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.19991087 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

On Thursday, July 20, 2017 at 7:47:04 PM UTC+2, Giovanni Di Ilio wrote:
>
> 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
>
> On Thursday, July 20, 2017 at 4:32:43 PM UTC+2, Wolfgang Bangerth wrote:
>>
>> 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/ 
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-20 Thread 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

On Thursday, July 20, 2017 at 4:32:43 PM UTC+2, Wolfgang Bangerth wrote:
>
> 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/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Mapping points of mesh generated by gmsh

2017-07-20 Thread 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: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Mapping points of mesh generated by gmsh

2017-07-20 Thread 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


-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

using namespace dealii;

class GmeshTest
{

public:
static constexpr int dim  = 2 ;

public:

GmeshTest(bool _use_hypercube)
:
triangulation(),
dof_handler(triangulation),
use_hypercube(_use_hypercube)
{
make_grid(); 
test_mapping();
}

private:
void make_grid()
{

std::string outmesh_name;
if(use_hypercube)
{
std::vector< unsigned  int > repetitions(dim);
repetitions[0] = 60;
repetitions[1] = 20;

std::array L = {6.0, 2.0};
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point(0.0,0.0), Point(L[0],L[1]), true);
outmesh_name="mesh_hypercube.vtk";
}
else
{
GridIn gridin;
gridin.attach_triangulation(triangulation);
std::ifstream f("test.msh"); 
gridin.read_msh(f);
outmesh_name="mesh_gmsh.vtk";
}

std::ofstream out(outmesh_name);
GridOut grid_out;
grid_out.write_vtk(triangulation, out);

}

void test_mapping()
{
auto cell = dof_handler.begin_active();
auto endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
auto cellN = cell->index();
if (cellN == 354)
{
Point test;
test(0) = 5.45;
test(1) = 0.55;

MappingQ mapping(1);

Point dp = mapping.transform_real_to_unit_cell(cell,test);
Point dp_real = mapping. transform_unit_to_real_cell(cell,dp);

std::cout << "cell=" << cell->index() << std::endl;
std::cout << "cell vertex(0): " << cell->vertex(0) << std::endl;
std::cout << "cell vertex(1): " << cell->vertex(1) << std::endl;
std::cout << "cell vertex(2): " << cell->vertex(2) << std::endl;
std::cout << "cell vertex(3): " << cell->vertex(3) << std::endl;
std::cout << "=" << std::endl;
std::cout << "test_point=   " << test << "  mapped_point=   " << dp <<" backmapping: "<
    
    
13 matches