Gmsh 4.4.1 in FEniCS? Meshio

mesh
1

Hi all, I understand there is a previous thread on a similar issue: https://fenicsproject.discourse.group/t/transitioning-from-mesh-xml-to-mesh-xdmf-from-dolfin-convert-to-meshio/412/18 but I’m not sure if this problem has to do with the version of gmsh.
I have 2 examples of mesh that I wish to work with: 1) Box.msh (2D) 2) Cylinder (3D)

For the 2D problem I write:

msh = meshio.read(“Box.msh”)

meshio.write(“Box.xdmf”, meshio.Mesh(points=msh.points, cells={“triangle”: msh.cells[“triangle”]}))

and for the 3D problem I write:

msh = meshio.read(“Cylinder.msh”)

meshio.write(“Cylinder.xdmf”, meshio.Mesh(points=msh.points, cells={“triangle”: msh.cells[“triangle”]}))

In both cases, I get a segmentation fault error:

[cbl-02:08311] *** Process received signal ***
[cbl-02:08311] Signal: Segmentation fault (11)
[cbl-02:08311] Signal code: Address not mapped (1)
[cbl-02:08311] Failing at address: 0x7ff500000720
[cbl-02:08311] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x3ef20)[0x7ff54f2a0f20]
[cbl-02:08311] [ 1] /lib/x86_64-linux-gnu/libc.so.6(+0x18e5c7)[0x7ff54f3f05c7]
[cbl-02:08311] [ 2] /home/karthic/.local/lib/python3.6/site-packages/h5py/_errors.cpython-36m-x86_64-linux-gnu.so(+0x43b7)[0x7ff50bdf73b7]
[cbl-02:08311] [ 3] /home/karthic/.local/lib/python3.6/site-packages/h5py/defs.cpython-36m-x86_64-linux-gnu.so(+0xc23e)[0x7ff50a82923e]
[cbl-02:08311] [ 4] /home/karthic/.local/lib/python3.6/site-packages/h5py/h5t.cpython-36m-x86_64-linux-gnu.so(+0x11b52)[0x7ff50a5bbb52]
[cbl-02:08311] [ 5] /home/karthic/.local/lib/python3.6/site-packages/h5py/h5t.cpython-36m-x86_64-linux-gnu.so(+0x9aef)[0x7ff50a5b3aef]
[cbl-02:08311] [ 6] /home/karthic/.local/lib/python3.6/site-packages/h5py/h5t.cpython-36m-x86_64-linux-gnu.so(PyInit_h5t+0x3a77)[0x7ff50a5f2a87]
[cbl-02:08311] [ 7] python3(_PyImport_LoadDynamicModuleWithSpec+0x178)[0x5e4268]
[cbl-02:08311] [ 8] python3[0x5e4522]
[cbl-02:08311] [ 9] python3(PyCFunction_Call+0x13e)[0x56246e]
[cbl-02:08311] [10] python3(_PyEval_EvalFrameDefault+0x58c6)[0x4fed26]
[cbl-02:08311] [11] python3[0x4f6128]
[cbl-02:08311] [12] python3[0x4f7d60]
[cbl-02:08311] [13] python3[0x4f876d]
[cbl-02:08311] [14] python3(_PyEval_EvalFrameDefault+0x467)[0x4f98c7]
[cbl-02:08311] [15] python3[0x4f7a28]
[cbl-02:08311] [16] python3[0x4f876d]
[cbl-02:08311] [17] python3(_PyEval_EvalFrameDefault+0x467)[0x4f98c7]
[cbl-02:08311] [18] python3[0x4f7a28]
[cbl-02:08311] [19] python3[0x4f876d]
[cbl-02:08311] [20] python3(_PyEval_EvalFrameDefault+0x467)[0x4f98c7]
[cbl-02:08311] [21] python3[0x4f7a28]
[cbl-02:08311] [22] python3[0x4f876d]
[cbl-02:08311] [23] python3(_PyEval_EvalFrameDefault+0x467)[0x4f98c7]
[cbl-02:08311] [24] python3[0x4f7a28]
[cbl-02:08311] [25] python3[0x4f876d]
[cbl-02:08311] [26] python3(_PyEval_EvalFrameDefault+0x467)[0x4f98c7]
[cbl-02:08311] [27] python3(_PyFunction_FastCallDict+0xf5)[0x4f4065]
[cbl-02:08311] [28] python3(_PyObject_FastCallDict+0x4f1)[0x57c8f1]
[cbl-02:08311] [29] python3(_PyObject_CallMethodIdObjArgs+0xee)[0x57cc5e]
[cbl-02:08311] *** End of error message ***
Segmentation fault (core dumped)

