# How to compute boundary-mesh and submesh from an halfdisk-mesh?

Hi everyone,

I am trying to compute both boundary mesh and a submesh of this boundary mesh from an halfdisk-mesh. I am using gmsh to build the halfdisk mesh.

But plotting the boundary mesh (fenics_halfdisk_bmesh) and part of the boundary mesh (fenics_halfdisk_bmesh_CurvedBoundary), transversal lines appear and I am not sure doing it properly:

• Why do I get transversal lines?
• Could you tell me how to get “well organized” submeshes ?
• Also, when I move the tolerance in the CurvedBoundary() class to lower (e.g. 1E-5) , the halfdisk_bmesh_CurvedBoundary does not contain anything anymore, could you tell me why? / How to choose idoine tolerance in that case?

Please find below the code to reproduce:

``````from fenics import *
import matplotlib.pyplot as plt
import gmsh
import meshio
import numpy as np
import os

gmsh.initialize()
gmsh.option.setNumber("Mesh.Algorithm", 5)
gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_resolution)
gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_resolution)
hd = [
]
gmsh.model.occ.intersect([(2, hd[0])], [(2, hd[1])])
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(2)
gmsh.write(input_MSHmesh_path)
#gmsh.fltk.run()
gmsh.finalize()

return mesh

def mesh3D_to_mesh2D(mesh, cell_type, prune_z=False):
# code source : https://fenicsproject.discourse.group/t/how-to-use-meshvaluecollection/5106

cells = np.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
cell_data = np.hstack([mesh.cell_data_dict["gmsh:geometrical"][key]
for key in mesh.cell_data_dict["gmsh:geometrical"].keys() if key==cell_type])
# Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate
points= mesh.points
if prune_z:
points = points[:,:2] # remove z
mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})

return mesh_new

def convert_msh_to_xml(MSHmesh_path):

A = str(os.path.splitext(MSHmesh_path)[0])
B = os.path.splitext(A)[0]
XMLmesh_output_path = str(B + ".xml")
print('xml path = {}'.format(XMLmesh_output_path))
meshio.write(XMLmesh_output_path, meshio.Mesh(points = msh.points, cells = {'triangle': msh.cells_dict['triangle']}))

return XMLmesh_output_path

########
# MAIN #
########

# 2D halfdisk mesh .msh
halfdisk_mesh3D = halfdisk_MSHmesh3D(1., 0.03, "./halfdisk_mesh3D.msh")

# Transform 3D .msh into 2D .msh
halfdisk_mesh2D = mesh3D_to_mesh2D(halfdisk_mesh3D, "triangle", prune_z=True)
meshio.write("./halfdisk_mesh2D.msh", halfdisk_mesh2D)

# Transform 2D .msh into FEniCS readable mesh formats (.xml and Mesh())
XMLmesh2D_output_path = convert_msh_to_xml("./halfdisk_mesh2D.msh")
fenics_halfdisk_mesh2D = Mesh(XMLmesh2D_output_path)

# Plot boundary submesh
fenics_halfdisk_bmesh = BoundaryMesh(fenics_halfdisk_mesh2D, "exterior")
plot(fenics_halfdisk_bmesh, title='Boundaries submesh')
plt.show()

# Plot the circular part of the boundary submesh
class CurvedBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-3
return near(sqrt(x[0] * x[0] + x[1] * x[1]), disk_radius, tol) #and on_boundary

fenics_halfdisk_bmesh_CurvedBoundary = SubMesh(fenics_halfdisk_bmesh, CurvedBoundary())
plot(fenics_halfdisk_bmesh_CurvedBoundary, title='Curved Boundary submesh')
plt.show()

``````

I think this is just an error in the `plot` function. If I replace the `plot` command with `plot_mesh` (where `plot_mesh` is defined as follows), there are no spurious lines crossing the mesh.

``````def plot_mesh(msh, title):
plt.figure()
ax = plt.gca()
coords = msh.coordinates()
cells = msh.cells()

vertices = [coords[k,:] for j in range(cells.shape[0]) for k in cells[j,:]]
codes = [code for j in range(cells.shape[0]) for code in [Path.MOVETO, Path.LINETO]]

path = Path(vertices, codes)
pc = PathCollection([path], edgecolors=["k"])

ax.set_xlim([-1.1, 1.1])
ax.set_ylim([-0.1, 1.1])
ax.set_title(title)
``````

EDIT: Note that you need to `from matplotlib.path import Path` and `from matplotlib.collections import PathCollection` for the above function to work.

DOLFIN checks the midpoint of the element as well as the endpoints. Most likely the midpoint lies more than `1e-5` inside the curved boundary. If you wanted to disable checking the midpoint, you could use the following to create the submesh on the curved boundary:

``````boundary_marker = MeshFunction('size_t', fenics_halfdisk_bmesh, fenics_halfdisk_bmesh.topology().dim(), 0)
CurvedBoundary().mark(boundary_marker, 1, check_midpoint=False)
fenics_halfdisk_bmesh_CurvedBoundary = SubMesh(fenics_halfdisk_bmesh, boundary_marker, 1)
``````
1 Like