Comment 12 for bug 1167036

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_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
>