Any advice would be greatly appreciated. An example of the Box.msh is also given. Also, is there an example that describes how to read in meshes from gmsh? Thanks

Box.msh:

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
4 4 1 0
1 1 0 0 1 1
2 1 1 0 1 1
3 0 1 0 1 1
4 0 0 0 1 1
1 0 0 0 0 1 0 1 2 2 3 -4
2 0 0 0 1 0 0 1 2 2 4 -1
3 1 0 0 1 1 0 1 2 2 1 -2
4 0 1 0 1 1 0 1 2 2 2 -3
1 0 0 0 1 1 0 1 3 4 1 2 3 4
$EndEntities
$Nodes
9 13 1 13
0 1 0 1
1
1 0 0
0 2 0 1
2
1 1 0
0 3 0 1
3
0 1 0
0 4 0 1
4
0 0 0
1 1 0 1
5
0 0.5000000000020591 0
1 2 0 1
6
0.499999999998694 0 0
1 3 0 1
7
1 0.499999999998694 0
1 4 0 1
8
0.5000000000020591 1 0
2 1 0 5
9
10
11
12
13
0.5000000000001412 0.5000000000001882 0
0.2500000000010295 0.7500000000010295 0
0.7500000000005501 0.7499999999997204 0
0.7499999999997088 0.2499999999997206 0
0.2499999999997088 0.2500000000005618 0
$EndNodes
$Elements
9 28 1 28
0 1 15 1
1 1
0 2 15 1
2 2
0 3 15 1
3 3
0 4 15 1
4 4
1 1 1 2
5 3 5
6 5 4
1 2 1 2
7 4 6
8 6 1
1 3 1 2
9 1 7
10 7 2
1 4 1 2
11 2 8
12 8 3
2 1 2 16
13 8 3 10
14 5 4 13
15 7 2 11
16 6 1 12
17 5 9 10
18 9 8 10
19 8 9 11
20 9 5 13
21 9 7 11
22 6 9 13
23 9 6 12
24 7 9 12
25 4 6 13
26 1 7 12
27 3 5 10
28 2 8 11
$EndElements

2

As far as i know, your .msh file has to be in GMSH Version 2 format, but yours i in Version 4 (first line in yout .msh file)
Edit: Also, for your 3D file, you have look at the type of finite element your gmsh file generated (look at this), i assume you used tetrahedrons, so your key word should be

meshio.write(“Cylinder.xdmf”, meshio.Mesh(points=msh.points, cells={“tetra”: msh.cells[“tetra”]}))
3

Hi,

Thanks for getting back to me. I’ve tried it using gmsh version 2 and 3 as well but the same problem persists. I have HDF5 version 1.10.0 installed on linux 18.04. This clashes with h5py version 2.9 (which requires HDF5 version 1.10.4) so I uninstalled and reinstalled h5py version 2.7. I’m not sure if that might be the problem? If it is I’m not sure how to resolve it either.

4

Ok that sounds really weird to me, unfortuantely, i don’t think i can help you with this type of problem.
How did you installed the meshio package? For me, when i installed it using

pip3 install meshio[all] --user

all dependencies necessary were installed correctly, so there is no conflict with HDF5 and h5py.

5

Hi,

