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!
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<
You can find it in dolfin/
Garth