dolfin-convert/gmsh incompatible with periodic boundaries

[Andrew McRae]

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!

[Garth Wells (garth-wells)]

I have had reports that gmsh periodic boundaries come with relatively large coordinate variations for points that should match. I'll add a user tolerance to the periodic bc computation at some point.

To test issues with periodic boundary conditions , the first step is to check master-slave pairs by visualising them using the function

      static MeshFunction<std::size_t> PeriodicBoundaryComputation:: masters_slaves(.......)

You can find it in dolfin/mesh/PeriodicBoundaryComputation.h. Send the computed MeshFunction for different dims (ie., vertices, edges, faces) to a VTK file. All entities will have index 0, 1 or 2, which corresponds to slave, master, nothing (I can't remember the precise order). You will easily spot any entities that do not match up.

Garth

[Andrew McRae (andymc)]

From eyeballing the .msh file, there are occasionally variations in the last decimal place - O(10^-15) for an O(1) number.

[Mikael Mortensen (mikael-mortensen)]

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

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 = Functi...

[Andrew McRae (andymc)]

Aha. That's interesting...
Aside: my gmsh version is 2.5.1 (the one that is installed by 'sudo apt-get install gmsh' on Precise)

====================
lc = 1/16;

Point(1) = {0, 0, 0, lc};
Point(2) = {1, 0, 0, lc};
Point(3) = {1, 1, 0, lc};
Point(4) = {0, 1, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Periodic Line {3} = {-1};
Periodic Line {4} = {-2};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(6) = {5};
====================

...is it my fault? :-P

[Jan Blechta (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

>
> 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 = Unit...

[Mikael Mortensen (mikael-mortensen)]

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 u...

[Andrew McRae (andymc)]

That's interesting, okay, thank you.

[Garth Wells (garth-wells)]

On 10 April 2013 14:53, Andrew McRae <email address hidden> wrote:
> >From eyeballing the .msh file, there are occasionally variations in the
> last decimal place - O(10^-15) for an O(1) number.
>

The code in PeriodicBoundaryComputation.cpp tests for matching
entities using DOLFIN_EPS, which is equal to 3.0e-16 (I don't know how
this number for DOLFIN_EPS came about).

Garth

> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/1167036
>
> Title:
> dolfin-convert/gmsh incompatible with periodic boundaries
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1167036/+subscriptions

[Andrew McRae (andymc)]

Mikael, can I ask what your exact procedure was? I got the same issue as before when I removed the periodic line.

[Mikael Mortensen (mikael-mortensen)]

Hi,

This was you geo file:

lc = 1/16;

Point(1) = {0, 0, 0, lc};
Point(2) = {1, 0, 0, lc};
Point(3) = {1, 1, 0, lc};
Point(4) = {0, 1, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Periodic Line {3} = {-1};
Periodic Line {4} = {-2};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(6) = {5};

I removed the two lines with Periodic such that the geo file was

lc = 1/16;

Point(1) = {0, 0, 0, lc};
Point(2) = {1, 0, 0, lc};
Point(3) = {1, 1, 0, lc};
Point(4) = {0, 1, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(6) = {5};

I opened this i gmsh, computed a mesh, used dolfin-convert on it and then
it worked. I retried it now just to make sure.

Mikael

On 11 April 2013 12:03, Andrew McRae <email address hidden> wrote:

> Mikael, can I ask what your exact procedure was? I got the same issue
> as before when I removed the periodic line.
>
> --
> 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_domai...

[Garth Wells (garth-wells)]

Could you test this again against the dev version on Bitbucket? I've pushed some changes that allow an optional tolerance argument, and I've reduced the default tolerance to 1e-10. To see how to change the tolerance, look at

    test/unit/fem/python/PeriodicBC.py

[Garth Wells (garth-wells)]

Note also that rather than

    on_boundary and near(x[0], 0)

you should use

    on_boundary and near(x[0], 0.0, 1.0e-10)

the function near uses a very strict tolerance.

[Andrew McRae (andymc)]

Mikael: I think that the periodicity isn't being (fully?) enforced when the "periodic" in the .geo file are removed, maybe because corresponding points are too far apart. The test I supplied then appears to work, but for the wrong reason. If I use the meshes for a dynamic simulation (shallow water equations), I get nonsense (blow-up from the edges) when I attempt to run a dynamic simulation with it. Indeed, I only get sensible results on my 16b mesh if the periodic lines are as they are written (e.g. changing -1 and -2 to 1 and 2 doesn't fix it!)

Garth: will try...

[Andrew McRae (andymc)]

Garth: no qualitative changes, "except" that (see my comment to Mikael) removing the Periodic specification from the .geo file now gives non-machine-precision errors, because the 'periodic DOF finder' now picks up the points that were previously too far apart!

Thanks for the advice regarding 'near'. near(x[0], 0) should be fine for my purposes, but near(x[0], 1e9) certainly wouldn't have been, and might explain some other issues I'd forgotten about.