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

mesh
1

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

def halfdisk_MSHmesh3D(disk_radius, mesh_resolution, input_MSHmesh_path): 

    gmsh.initialize()
    gmsh.model.add("halfdisk")
    gmsh.option.setNumber("Mesh.Algorithm", 5) 
    gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_resolution) 
    gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_resolution) 
    hd = [
          gmsh.model.occ.addDisk(0., 0., 0., disk_radius, disk_radius), 
          gmsh.model.occ.addRectangle(-disk_radius, 0., 0., 2*disk_radius, disk_radius) 
         ]
    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()
    mesh = meshio.read(input_MSHmesh_path)

    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): 

    msh = meshio.read(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 
        disk_radius = 1.
        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()
2

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.add_collection(pc)
    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
3

Many thanks for your detailed answer conpierce8, it works.