Cracked Domain mesh

I understand well this is not a PyGmsh forum, but I am taking the liberty to ask a meshing question as some users here might have the solution ready on their desk.

I would like to use Dolfinx to solve a Laplacian-type equation on a cracked domain like the one pictured below.

import gmsh
gmsh.initialize()
gmsh.model.add("Cracked_plate")
## SQ stands for size quarter
SQ = 1

p1 = gmsh.model.occ.addPoint(0, 0, 0, 0, 1)
p2 = gmsh.model.occ.addPoint(0, SQ, 0, 0, 2)
p3 = gmsh.model.occ.addPoint(SQ, SQ, 0, 0, 3)
p4 = gmsh.model.occ.addPoint(SQ, - SQ, 0, 0, 4)
p5 = gmsh.model.occ.addPoint(-SQ, - SQ, 0, 0, 5)
p6 = gmsh.model.occ.addPoint(-SQ, SQ, 0, 0, 6)
## Reusing p2 instead of p7
#p7 = gmsh.model.occ.addPoint(0, SQ, 0, 0, 7)
#p8 = gmsh.model.occ.addPoint(0, 0, 0, 8)


# Creating connection lines
l1 = gmsh.model.occ.addLine(1, 2, 10)
l2 = gmsh.model.occ.addLine(2, 3, 11)
l3 = gmsh.model.occ.addLine(3, 4, 12)
l4 = gmsh.model.occ.addLine(4, 5, 13)

l5 = gmsh.model.occ.addLine(5, 6, 14)
l6 = gmsh.model.occ.addLine(6, 2, 15)
l7 = gmsh.model.occ.addLine(2, 1, 16)
#l8 = gmsh.model.occ.addLine(8, 1, 17)


loop_outer = gmsh.model.occ.addCurveLoop([l1,l2,l3,l4, l5, l6, l7], 30)
gmsh.model.occ.addPlaneSurface([30], 100)

gmsh.model.occ.synchronize()

gdim = 2
gmsh.model.addPhysicalGroup(100, [20], 400)

For the life of me I cannot get it to work. Points 2 and 7 share the same physical location, I wondered if that could be the problem, tried both options, to have them as separate points or re-use the first defined.

I looked all over the internet to no avail. I found mentions to an “embed” method, as well as references to a crack plugin on this forum, but it seemed incredibly laborious to set up!

I could get it to "work” by adding a small separation between points 2 and 7

eps = 0.01
p1 = gmsh.model.occ.addPoint(0, 0, 0, 0, 1)
p2 = gmsh.model.occ.addPoint(0+eps, SQ, 0, 0, 2)
p3 = gmsh.model.occ.addPoint(SQ, SQ, 0, 0, 3)
p4 = gmsh.model.occ.addPoint(SQ, - SQ, 0, 0, 4)
p5 = gmsh.model.occ.addPoint(-SQ, - SQ, 0, 0, 5)
p6 = gmsh.model.occ.addPoint(-SQ, SQ, 0, 0, 6)
## Reusing p2 instead of p7
p7 = gmsh.model.occ.addPoint(0-eps, SQ, 0, 0, 7)

but I am interested in the singular field at the tip and would prefer to have the exact geometry.

If anybody had the patience to point me in the right direction that would be great, thanks a lot