def METHOD_NAME(self, plot=True):
    vertices = np.array([
        ( 0.0, 0.0), ( 1.0, 0.0), ( 1.0,  1.0), (0.0,  1.0),
        (-1.0, 1.0), (-1.0, 0.0), (-1.0, -1.0), (0.0, -1.0)], dtype=np.float)
    facets = np.array([
        (0, 1), (1, 2), (2, 3), (3, 4), 
        (4, 5), (5, 6), (6, 7), (7, 0)], dtype=np.int)
    subdomain = np.array([
        (1, 0), (1, 0), (1, 0), (1, 0),
        (1, 0), (1, 0), (1, 0), (1, 0)], dtype=np.int)
    domain = HalfEdgeDomain.METHOD_NAME(vertices, facets, subdomain)
    domain.boundary_uniform_refine(n=1)
    mesh = domain.to_halfedgemesh()
    if plot:
        fig = plt.figure()
        axes = fig.gca()
        mesh.add_plot(axes)
        mesh.find_node(axes, color='k', showindex=True)
        print("halfede:")
        for i, val in enumerate(mesh.entity('halfedge')):
            print(i, ":", val)
        fig = plt.figure()
        axes = fig.gca()
        mesh.add_halfedge_plot(axes, showindex=True)
        mesh.find_node(axes, color='g', showindex=True)
        plt.show()