Reinstalled again since I did not add [all] previously. Also, reinstalled h5py to version 2.9 but I got the same error:

local/lib/python3.6/site-packages/h5py/init.py:75: UserWarning: h5py is running against HDF5 1.10.0 when it was built against 1.10.4, this may cause problems
‘{0}.{1}.{2}’.format(*version.hdf5_built_version_tuple)
Traceback (most recent call last):
File “/home/karthic/Documents/Fenics Scripts/DiffusionTest.py”, line 12, in
meshio.write(“Box.xdmf”, meshio.Mesh(points=msh.points, cells={“triangle”: msh.cells[“triangle”]}))
File “/home/karthic/.local/lib/python3.6/site-packages/meshio/_helpers.py”, line 240, in write
return interface.write(filename, mesh, *args, **_kwargs)
File “/home/karthic/.local/lib/python3.6/site-packages/meshio/_xdmf/main.py”, line 465, in write
XdmfWriter(*args, **kwargs)
File “/home/karthic/.local/lib/python3.6/site-packages/meshio/_xdmf/main.py”, line 275, in init
self.h5_file = h5py.File(self.h5_filename, “w”)
File “/home/karthic/.local/lib/python3.6/site-packages/h5py/_hl/files.py”, line 391, in init
fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0, **kwds)
File “/home/karthic/.local/lib/python3.6/site-packages/h5py/_hl/files.py”, line 108, in make_fapl
plist.set_libver_bounds(low, high)
File “h5py/_objects.pyx”, line 54, in h5py._objects.with_phil.wrapper
File “h5py/_objects.pyx”, line 55, in h5py._objects.with_phil.wrapper
File “h5py/h5p.pyx”, line 1140, in h5py.h5p.PropFAID.set_libver_bounds
ValueError: Invalid high library version bound (invalid high library version bound)

6

Ok this seems to be very strange to me. Are you using anaconda?
If so, maybe creating a new environment may solve the problem, otherwise i can’t help you any further, sorry!

7

Hi,

I’m using vscode not anaconda. Cannot seem to find any similar problem either. Thanks anyways, appreciate the help :slight_smile:

8

Ok then i’m really not able to help you any further, i’m sorry. For me, anaconda works perfectly fine.

9

Hi Kart,

I don’t identify the error message. But I suggest you to uninstall meshio and h5py again and try the following commands:

> pip3 install meshio --user
#Install h5py using your local version of HDF5 with:
> pip3 install --no-binary=h5py h5py --user

Also, about your request, you can replicate the following example and apply it to your specific case so you may read GMSH meshes correctly:

GMSH 4.4.1 code (2 Dimensional mesh):

Point(6) = {0, 0, 0, 1};
Point(7) = {2, 0, 0, 1};
Point(8) = {0, 2, 0, 1};
Point(9) = {2, 2, 0, 1};
Line(10) = {6, 7};               
Line(20) = {7, 9};           
Line(30) = {9, 8};       
Line(40) = {8, 6};       
Line Loop(1) = {10, 20, 30, 40};
//+
Plane Surface(1) = {1};
Physical Surface(0) = {1};
//+
Physical Curve(1) = {40};
Physical Curve(2) = {20};
Physical Curve(3) = {30};
Physical Curve(4) = {10};

Simple Poisson test problem:

import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp, Measure,\
                   DirichletBC, FunctionSpace, Constant, TrialFunction, \
                   TestFunction, dot, grad, dx, Function, solve

msh = meshio.read("./Square.msh")
Wri_path = './Dolfin_mesh_functions/'
meshio.write(Wri_path+"mesh.xdmf",
             meshio.Mesh(points=msh.points,
                         cells={"triangle": msh.cells["triangle"]}))

meshio.write(Wri_path+"mf.xdmf",
             meshio.Mesh(points=msh.points,
                         cells={"line": msh.cells["line"]},
                         cell_data={"line": {"name_to_read": msh.cell_data["line"]["gmsh:physical"]}}))

