simple mesh doesn't work with fenics at certain characteristic length

[quan4444]

I have an issue presented below. Thank you for helping!!

Mesh: a box with a cylinder inside. The box and cylinder are in 2 different physical groups. Mesh size is determined by changing characteristic length. The code is posted below

Process: Generate mesh in pygmsh, export in '.xdmf' and '.h5' file type. Import to fenics and perform equibiaxial extension.

Issue: This process works fine for certain characteristic length (e.g. 3,2,1.5,1,0.85,0.76). However, fenics gives the error "Unable to compute matrix factorization | The provided data did not satisfy the prerequisites" for other characteristic length (e.g. 0.80,0.79,0.65,0.64,0.5).

import pygmsh
import gmsh
import meshio
import numpy as np

# inputs
L=W=40.0 # length and width of box
D=2.0 # depth of box
R=10.0 # radius of circle
mesh_size = 0.79 # change the characteristic length

# initialize pygmsh.occ.Geometry()
geom = pygmsh.occ.Geometry()
model = geom.__enter__()
model.characteristic_length_min=mesh_size
model.characteristic_length_max=mesh_size

# create 3D shapes
cyl = model.add_cylinder(x0=[L/2,W/2,0.0],axis=[0.0,0.0,D],radius=R)
thin_samp = model.add_box(x0=[0.0,0.0,0.0],extents=[L,W,D])
thin_samp_cyl = model.boolean_fragments(thin_samp,cyl)

# label the domains
domain0 = geom.add_physical(thin_samp_cyl[1],label='0')
domain1 = geom.add_physical(thin_samp_cyl[0],label='1')

mesh = geom.generate_mesh(dim=3)
mesh_name = 'output_pygmsh/thin_samp_cyl_ms'+str(mesh_size)
gmsh.write(mesh_name+'.msh')
mesh.write(mesh_name+'.vtk')

# export to .xdmf
msh = meshio.read(mesh_name+'.msh')

tetra_cells = np.array([None])
for cell in msh.cells:
    if cell.type == 'tetra':
        # collect individual meshes
        if tetra_cells.all() == None:
            tetra_cells = cell.data
        else:
            tetra_cells = np.concatenate((tetra_cells,cell.data))
        
tetra_data = []
for key in msh.cell_data_dict['gmsh:physical'].keys():
    if key == 'tetra':
        tetra_data = msh.cell_data_dict['gmsh:physical'][key]

tetra_mesh = meshio.Mesh(points=msh.points, cells={'tetra':tetra_cells},\
                        cell_data={'regions':[tetra_data]})

meshio.write(mesh_name+'.xdmf',tetra_mesh)

[wraith1995]

Have you tried generating a mesh like this from the usual gmsh API and seeing what happens? I think you need to determine if this a gmsh problem or a fenics problem.