Unexpected distribution calculating CellDiameter using GMSH mesh

Hi all,

When trying to calculate the CellDiameter from a mesh made using GMSH, I get a unexpected distribution.

My mesh is made using the following script:

# Import modules:
import gmsh
import sys
nodes = 102
nodesy = int(102/3)

# Initialize gmsh:
# cube points:
lc = 0.25
point1 = gmsh.model.geo.add_point(0, 0, 0, lc) # Lower left corner
point2 = gmsh.model.geo.add_point(1, 0, 0, lc) # Lower right corner
point3 = gmsh.model.geo.add_point(1, 0.33, 0, lc) # Outlet start
point4 = gmsh.model.geo.add_point(1, 0.66, 0 ,lc) # Outlet end
point5 = gmsh.model.geo.add_point(1, 1, 0, lc) # Upper right corner
point6 = gmsh.model.geo.add_point(0, 1, 0, lc) # Upper left corner
point7 = gmsh.model.geo.add_point(0, 0.66, 0, lc) # Upper inlet
point8 = gmsh.model.geo.add_point(0, 0.33, 0, lc) # Lower inlet

line1 = gmsh.model.geo.add_line(point1, point2) # Bottom wall
line2 = gmsh.model.geo.add_line(point2, point3) # Lower right wall
line3 = gmsh.model.geo.add_line(point3, point4) # Outlet
line4 = gmsh.model.geo.add_line(point4, point5) # Upper right wall
line5 = gmsh.model.geo.add_line(point5, point6) # Top wall
line6 = gmsh.model.geo.add_line(point6, point7) # Upper left wall
line7 = gmsh.model.geo.add_line(point7, point8) # Inlet
line8 = gmsh.model.geo.add_line(point8, point1) # Lower left wall

face1 = gmsh.model.geo.add_curve_loop([line1, line2, line3, line4, line5, line6, line7, line8]) 

surface1 = gmsh.model.geo.add_plane_surface([face1])

gmsh.model.geo.add_physical_group(1, [line1, line2, line4, line5, line6, line8], 1, "Walls")
gmsh.model.geo.add_physical_group(1, [line3], 3, "Outlet")
gmsh.model.geo.add_physical_group(1, [line7], 2, "Inlet")

gmsh.model.geo.add_physical_group(2, [surface1], 4, "Surface")

gmsh.model.geo.mesh.setTransfiniteCurve(line1, nodes)
gmsh.model.geo.mesh.setTransfiniteCurve(line2, nodesy)
gmsh.model.geo.mesh.setTransfiniteCurve(line3, nodesy)
gmsh.model.geo.mesh.setTransfiniteCurve(line4, nodesy)
gmsh.model.geo.mesh.setTransfiniteCurve(line5, nodes)
gmsh.model.geo.mesh.setTransfiniteCurve(line6, nodesy)
gmsh.model.geo.mesh.setTransfiniteCurve(line7, nodesy)
gmsh.model.geo.mesh.setTransfiniteCurve(line8, nodesy)

gmsh.model.geo.mesh.setTransfiniteSurface(surface1, "left", [point1, point2, point5, point6])

# Create the relevant Gmsh data structures
# from Gmsh model.
# Generate mesh:
# Write mesh data:
# Creates  graphical user interface
if 'close' not in sys.argv:
# It finalize the Gmsh API

Which generates this mesh (from the picture the mesh looks uneven, but is a mapped mesh):

My MWE is as follows:

from fenics import *
import matplotlib.pyplot as plt
import meshio

msh = meshio.read("test.msh")

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)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
h = CellDiameter(mesh) # Diameter of each element

Which results in this plot of the celldiameter:

Were the celldiameter is not distributed as expected.

The distribution is similar if you use the CellSize filter in Paraview on your mesh.xdmf file. Note that the interface is sharper in Paraview as it computes it per cell (DG-0), while the plot function projects into first order Lagrange elements, making a sharp interface.

If you use the following code:

h = CellDiameter(mesh) # Diameter of each element
H = FunctionSpace(mesh, "DG", 0)
h_ = project(h, H)
with XDMFFile("cd.xdmf") as xdmf:
    xdmf.write_checkpoint(h_, "h", 0.0, append=False)

you get a sharp interface.