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(a0grad(u), grad(v))dx(0) + inner(a1grad(u), grad(v))dx(1)
- fvdx(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
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
meshio.write(‘unitsquare.xdmf’, meshio.Mesh(points=mesh.points, cells={“triangle”: mesh.cells[“triangle”]}))
`