Comment 8 for bug 1167036

Hi,

Andrew, your gmsh file works just fine if you just remove the two lines with periodic input. Not quite sure why it messes things up, but like you said there is some directional stuff on the edges there. If you leave it to dolfin to compute the periodicity of facets and dofs, then it works.

Mikael

Den Apr 10, 2013 kl. 7:44 PM skrev Jan Blechta:

> On Wed, 10 Apr 2013 17:07:55 -0000
> Mikael Mortensen <email address hidden> wrote:
>> Hi,
>>
>> I get the same errors as Andrew report, but only for the mesh posted.
>> A similar mesh that I create myself from gmsh is ok. Could you please
>> post the 16.geo file so that I can run that as well?
>>
>> All master_slave entities seem to be correctly found. Seems to me like
>> an ordering or edge direction issue, a regular Lagrange element is ok
>> with all meshes. I thought the edge direction would be ok if the
>> ordering was correct, but mesh.order() does not help.
>>
>> By the way, I just had a quick look at dolfin-convert. Is there really
>> need for an XmLHandler/xml_writer anymore (with the mesh domains
>> incorporated in the mesh)? Wouldn't it be better to simply use dolfin
>> functions to write these mesh files? Like
>>
>> ifile = File('filename.xml.gz')
>> ifile << mesh
>
> Most (if not all) converters (every format has different implementation)
> don't use Mesh class. They simply manipulate entries between files by
> floats, numpy arrays, lists, ... As a results there is no point running
> it in parallel by the way.
>
> Jan

Thanks for clearing that out:-)

