Comment 7 for bug 1167036

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