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.


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.0)

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

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

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

mesh = UnitSquareMesh(64, 64)

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

boundaries = MeshFunction(“size_t”, mesh, 1)
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 = (inner(a0grad(u), grad(v))dx(0) + inner(a1grad(u), grad(v))dx(1)
- f
dx(0) - fvdx(1))

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

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


Mesh generation

import pygmsh
import meshio

pad_w = 0.5
pad_h = 1
geom = pygmsh.opencascade.Geometry(characteristic_length_min=0.02, characteristic_length_max=0.02)

pad_points = [[0, pad_h, 0], [pad_w, pad_h, 0], [pad_w, 0, 0], [0, 0, 0]]
pad_points1 = [[pad_w, pad_h, 0], [2pad_w, pad_h, 0], [pad_w2, 0, 0], [pad_w, 0, 0]]

cooling_pad_buf = geom.add_polygon(pad_points, lcar=4)
cooling_pad = geom.boolean_union([cooling_pad_buf, cooling_pad_buf])

cooling_pad_buf1 = geom.add_polygon(pad_points1, lcar=4)
cooling_pad1 = geom.boolean_union([cooling_pad_buf1, cooling_pad_buf1])

geom.add_physical(cooling_pad, label=“cooling_pad”)
geom.add_physical(cooling_pad1, label=“heat”)

cooling_pad_buf.line_loop.lines[3]], label=“sides_1”)
geom.add_physical(cooling_pad_buf1.line_loop.lines[1], label=“sides_1”)
geom.add_physical(cooling_pad_buf.line_loop.lines[3], label=“sides_2”)

cooling_pad_buf1.line_loop.lines[1]], label=“sides_4”)

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


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)

cooling_pad_buf1 = geom.add_polygon(pad_points1, lcar=4) #, make_surface=True)
cooling_pad1_surf = geom.add_plane_surface(cooling_pad_buf1.line_loop)

cooling_pad1 = geom.boolean_difference([cooling_pad1_surf], [cooling_pad], delete_first=True, delete_other=False)