 # Meshio xdmf mesh doesn't give same results as with unitsquare-mesh

I was switching over to create my meshes in xdmf format (pygmsh -> meshio).
However, using meshio xdmf meshes doesn’t give me correct looking results. For replicating the problem, I took the two-materials demo (two materials with different conduction placed next to each other; DC boundary condition on left = 1, DC boundary condition on right = 0). When replacing the unitsquare mesh with the xdmf mesh

`
mesh = UnitSquareMesh(64, 64)

replaced by:

mesh = Mesh()
with XDMFFile(“unitsquare.xdmf”) as infile:
`
. The result changes from (a) to (b).

Is there something I have to change how I import the xdmf-mesh, or is there something off in the mesh generation itself already?

If someone could shed some light on this, that would be great.
Steffen

`

# Fenics code

from dolfin import *
import meshio
import numpy as np
import mshr

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x, 0.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x, 1.0)

class Material2(SubDomain):
def inside(self, x, on_boundary):
return between(x, (0.5, 1.0))

left = Left()
right = Right()
material = Material2()

mesh = UnitSquareMesh(64, 64)

domains = MeshFunction(“size_t”, mesh, 2)
domains.set_all(0)
material.mark(domains, 1)

boundaries = MeshFunction(“size_t”, mesh, 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 3)

a0 = Constant(1.0)
a1 = Constant(10.0)
f = Constant(0.0)

V = FunctionSpace(mesh, “CG”, 2)
u = TrialFunction(V)
v = TestFunction(V)

bcs = [DirichletBC(V, 1.0, boundaries, 1),
DirichletBC(V, 0.0, boundaries, 3)]

dx = Measure(‘dx’, domain=mesh, subdomain_data=domains)
ds = Measure(‘ds’, domain=mesh, subdomain_data=boundaries)

- f
v
dx(0) - fvdx(1))

a, L = lhs(F), rhs(F)

u = Function(V)
solve(a == L, u, bcs)

XDMFFile(‘two_material_no_external_meshing.xdmf’).write(u)

# Mesh generation

import pygmsh
import meshio

mesh = pygmsh.generate_mesh(geom, extra_gmsh_arguments=["-format", “msh2”])

# Meshio

meshio.write(‘unitsquare.xdmf’, meshio.Mesh(points=mesh.points, cells={“triangle”: mesh.cells[“triangle”]}))
`

Found it, in pygmsh I was basically creating two rectangles that were not sharing the middle line, but had a boundary between them instead. In Fenics that boiled down to having an implicit Neumann boundary condition at the intersection.

the pygmsh file got fixed as such

``````cooling_pad = geom.add_plane_surface(cooling_pad_buf.line_loop)