mesh = Mesh()
with XDMFFile(Wri_path+"mesh.xdmf") as infile:
    infile.read(mesh)
File(Wri_path+"Dolfin_circle_mesh.pvd").write(mesh)

mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile(Wri_path+"mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File(Wri_path+"Dolfin_circle_facets.pvd").write(mf)

V = FunctionSpace(mesh, 'P', 1)

# Define boundary conditions base on GMSH mesh marks [Physical Curves: 1, 2, 3, 4]
bc1 = DirichletBC(V, Constant(0.0), mf, 1)
bc2 = DirichletBC(V, Constant(10.0), mf, 2)
bc3 = DirichletBC(V, Constant(6.0), mf, 3)
bc4 = DirichletBC(V, Constant(3.0), mf, 4)

bc = [bc1, bc2, bc3, bc4]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-15.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)
File("Solution.pvd").write(u)

I hope you find it useful.
Best regards,

Santiago

4 Likes
10

Hi

Thank you for your example. It was very useful. I just wanted to point out that I was getting the error:

KeyError: ‘line’

I had to made the following change in the gmsh geo file

Physical Curve -> Physical Line

Also, just to make all the steps explicit, I created the msh file with

gmsh -2 Square.geo -p

I think that the flag -p is not strictly necessary

Best regards,
Pedro A. Vazquez

11

Sure!
It depends on the GMSH version one is using, I have the latest 4.4.1.

All the best,
Santiago

12

Hi Santiago,

Very sorry for the late reply, I have been away for some time. I have not had the chance to run the sample test problem you provided but your solution worked like a charm! Thank you so much! Appreciate the help

Best Regards,
Kart

13

Hi,

This may be a bit late, but I had some clashes like yours that got fixed by importing h5py before importing meshio (don’t really remember if the order matters) even if I didn’t directly need h5py in my program.
Looks like you already found a solution but I’m posting this for anyone with a similar issue who may run into this thread

Cheers,
Miguel

1 Like
14

I would like to point out that this answer doesn’t work for me. After uninstalling meshio and h5py, reinstalling them using the above comands yield a gigantic error message for the h5py package.

Here’s the last part of that gigantic error message:

h5py/h5a.pyx:184:36: Cannot convert 'void *' to Python object
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/tmp/pip-install-m_w5y09w/h5py/setup.py", line 159, in <module>
    cmdclass = CMDCLASS,
  File "/usr/local/lib/python3.6/dist-packages/setuptools/__init__.py", line 145, in setup
    return distutils.core.setup(**attrs)
  File "/usr/lib/python3.6/distutils/core.py", line 148, in setup
    dist.run_commands()
  File "/usr/lib/python3.6/distutils/dist.py", line 955, in run_commands
    self.run_command(cmd)
  File "/usr/lib/python3.6/distutils/dist.py", line 974, in run_command
    cmd_obj.run()
  File "/usr/local/lib/python3.6/dist-packages/setuptools/command/install.py", line 61, in run
    return orig.install.run(self)
  File "/usr/lib/python3.6/distutils/command/install.py", line 589, in run
    self.run_command('build')
  File "/usr/lib/python3.6/distutils/cmd.py", line 313, in run_command
    self.distribution.run_command(command)
  File "/usr/lib/python3.6/distutils/dist.py", line 974, in run_command
    cmd_obj.run()
  File "/usr/lib/python3.6/distutils/command/build.py", line 135, in run
    self.run_command(cmd_name)
  File "/usr/lib/python3.6/distutils/cmd.py", line 313, in run_command
    self.distribution.run_command(command)
  File "/usr/lib/python3.6/distutils/dist.py", line 974, in run_command
    cmd_obj.run()
  File "/tmp/pip-install-m_w5y09w/h5py/setup_build.py", line 209, in run
    language_level=2)
  File "/tmp/pip-install-m_w5y09w/h5py/.eggs/Cython-3.0a3-py3.6-linux-x86_64.egg/Cython/Build/Dependencies.py", line 1105, in cythonize
    cythonize_one(*args)
  File "/tmp/pip-install-m_w5y09w/h5py/.eggs/Cython-3.0a3-py3.6-linux-x86_64.egg/Cython/Build/Dependencies.py", line 1263, in cythonize_one
    raise CompileError(None, pyx_file)
