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:

infile.read(mesh)

`

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

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 = (inner(a0*grad(u), grad(v)) dx(0) + inner(a1grad(u), grad(v))dx(1)*dx(0) - f

- fv

*v*dx(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

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], [2*pad_w, pad_h, 0], [pad_w*2, 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

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

`