Circular mesh looks not okay

[leonardo]

I am trying to generate a circular mesh with pygmsh. This is the code:

import pygmsh
import gmsh
import meshio

path = "/home/sphere/"

resolution = .0125

# An empty geometry
geometry = pygmsh.geo.Geometry()
# Create a model to add data to
model = geometry.__enter__()

# A circle centered at the origin and radius 1
circle = model.add_circle([0.0, 0.0, 0.0], 1.0, mesh_size=5*resolution)  # meshes are always 3D, I will suppress the third component in case

# Sinchronize, before adding physical entities
model.synchronize()

# Tagging
model.add_physical(circle.curve_loop.curves, "circle")

# Generate the mesh
geometry.generate_mesh(dim=1)
gmsh.write(path+"circle_mesh.msh")
gmsh.clear()
geometry.__exit__()


def create_mesh(mesh, cell_type, prune_z = False):
    cells = mesh.get_cells_type(cell_type)  # get the cells of some type: it will change!
    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

# Using the above function, create line and "plane" mesh
mesh_from_file = meshio.read(path+"circle_mesh.msh")
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write(path+"circle_mesh.xdmf", line_mesh)

What Fenics sees is then:

What ParaView sees is:

How can I get rid of the “inner triangle” and have just a regular circle like Paraview displays it?

[dokken]

What do you mean by

@leonardo

What Fenics sees is then:

Do you mean that if you call plot(mesh) you get this plot?
How does a solution on the geometry look like when exported to Paraview as xdmf or pvd file?

[leonardo]

Yes, I used plot(mesh_circle).

Here is an image of a constantly one function, using paraview, including the mesh too:

It comes from the code:

mesh_path = ...

mesh_circle = Mesh()
with XDMFFile(mesh_path + "circle_mesh.xdmf") as infile:
    infile.read(mesh_circle)

L1_S = FiniteElement("Lagrange", mesh_circle.ufl_cell(), 1)
V_S = FunctionSpace(mesh_circle, L1_S)

u = Function(V_S)
u.vector()[:]=1


vtkfile = File("u.pvd")

vtkfile << u

[leonardo]

As a follow up, here is a coarser version of the mesh:

[dokken]

I think this is simply a bug in plotting the mesh, as both the integration area (circumference) and the function u is correctly projected with the following code:

from pickletools import uint2
import pygmsh
import gmsh
import meshio

path = "./"

resolution = .0125

# An empty geometry
geometry = pygmsh.geo.Geometry()
# Create a model to add data to
model = geometry.__enter__()

# A circle centered at the origin and radius 1
circle = model.add_circle([0.0, 0.0, 0.0], 1.0, mesh_size=5*resolution)  # meshes are always 3D, I will suppress the third component in case

# Sinchronize, before adding physical entities
model.synchronize()

# Tagging
model.add_physical(circle.curve_loop.curves, "circle")

# Generate the mesh
geometry.generate_mesh(dim=1)
gmsh.write(path+"circle_mesh.msh")
gmsh.clear()
geometry.__exit__()


def create_mesh(mesh, cell_type, prune_z = False):
    cells = mesh.get_cells_type(cell_type)  # get the cells of some type: it will change!
    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

# Using the above function, create line and "plane" mesh
mesh_from_file = meshio.read(path+"circle_mesh.msh")
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write(path+"circle_mesh.xdmf", line_mesh)


from dolfin import *

mesh = Mesh()
with XDMFFile("circle_mesh.xdmf") as xdmf:
    xdmf.read(mesh)

x, y = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, "CG", 1)
u =project(x+y, V)
File("u.pvd")<< u
import matplotlib.pyplot as plt
plt.savefig("mesh.png")
print(assemble(1*dx(domain=mesh)))