Cython.Compiler.Errors.CompileError: /tmp/pip-install-m_w5y09w/h5py/h5py/h5a.pyx
----------------------------------------

ERROR: Command errored out with exit status 1: /usr/bin/python3 -u -c ‘import sys, setuptools, tokenize; sys.argv[0] = ‘"’"’/tmp/pip-install-m_w5y09w/h5py/setup.py’“'”‘; file=’“'”‘/tmp/pip-install-m_w5y09w/h5py/setup.py’“'”‘;f=getattr(tokenize, ‘"’“‘open’”’“‘, open)(file);code=f.read().replace(’”‘"’\r\n’“'”‘, ‘"’"’\n’“'”‘);f.close();exec(compile(code, file, ‘"’“‘exec’”’"‘))’ install --record /tmp/pip-record-va1qs7va/install-record.txt --single-version-externally-managed --user --prefix= --compile --install-headers /root/.local/include/python3.6m/h5py Check the logs for full command output.

Thus, I am totally unable to work with any mesh generated with Gmsh.

15

I had the same problem in the FEnics Docker container that I downloaded from the downloads page. To install meshio one would create a different docker image using the process mentioned in the FEnics Containers page.The contents of the Dockerfile using instructions here would look like,

FROM quay.io/fenicsproject/stable:current
USER root
RUN apt-get -qq update && \
    apt-get -y upgrade && \
    apt-get clean && \
    pip install --upgrade pip && \
    pip install meshio[all] && \
    rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*
USER root

However this would give the error

ValueError: Invalid high library version bound (invalid high library version bound)

This is because the HDF5 version that ubuntu-18(container used by Fenics) installs by default is 1.10.0, which can be easily verified by the following command

dpkg -l | grep hdf5
ii  libhdf5-100:amd64                      1.10.0-patch1+docs-4                amd64        Hierarchical Data Format 5 (HDF5) - runtime files - serial version
ii  libhdf5-mpich-100:amd64                1.10.0-patch1+docs-4                amd64        Hierarchical Data Format 5 (HDF5) - runtime files - MPICH2 version
ii  libhdf5-mpich-dev                      1.10.0-patch1+docs-4                amd64        Hierarchical Data Format 5 (HDF5) - development files - MPICH version

But pip install meshio[all] would install h5py-2.9 which would require HDF5-1.10.3 or above.The one h5py version that would work with this is 2.7.1.But installing this using pip gives the following error ,reasons for which I am unable to find out now.

Segmentation fault (core dumped)

So the workaround is installing h5py from ubuntu repository only to maintain continuity of versions.The Dockerfile could look like

FROM quay.io/fenicsproject/stable:current
USER root
RUN apt-get -qq update && \
    apt-get -y upgrade && \
    apt-get clean && \
    apt-get -y install python3-h5py && \
    pip install --upgrade pip && \
    pip install meshio[all] && \
    rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*
USER root

That solved the problem.Obviously, the workaround can be used in cases of direct installation .Just install h5py from official repository and follow usual steps as mentioned

1 Like
16

Dear @Sambeet_Panigrahi,

Thank you very much for clarifying and providing a solution.
I just want to point out that the following command can also serve as a quick but somewhat abrupt solution to the problem (since you avoid checking for the desired version).

Try:

HDF5_DISABLE_VERSION_CHECK=1 python filename.py

It worked just fine for me!

Regards,
Santiago

17

Hello!

I am trying to perform this test but it generates an error that does not recognize the Mesh (). See the code and error below:

import meshio
import numpy as np
msh = meshio.read(“Square.msh”)