>
>>
>> Just curious
>>
>> Mikael
>>
>> Den Apr 10, 2013 kl. 8:53 AM skrev Andrew McRae:
>>
>>>> From eyeballing the .msh file, there are occasionally variations
>>>> in the
>>> last decimal place - O(10^-15) for an O(1) number.
>>>
>>> --
>>> You received this bug notification because you are a member of
>>> DOLFIN Team, which is subscribed to DOLFIN.
>>> https://bugs.launchpad.net/bugs/1167036
>>>
>>> Title:
>>> dolfin-convert/gmsh incompatible with periodic boundaries
>>>
>>> Status in DOLFIN:
>>> New
>>>
>>> Bug description:
>>> I think there might be some edge directionality problems. I can't
>>> explain it too easily, so I'll provide step-by-step instructions
>>> instead.
>>>
>>> Start with the attached gmsh files (unit square, characteristic
>>> length 1/16 and 1/32, periodic lines used where appropriate). I
>>> run 'dolfin-convert 16.msh 16.xml', 'dolfin-convert 32.msh
>>> 32.xml'. I then run the following python code to generate two
>>> more meshes from this: ======================== from dolfin import
>>> * mesh = Mesh('16.xml') File('16a.xml') << mesh File('16b.xml') <<
>>> refine(mesh) mesh = Mesh('32.xml')
>>> File('32a.xml') << mesh
>>> ========================
>>> The main difference between 16.xml and 16a.xml is: in 16a, the
>>> triangle definitions contain the vertices in ascending order by
>>> vertex number, whereas the order is 'random' in 16.xml (the same as
>>> in the .msh file). In addition, the useless third coordinate has
>>> been stripped out. 16b.xml is related to mesh 16a; each triangle
>>> has been subdivided into 4 (I think).
>>>
>>> I now run the following code, which sets up a periodic BDM function
>>> space on each of these four meshes (plus two others, to act as a
>>> reference), projects the simple field (y,0) into this space, and
>>> displays a plot and the errornorm.
>>>
>>> ==================================
>>> from dolfin import *
>>> mesh1 = Mesh('16.xml')
>>> mesh2 = Mesh('16a.xml')
>>> mesh3 = UnitSquareMesh(16,16)
>>> mesh4 = Mesh('32a.xml')
>>> mesh5 = Mesh('16b.xml')
>>> mesh6 = UnitSquareMesh(32,32)
>>> class PeriodicBoundaryX(SubDomain):
>>> def inside(self, x, on_boundary):
>>> return on_boundary and near(x[0], 0)
>>> def map(self, x, y):
>>> y[0] = x[0] - 1.0
>>> y[1] = x[1]
>>> pbc = PeriodicBoundaryX()
>>> S1 = FunctionSpace(mesh1,'BDM',1,constrained_domain=pbc)
>>> S2 = FunctionSpace(mesh2,'BDM',1,constrained_domain=pbc)
>>> S3 = FunctionSpace(mesh3,'BDM',1,constrained_domain=pbc)
>>> S4 = FunctionSpace(mesh4,'BDM',1,constrained_domain=pbc)
>>> S5 = FunctionSpace(mesh5,'BDM',1,constrained_domain=pbc)
>>> S6 = FunctionSpace(mesh6,'BDM',1,constrained_domain=pbc)
>>> expr = Expression(("x[1]","0.0"))
>>> u1 = project(expr,S1,solver_type="lu")
>>> u2 = project(expr,S2,solver_type="lu")
>>> u3 = project(expr,S3,solver_type="lu")
>>> u4 = project(expr,S4,solver_type="lu")
>>> u5 = project(expr,S5,solver_type="lu")
>>> u6 = project(expr,S6,solver_type="lu")
>>> print errornorm(expr,u1)
>>> print errornorm(expr,u2)
>>> print errornorm(expr,u3)
>>> print errornorm(expr,u4)
>>> print errornorm(expr,u5)
>>> print errornorm(expr,u6)
>>> plot(u1)
>>> plot(u2)
>>> plot(u3)
>>> plot(u4)
>>> plot(u5)
>>> plot(u6)
>>> interactive()
>>> ==================================
>>>
>>> We see that meshes 1, 2 and 4 all look broken - these are either the
>>> meshes output by dolfin-convert, or these meshes loaded in and
>>> written back out to file (sorting the vertices). In particular,
>>> the plots suggest that there are problems at the periodic sides of
>>> the domain. [This can be supported by changing the expression from
>>> (y,0) to (0,x), in which there is no flow across the periodic
>>> boundary]
>>>
>>> Meshes 3 and 6, the inbuilt UnitSquare meshes, are fine. Crucially,
>>> so is Mesh 5, which was generated by loading in the dolfin-converted
>>> mesh, refining once, and outputting the result to file.
>>>
>>> The errornorm quantifies what we see in the plots - nearly 0.1 for
>>> meshes 1, 2 and 4, but machine-precision for meshes 3, 5 and 6.
>>>
>>> I previously found the same problem and reported it as Bug #1082370,
>>> but I wasn't able to make such a strong case. The output here is
>>> far more damning!
>>>
>>> To manage notifications about this bug go to:
>>> https://bugs.launchpad.net/dolfin/+bug/1167036/+subscriptions
>>
>
> --
> You received this bug notification because you are a member of DOLFIN
> Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/1167036
>
> Title:
> dolfin-convert/gmsh incompatible with periodic boundaries
>
> Status in DOLFIN:
> New
>
> Bug description:
> I think there might be some edge directionality problems. I can't
> explain it too easily, so I'll provide step-by-step instructions
> instead.
>
> Start with the attached gmsh files (unit square, characteristic length 1/16 and 1/32, periodic lines used where appropriate). I run 'dolfin-convert 16.msh 16.xml', 'dolfin-convert 32.msh 32.xml'. I then run the following python code to generate two more meshes from this:
> ========================
> from dolfin import *
> mesh = Mesh('16.xml')
> File('16a.xml') << mesh
> File('16b.xml') << refine(mesh)
> mesh = Mesh('32.xml')
> File('32a.xml') << mesh
> ========================
> The main difference between 16.xml and 16a.xml is: in 16a, the triangle definitions contain the vertices in ascending order by vertex number, whereas the order is 'random' in 16.xml (the same as in the .msh file). In addition, the useless third coordinate has been stripped out. 16b.xml is related to mesh 16a; each triangle has been subdivided into 4 (I think).
>
> I now run the following code, which sets up a periodic BDM function
> space on each of these four meshes (plus two others, to act as a
> reference), projects the simple field (y,0) into this space, and
> displays a plot and the errornorm.
>
> ==================================
> from dolfin import *
> mesh1 = Mesh('16.xml')
> mesh2 = Mesh('16a.xml')
> mesh3 = UnitSquareMesh(16,16)
> mesh4 = Mesh('32a.xml')
> mesh5 = Mesh('16b.xml')
> mesh6 = UnitSquareMesh(32,32)
> class PeriodicBoundaryX(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and near(x[0], 0)
> def map(self, x, y):
> y[0] = x[0] - 1.0
> y[1] = x[1]
> pbc = PeriodicBoundaryX()
> S1 = FunctionSpace(mesh1,'BDM',1,constrained_domain=pbc)
> S2 = FunctionSpace(mesh2,'BDM',1,constrained_domain=pbc)
> S3 = FunctionSpace(mesh3,'BDM',1,constrained_domain=pbc)
> S4 = FunctionSpace(mesh4,'BDM',1,constrained_domain=pbc)
> S5 = FunctionSpace(mesh5,'BDM',1,constrained_domain=pbc)
> S6 = FunctionSpace(mesh6,'BDM',1,constrained_domain=pbc)
> expr = Expression(("x[1]","0.0"))
> u1 = project(expr,S1,solver_type="lu")
> u2 = project(expr,S2,solver_type="lu")
> u3 = project(expr,S3,solver_type="lu")
> u4 = project(expr,S4,solver_type="lu")
> u5 = project(expr,S5,solver_type="lu")
> u6 = project(expr,S6,solver_type="lu")
> print errornorm(expr,u1)
> print errornorm(expr,u2)
> print errornorm(expr,u3)
> print errornorm(expr,u4)
> print errornorm(expr,u5)
> print errornorm(expr,u6)
> plot(u1)
> plot(u2)
> plot(u3)
> plot(u4)
> plot(u5)
> plot(u6)
> interactive()
> ==================================
>
> We see that meshes 1, 2 and 4 all look broken - these are either the
> meshes output by dolfin-convert, or these meshes loaded in and written
> back out to file (sorting the vertices). In particular, the plots
> suggest that there are problems at the periodic sides of the domain.
> [This can be supported by changing the expression from (y,0) to (0,x),
> in which there is no flow across the periodic boundary]
>
> Meshes 3 and 6, the inbuilt UnitSquare meshes, are fine. Crucially,
> so is Mesh 5, which was generated by loading in the dolfin-converted
> mesh, refining once, and outputting the result to file.
>
> The errornorm quantifies what we see in the plots - nearly 0.1 for
> meshes 1, 2 and 4, but machine-precision for meshes 3, 5 and 6.
>
> I previously found the same problem and reported it as Bug #1082370,
> but I wasn't able to make such a strong case. The output here is far
> more damning!
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1167036/+subscriptions