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:
gmsh.initialize()
# 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.
gmsh.model.geo.synchronize()
# Generate mesh:
gmsh.model.mesh.generate(dim=2)
# Write mesh data:
gmsh.write("test.msh")
# Creates graphical user interface
if 'close' not in sys.argv:
gmsh.fltk.run()
# It finalize the Gmsh API
gmsh.finalize()
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(mesh)
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
plt.colorbar(plot(h))
plt.show()
Which results in this plot of the celldiameter:
Were the celldiameter is not distributed as expected.