line_cells = []
for cell in msh.cells:
if cell.type == “triangle”:
triangle_cells = cell.data
elif cell.type == “line”:
if len(line_cells) == 0:
line_cells = cell.data
else:
line_cells = np.vstack([line_cells, cell.data])

line_data = []
for key in msh.cell_data_dict[“gmsh:physical”].keys():
if key == “line”:
if len(line_data) == 0:
line_data = msh.cell_data_dict[“gmsh:physical”][key]
else:
line_data = np.vstack([line_data, msh.cell_data_dict[“gmsh:physical”][key]])
elif key == “triangle”:
triangle_data = msh.cell_data_dict[“gmsh:physical”][key]

triangle_mesh = meshio.Mesh(points=msh.points, cells={“triangle”: triangle_cells})
line_mesh =meshio.Mesh(points=msh.points,
cells=[(“line”, line_cells)],
cell_data={“name_to_read”:[line_data]})
meshio.write(“mesh.xdmf”, triangle_mesh)

meshio.xdmf.write(“mf.xdmf”, line_mesh)

mesh = Mesh()
with XDMFFile(“mesh.xdmf”) as infile:
infile.read(mesh)
File(“Dolfin_circle_mesh.pvd”).write(mesh)

mvc = MeshValueCollection(“size_t”, mesh, 1)
with XDMFFile(“mf.xdmf”) as infile:
infile.read(mvc, “name_to_read”)
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File(“Dolfin_circle_facets.pvd”).write(mf)

V = FunctionSpace(mesh, ‘P’, 1)

bc1 = DirichletBC(V, Constant(0.0), mf, 1)
bc2 = DirichletBC(V, Constant(10.0), mf, 2)
bc3 = DirichletBC(V, Constant(6.0), mf, 3)
bc4 = DirichletBC(V, Constant(3.0), mf, 4)

bc = [bc1, bc2, bc3, bc4]

Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-15.0)
a = dot(grad(u), grad(v))dx
L = fv*dx

Compute solution

u = Function(V)
solve(a == L, u, bc)
File(“Solution.pvd”).write(u)

Error:

File “/Users/XYZ/Desktop/Gmsh/test/test.py”, line 42, in
mesh = Mesh()

NameError: name ‘Mesh’ is not defined

I found the error:

We needed to import Dolfin!

from dolfin import *

Now it worked!

1 Like
18

I have been having a similar issue not yet solved by current recommendations. I found a workaround I wanted to post here in case others have the same issue.

I have generated a basic .msh file from gmsh and pygmsh according to many tutorials found online. Any .msh file produces this error, but I am specifically using the 2D example from this page.

I have discovered that I will get the h5py / HDF5 version mismatch error whenever I attempt to use meshio.read('fileName.msh') if there was a previous instance of import dolfin or import fenics . While obviously not an ideal solution, my current workaround is to generate .xdmf files in a separate script which uses meshio but does not use fenics or dolfin.

For example, I generate my mesh with the following (taken from the above tutorial):

import pygmsh
resolution = 0.01

# Channel parameters
L = 2.2
H = 0.41
c = [0.2, 0.2, 0]
r = 0.05

# Initialize empty geometry using the build in kernel in GMSH
geometry = pygmsh.geo.Geometry()
# Fetch model we would like to add data to
model = geometry.__enter__()
# Add circle
circle = model.add_circle(c, r, mesh_size=resolution)

# Add points with finer resolution on left side
points = [model.add_point((0, 0, 0), mesh_size=resolution),
          model.add_point((L, 0, 0), mesh_size=5*resolution),
          model.add_point((L, H, 0), mesh_size=5*resolution),
          model.add_point((0, H, 0), mesh_size=resolution)]

# Add lines between all points creating the rectangle
channel_lines = [model.add_line(points[i], points[i+1])
                 for i in range(-1, len(points)-1)]

# Create a line loop and plane surface for meshing
channel_loop = model.add_curve_loop(channel_lines)
plane_surface = model.add_plane_surface(
    channel_loop, holes=[circle.curve_loop])

