Error when importing msh from gmsh

HI, I was using a code written in 2021 by importing mesh from fenics and converting it to xml format using meshio-convert. Since I wanted to move to using xdmf format I used the following code from this forum.

import meshio
msh = meshio.read("Rec_geo.msh")
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif  cell.type == "tetra":
        tetra_cells = cell.data

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
#tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
#meshio.write("Rec_geo.xdmf", tetra_mesh)

meshio.write("Rgeo.xdmf", triangle_mesh)

Now I used the following code to import the mesh and tried to test if I can run.

from math import degrees
from operator import not_
import dolfin as df
from dolfin.fem.interpolation import interpolate
from dolfin.function.constant import Constant
import matplotlib.pyplot as plt
import numpy as np
import os

from numpy.core.fromnumeric import diagonal
from numpy.core.numeric import True_
from numpy.lib.twodim_base import _trilu_indices_form_dispatcher
import ufl as ufl


model = "plane_stress"
E = df.Constant(110000.)
nu = df.Constant(0.3)
mu = E / (2 * (1 + nu))
lmbda = (E * nu) / ((1 + nu) * (1 - 2 * nu))
if model == "plane_stress":
    lmbda = 2 * mu * lmbda / (lmbda + 2 * mu)
def eps(v):
    return df.sym(df.grad(v))
def sigmaE(v, lmbda, mu):
    sigE = lmbda * df.tr(eps(v)) * df.Identity(2) + 2.0 * mu * eps(v)
    return sigE
def sigmaA(v, sigA, lmbda, mu):
    sig = sigA + lmbda * df.tr(eps(v)) * df.Identity(2) + 2.0 * mu * eps(v)
    return sig
def load_mesh(mesh_file):
    mesh = df.Mesh()
    with df.XDMFFile(mesh_file) as infile:
        infile.read(mesh)
        # These are the markers
        ffun = df.MeshFunction("size_t", mesh, 2)
        infile.read(ffun, "name_to_read")
    return mesh, ffun
#mesh, ffun = load_mesh('Rgeo.xdmf')
mesh = df.RectangleMesh(df.Point(0,0),df.Point(10,5),10,5,diagonal="right")
V0 = df.VectorFunctionSpace(mesh, "Lagrange", degree=2)#,dim=2)#,form_degree=3)
V0sig = df.TensorFunctionSpace(mesh, 'DG', degree=2)#,shape=(2,2))
def sigA_init_turning(V,Coef_sig_A):
    ep = np.array(V.tabulate_dof_coordinates())[:, 1].max() - np.array(V.tabulate_dof_coordinates())[:, 1].min()
    x_1 = [0,0.015,0.025,0.035,0.05,0.065,0.075,0.08,0.085,0.09,0.095,0.1,0.105,0.11,0.115,0.12,0.125,0.13,0.135]
    y = [500,50,-250,-330,-300,-210,-210,-210,-210,-210,-210,-210,-210,-210,-210,-210,-210,-210,-210]
    fL= np.polyfit(x_1,y,7)
    p = np.poly1d(fL)
    s1=fL[0]
    s2=fL[1]
    s3=fL[2]
    s4=fL[3]
    s5=fL[4]
    s6=fL[5]
    s7=fL[6]
    s8=fL[7]
    depth = df.Expression('ep - x[1]', degree=1, ep=ep)
    #sig_xx = df.Expression('x[1]^7 + s1',degree = 1,s1=s1)
    sig_xx = df.Expression('s1 * y * y*y*y*y*y*y+ s2 * y*y*y*y*y*y\
     + s3 * y*y*y*y*y + s4 * y*y*y*y + s5 * y*y*y + s6* y*y + s7 * y + s8'\
    ,degree=1,s1=s1,s2=s2,s3=s3,s4=s4,s5=s5,s6=s6,s7=s7,s8=s8,y=depth)
    sigA = df.Expression((('d <= d1 ? sigA_xx : 0.0', 'd <= d1 ? 0. : 0.'), \
                          ('d <= d1 ? 0. : 0.', 'd <= d1 ? 0. : 0.')), \
                         degree=1, ep=ep, d=depth, d1=0.135,sigA_xx=sig_xx)
    
    return sigA
sigfun = sigA_init_turning(V0,0)
df.project(sigfun,V0sig)
du = df.TrialFunction(V0)
u_ = df.TestFunction(V0)
val=sigmaE(du, lmbda, mu)
#a = df.inner(sigmaE(du, lmbda, mu), eps(u_)) * df.dx
value=(df.grad(du))
print((value))

Now I receive the following error.

Invalid MIT-MAGIC-COOKIE-1 keyShapes do not match: <Argument id=140057592358176> and <Coefficient id=140057551284608>.
Traceback (most recent call last):
  File "testing.py", line 71, in <module>
    df.project(sigfun,V0sig)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py", line 129, in project
    L = ufl.inner(w, v) * dx
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 158, in inner
    return Inner(a, b)
  File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 147, in __new__
    error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: <Argument id=140057592358176> and <Coefficient id=140057551284608>.

I also tried enforcing the shape of the element.(these portions are commented), but they don’t seem to work either. If I use the mesh generated using df.rectangularMesh() there seems to be no error.
I found a workaround by defining the gradient in 2D which I would like to avoid to make the program more abstract. It would be helpful if someone can explain why the error happenes when I use *.xdmf for mesh input.

PS : mesh files generated in 2021 using gmsh the same method worked very well.

You need to remove the third dimension of msh.points as done in:

1 Like