Error exporting field to XDMF

dolfinx
1

The following program solves the equation fine when using dolfinx.mesh.create_rectangle and crushes when importing XDMF mesh generated by gmsh. Is it a bug, or did I do something wrong?

Steps to reproduce:

  1. Uncomment create_rectangle code section in program.py and run
    python3 program.py
    Observe it finishing successfuly.
  2. Generate and convert gmsh mesh
    gmsh mesh.geo -2
python3 convert_msh.py mesh.msh mesh.xdmf
  3. Uncomment gmsh code section in program.py and run
    python3 program.py

Observe it fail with corrupted double-linked list error.

Build commits:
DolfinX: Simplify Spack instructions (#2146) · FEniCS/dolfinx@2f63364 · GitHub
UFL: Signature computation speedup (#94) · FEniCS/ufl@2ad4ef1 · GitHub
FFCX: Update variants (#476) · FEniCS/ffcx@3d440ca · GitHub
BasiX: Make custom element demo clearer (#440) · FEniCS/basix@83aed93 · GitHub

Full error message:

corrupted double-linked list
[mypc:26136] *** Process received signal ***
[mypc:26136] Signal: Aborted (6)
[mypc:26136] Signal code:  (-6)
[mypc:26136] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0x12200)[0x7f4073053200]
[mypc:26136] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x141)[0x7f4072e9e8a1]
[mypc:26136] [ 2] /lib/x86_64-linux-gnu/libc.so.6(abort+0x112)[0x7f4072e88546]
[mypc:26136] [ 3] /lib/x86_64-linux-gnu/libc.so.6(+0x79eb8)[0x7f4072edfeb8]
[mypc:26136] [ 4] /lib/x86_64-linux-gnu/libc.so.6(+0x8191a)[0x7f4072ee791a]
[mypc:26136] [ 5] /lib/x86_64-linux-gnu/libc.so.6(+0x827fc)[0x7f4072ee87fc]
[mypc:26136] [ 6] /lib/x86_64-linux-gnu/libc.so.6(+0x82bf5)[0x7f4072ee8bf5]
[mypc:26136] [ 7] /lib/x86_64-linux-gnu/libc.so.6(cfree+0x64)[0x7f4072eec9b4]
[mypc:26136] [ 8] /home/user/Projects/DolfinX/build-dolfinx/lib/lib/libdolfinx.so.0.3(_ZN2xt17xtensor_containerINS_7uvectorIdSaIdEEELm4ELNS_11layout_typeE1ENS_22xtensor_expression_tagEED1Ev+0x24)[0x7f407245c044]
[mypc:26136] [ 9] /home/user/Projects/DolfinX/build-dolfinx/lib/lib/libdolfinx.so.0.3(_ZNK7dolfinx3fem8FunctionISt7complexIdEE4evalERKN2xt17xtensor_containerINS5_7uvectorIdSaIdEEELm2ELNS5_11layout_typeE1ENS5_22xtensor_expression_tagEEERKN3tcb4spanIKiLln1EEERNS6_INS7_IS3_SaIS3_EEELm2ELSA_1ESB_EE+0xa25)[0x7f4072548a15]
[mypc:26136] [10] /home/user/Projects/DolfinX/build-dolfinx/lib/lib/libdolfinx.so.0.3(+0x1b9a40)[0x7f407253aa40]
[mypc:26136] [11] /home/user/Projects/DolfinX/build-dolfinx/lib/lib/libdolfinx.so.0.3(_ZN7dolfinx2io10xdmf_utils21get_point_data_valuesERKNS_3fem8FunctionISt7complexIdEEE+0xa)[0x7f407253b2fa]
[mypc:26136] [12] /home/user/Projects/DolfinX/build-dolfinx/lib/lib/libdolfinx.so.0.3(+0x1a17c3)[0x7f40725227c3]
[mypc:26136] [13] /home/user/Projects/DolfinX/build-dolfinx/lib/lib/libdolfinx.so.0.3(+0x1948a7)[0x7f40725158a7]
[mypc:26136] [14] /home/user/.local/lib/python3.9/site-packages/dolfinx/cpp.cpython-39-x86_64-linux-gnu.so(+0x191841)[0x7f407278e841]
[mypc:26136] [15] /home/user/.local/lib/python3.9/site-packages/dolfinx/cpp.cpython-39-x86_64-linux-gnu.so(+0x71583)[0x7f407266e583]
[mypc:26136] [16] python3[0x5414f7]
[mypc:26136] [17] python3(_PyObject_MakeTpCall+0x26b)[0x52241b]
[mypc:26136] [18] python3[0x53e216]
[mypc:26136] [19] python3(_PyEval_EvalFrameDefault+0x5441)[0x51bf71]
[mypc:26136] [20] python3[0x5158a1]
[mypc:26136] [21] python3(_PyFunction_Vectorcall+0x3e3)[0x52d243]
[mypc:26136] [22] python3(_PyEval_EvalFrameDefault+0x72b)[0x51725b]
[mypc:26136] [23] python3[0x5158a1]
[mypc:26136] [24] python3(_PyEval_EvalCodeWithName+0x47)[0x515617]
[mypc:26136] [25] python3(PyEval_EvalCode+0x23)[0x5155c3]
[mypc:26136] [26] python3[0x62fc37]
[mypc:26136] [27] python3[0x62c4d0]
[mypc:26136] [28] python3[0x62f3b9]
[mypc:26136] [29] python3(PyRun_SimpleFileExFlags+0x193)[0x62ede3]
[mypc:26136] *** End of error message ***
 26136 IOT instruction  python3 program.py

DolfinX program program.py:

from dolfinx import (io, fem, mesh, cpp)
import ufl
import numpy as np
from mpi4py import MPI
from slepc4py import SLEPc
from petsc4py.PETSc import ScalarType


## create_rectangle version

# msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
#                             points=((0, 0), (1, 1)),
#                             n=(40, 40),
#                             cell_type=mesh.CellType.triangle)

#######

## gmsh version

with io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    msh = xdmf.read_mesh(name="Grid")

#######

msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)

N1curl = ufl.FiniteElement("Nedelec 1st kind H(curl)", msh.ufl_cell(), 1)
V = fem.FunctionSpace(msh, N1curl)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

a = ufl.inner(ufl.curl(u), ufl.curl(v)) * ufl.dx
b = ufl.inner(u, v) * ufl.dx

bc_facets = np.where(np.array(cpp.mesh.compute_boundary_facets(msh.topology)) == 1)[0]
bc_dofs = fem.locate_dofs_topological(V, msh.topology.dim-1, bc_facets)

u_bc = fem.Function(V)
u_bc.x.array[:] = 0
bc = fem.dirichletbc(u_bc, bc_dofs)

A = fem.petsc.create_matrix(fem.form(a))
fem.petsc.assemble_matrix(A, fem.form(a), bcs=[bc])
A.assemble()
B = fem.petsc.create_matrix(fem.form(b))
fem.petsc.assemble_matrix(B, fem.form(b), bcs=[bc])
B.assemble()

eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
st = eps.getST()
st.setType(SLEPc.ST.Type.SINVERT)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
eps.setTarget(5.5)
eps.setDimensions(20)
eps.solve()

vals = [(i, eps.getEigenvalue(i)) for i in range(eps.getConverged())]
vals.sort(key=lambda x: x[1].real)

j = 0
E = fem.Function(V)
with io.XDMFFile(MPI.COMM_WORLD, "out.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
    for i, _ in vals:
        print(f"eigenvalue: {eps.getEigenpair(i, E.vector).real:.12f}")
        E.name = f"E-{j+1:03d}-{eps.getEigenpair(i, E.vector).real:.4f}"
        xdmf.write_function(E)
        j += 1

msh to xdmf conversion script convert_msh.py:

from sys import argv
import meshio

def convert_msh(input_file: str, output_file: str,
                cell_type="triangle", cell_data_tag="gmsh:physical"):
    msh = meshio.read(input_file)
    cells = msh.get_cells_type(cell_type)
    if cells.size == 0:
        raise AttributeError(f"There's no {cell_type} data in {input_file}")
    regions = msh.get_cell_data(cell_data_tag, cell_type)
    out_mesh = meshio.Mesh(points=msh.points,
                           cells={cell_type: cells},
                           cell_data={"regions": [regions]})
    meshio.write(output_file, out_mesh, file_format="xdmf")

if __name__ == "__main__":
    if len(argv) < 3 or len(argv) > 4:
        exit("Usage:\t./plot mesh.msh mesh.xdmf [[line|triangle]]")

    region = "triangle"
    if len(argv) == 4:
        region = argv[3].lower()

    convert_msh(argv[1], argv[2], cell_type=region)

gmsh script mesh.geo:

val_pi = 3.1415926535897932384626433832795;

Point(1) = {0, 0, 0, 0.1};
Point(2) = {0, val_pi, 0, 0.1};
Point(3) = {val_pi, val_pi, 0, 0.1};
Point(4) = {val_pi, 0, 0, 0.1};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line loop(1) = {1,2,3,4};
Plane Surface(1) = {1};
Physical Surface(1) = {1};
Physical Line(1) = {1,2,4};
Physical Line(2) = {3};
2

The problem is that you are sending in a 3D mesh, as you are not pruning the z-values when converting the msh to xdmf.
I.e.

should be

  out_mesh = meshio.Mesh(points=msh.points[:, :2],
                           cells={cell_type: cells},
                           cell_data={"regions": [regions]})

and then the code runs as expected.
Note that you can check the geometrical dimension of your mesh by printing

print(msh.geometry.dim)

I’ll have a look at the xdmf issue with manifolds tomorrow.

3

Thanks, now it works!