Create periodic mesh of a sector geometry

mesh
1

I have a sector geometry where the sector boundaries should be periodic. I am using routines defined in script to turn the mesh periodic. First challenge I encounter is that script.create_periodic_mesh does not work with 0.10.0 (conda). Switching to 0.9.0 allows me to create periodic mesh. However, then I am unable to port cell tags from the original to periodic mesh because of duplicate mesh entity. Facet tags get ported fine though.

I will appreciate if someone could suggest a fix ideally for 0.10. Below is the MWE that reproduces the problem on both 0.9 and 0.10:

import dolfinx, gmsh, ufl, mpi4py, petsc4py, script
import numpy as np

sec_angle = np.pi/6
R = 1

gmsh.finalize()
gmsh.initialize()
# gmsh.option.setNumber("General.Terminal", 0)
gmsh.clear()

occ = gmsh.model.occ
gdim = 2

sec = occ.addCircle(0, 0, 0, R, angle1=np.pi/2-sec_angle/2, angle2=np.pi/2+sec_angle/2)
sec_pnts = occ.get_entities(0)
p0 = occ.addPoint(0, 0, 0)
p1 = occ.addPoint(0, R, 0)
l1 = occ.addLine(sec_pnts[0][1], p0)
l2 = occ.addLine(sec_pnts[1][1], p0)
curve = occ.addCurveLoop([sec, l1, l2])
surf = occ.addPlaneSurface([curve])
l3 = occ.addLine(p0, p1)
frag = occ.fragment([(gdim, surf)], [(gdim-1, l3)])
occ.removeAllDuplicates()
occ.synchronize()

all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

# sector boundary tags by visual inspection
left_edge_tag = (7,)
right_edge_tag = (5,)
gmsh.model.mesh.setPeriodic(gdim-1, left_edge_tag, right_edge_tag,
                            [np.cos(sec_angle), -np.sin(sec_angle), 0, 0, 
                            np.sin(sec_angle), np.cos(sec_angle), 0, 0, 
                            0, 0, 1, 0,
                            0, 0, 0, 1])

gmsh.model.mesh.generate(gdim)

mpi_rank = 0
# dolfinx_model = dolfinx.io.gmsh.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim)
if dolfinx.__version__ == "0.9.0":
    from dolfinx.io import gmshio
    dolfinx_model = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim)
    mesh0, ct0, ft0 = dolfinx_model[0], dolfinx_model[1], dolfinx_model[2]
if dolfinx.__version__ == '0.10.0':
    dolfinx_model = dolfinx.io.gmsh.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim)
    mesh0, ct0, ft0 = dolfinx_model.mesh, dolfinx_model.cell_tags, dolfinx_model.facet_tags

# create periodic mesh
def indicator(x):
    angle = np.arctan2(x[1], x[0])
    edge_indices = np.isclose(angle, np.pi/2-sec_angle/2) # right edge
    return edge_indices

def mapping(x):
    rot_matrix = np.array([[np.cos(sec_angle), -np.sin(sec_angle)], 
                           [np.sin(sec_angle), np.cos(sec_angle)]])
    values = np.zeros(x.shape)
    values[:2,:] = rot_matrix @ x[:2]
    return values

# works with 0.9.0 but not 0.10.0
mesh, replaced_vertices, replacement_map = script.create_periodic_mesh(mesh0, indicator, mapping)

ct = script.transfer_meshtags_to_periodic_mesh(
        mesh0, mesh, replaced_vertices, ct0
    )
ft = script.transfer_meshtags_to_periodic_mesh(
        mesh0, mesh, replaced_vertices, ft0
    )