Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Terribly hacky MVC workaround for pygmsh-meshio-dolfin boundary tagging
from pygmsh.built_in.geometry import Geometry
from pygmsh import generate_mesh
import meshio
import dolfin
geom = Geometry()
lc = .1
p0 = geom.add_point([0, 0, 0], lcar=lc)
p1 = geom.add_point([1, 0, 0], lcar=lc)
p2 = geom.add_point([1, 1, 0], lcar=lc)
p3 = geom.add_point([0, 1, 0], lcar=lc)
l0 = geom.add_line(p0, p1)
l1 = geom.add_line(p1, p2)
l2 = geom.add_line(p2, p3)
l3 = geom.add_line(p3, p0)
ll = geom.add_line_loop(lines=[l0, l1, l2, l3])
ps = geom.add_plane_surface(ll)
# Tag line and surface
geom.add_physical_line(lines=l3, label=444)
geom.add_physical_surface(surfaces=ps, label=456)
print("\n".join(geom._GMSH_CODE))
points, cells, point_data, cell_data, field_data = generate_mesh(geom)
print("Writing 2d mesh for dolfin Mesh")
meshio.write("mesh_2d.xdmf", meshio.Mesh(
points=points,
cells={"triangle": cells["triangle"]}))
print("Writing 1d line data for dolfin MVC")
meshio.write("mvc_1d.xdmf", meshio.Mesh(
points=points,
cells={"line": cells["line"]},
cell_data={"line": {"name_to_read": cell_data["line"]["gmsh:physical"]}}
))
mesh_2d = dolfin.Mesh()
with dolfin.XDMFFile("mesh_2d.xdmf") as infile:
print("Reading 2d mesh into dolfin")
infile.read(mesh_2d)
mvc = dolfin.MeshValueCollection("size_t", mesh_2d, 1)
with dolfin.XDMFFile("mvc_1d.xdmf") as infile:
print("Reading 1d line data into dolfin mvc")
infile.read(mvc, "name_to_read")
print("Constructing MeshFunction from MeshValueCollection")
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh_2d, mvc)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.