# Call gmsh kernel before add physical entities
model.synchronize()

volume_marker = 6
model.add_physical([plane_surface], "Volume")
model.add_physical([channel_lines[0]], "Inflow")
model.add_physical([channel_lines[2]], "Outflow")
model.add_physical([channel_lines[1], channel_lines[3]], "Walls")
model.add_physical(circle.curve_loop.curves, "Obstacle")


# We generate the mesh using the pygmsh function `generate_mesh`. Generate mesh returns a `meshio.Mesh`. However, this mesh is tricky to extract physical tags from. Therefore we write the mesh to file using the `gmsh.write` function.
geometry.generate_mesh(dim=2)
import gmsh
gmsh.write("mesh.msh")
gmsh.clear()
geometry.__exit__()


import meshio
mesh_from_file = meshio.read("mesh.msh")

import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh


line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

and run it in the following basic example:

# start
print('start')

# import 
import fenics as fn 
import dolfin
import sys
import numpy
import matplotlib.pyplot as plt
import pygmsh
import gmsh

fullDim = 2 #2D model

mesh = fn.Mesh()
mvc = fn.MeshValueCollection("size_t", mesh, fullDim)
with fn.XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
mf = fn.cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc2 = fn.MeshValueCollection("size_t", mesh, fullDim-1)
with fn.XDMFFile("facet_mesh.xdmf") as infile:
   infile.read(mvc2, "name_to_read")
cf = fn.cpp.mesh.MeshFunctionSizet(mesh, mvc2)


# Define function space
V = fn.FunctionSpace(mesh, 'Lagrange', 1)
#markers = fn.MeshFunction('size_t', mesh, fullDim, cf)
#dx = fn.Measure('dx', domain=mesh, subdomain_data=markers)
#dx = fn.Measure('dx', domain=mesh, subdomain_data=cf)
dx = fn.Measure('dx', domain=mesh, subdomain_data=mf)

# Define boundary condition
bc = fn.DirichletBC(V, fn.Constant(0), 'on_boundary')
bcs = [bc]
bc = fn.DirichletBC(V, fn.Constant(1), cf, 5)
bcs.append(bc)

# some plot setup
plt.figure()
plt.subplot(1,3,1)
fn.plot(mf, title='Sub Domains')


# equation solving
u = fn.TrialFunction(V)
v = fn.TestFunction(V)
a = fn.dot(fn.grad(u), fn.grad(v)) * fn.dx
L = fn.Constant('0') * v * fn.dx
u = fn.Function(V)
fn.solve(a == L, u, bcs)

electric_field = fn.project(-fn.grad(u))
plt.subplot(1,3,2)
fn.plot(u, title='Fields')
plt.subplot(1,3,3)
fn.plot(electric_field, title='Solution')

plt.show()

This solution works because the two steps are done in different files. However, if I attempt to modify the start of the second file so that the .msh to .xdmf conversion occurs in the same file as the fenics computations, the following code will fail with the h5py HDF5 errors:

# start
print('start')

# import 
import fenics as fn 
import dolfin
import sys
import numpy
import matplotlib.pyplot as plt
import pygmsh
import gmsh

fullDim = 2 #2D model

##### this is new ################
import meshio
mesh_from_file = meshio.read("mesh.msh")

import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

######### the rest is unchanged ###############

mesh = fn.Mesh()
mvc = fn.MeshValueCollection("size_t", mesh, fullDim)
with fn.XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
mf = fn.cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc2 = fn.MeshValueCollection("size_t", mesh, fullDim-1)
with fn.XDMFFile("facet_mesh.xdmf") as infile:
   infile.read(mvc2, "name_to_read")
cf = fn.cpp.mesh.MeshFunctionSizet(mesh, mvc2)

## and then the rest

I hope this helps someone else with a similar issue.

1 Like
19

Useful for me! Thanks a lot!