How to use the meshzoo.rectangle.create_mesh function in meshzoo

def test_singular_perturbation(self):
        import meshzoo
        from scipy.sparse import linalg

        # Define the problem
        class Poisson(LinearFvmProblem):
            @staticmethod
            def apply(u):
                return integrate(lambda x: - 1.0e-2 * n_dot_grad(u(x)), dS) \
                       + integrate(lambda x: u(x), dV) \
                       - integrate(lambda x: 1.0, dV)
            dirichlet = [(lambda x: 0.0, ['Boundary'])]

        # Create mesh using meshzoo
        vertices, cells = meshzoo.rectangle.create_mesh(
                0.0, 1.0, 0.0, 1.0,
                21, 21,
                zigzag=True
                )
        mesh = pyfvm.meshTri.meshTri(vertices, cells)

        linear_system = pyfvm.discretize(Poisson, mesh)

        x = linalg.spsolve(linear_system.matrix, linear_system.rhs)

        k0 = -1
        for k, coord in enumerate(mesh.node_coords):
            # print(coord - [0.5, 0.5, 0.0])
            if numpy.linalg.norm(coord - [0.5, 0.5, 0.0]) < 1.0e-5:
                k0 = k
                break

def is_inside(self, x): return x[1] >= 0.5
            is_boundary_only = True

        # Define the problem
        class Poisson(LinearFvmProblem):
            @staticmethod
            def apply(u):
                return integrate(lambda x: -n_dot_grad(u(x)), dS) \
                       - integrate(lambda x: 1.0, dV)
            dirichlet = [
                    (lambda x: 0.0, ['Gamma0']),
                    (lambda x: 1.0, ['Gamma1'])
                    ]

        # Create mesh using meshzoo
        vertices, cells = meshzoo.rectangle.create_mesh(
                0.0, 1.0, 0.0, 1.0,
                21, 21,
                zigzag=True
                )
        mesh = pyfvm.meshTri.meshTri(vertices, cells)
        mesh.mark_subdomains([Gamma0(), Gamma1()])

        linear_system = pyfvm.discretize(Poisson, mesh)

        x = linalg.spsolve(linear_system.matrix, linear_system.rhs)

        k0 = -1
        for k, coord in enumerate(mesh.node_coords):
            # print(coord - [0.5, 0.5, 0.0])
            if numpy.linalg.norm(coord - [0.5, 0.5, 0.0]) < 1.0e-5:
                k0 = k

def test_poisson(self):
        import meshzoo
        from scipy.sparse import linalg

        # Define the problem
        class Poisson(LinearFvmProblem):
            @staticmethod
            def apply(u):
                return integrate(lambda x: - n_dot_grad(u(x)), dS) \
                       - integrate(lambda x: 1.0, dV)
            dirichlet = [(lambda x: 0.0, ['Boundary'])]

        # Create mesh using meshzoo
        vertices, cells = meshzoo.rectangle.create_mesh(
                0.0, 1.0, 0.0, 1.0,
                21, 21,
                zigzag=True
                )
        mesh = pyfvm.meshTri.meshTri(vertices, cells)

        linear_system = pyfvm.discretize(Poisson, mesh)

        x = linalg.spsolve(linear_system.matrix, linear_system.rhs)

        k0 = -1
        for k, coord in enumerate(mesh.node_coords):
            # print(coord - [0.5, 0.5, 0.0])
            if numpy.linalg.norm(coord - [0.5, 0.5, 0.0]) < 1.0e-5:
                k0 = k
                break

edge_ce_ratio = mesh.ce_ratios[edge_ids]

        # project the magnetic potential on the edge at the midpoint
        magnetic_potential = \
            0.5 * numpy.cross(self.magnetic_field, edge_midpoint.T).T

        # The dot product , executed for many points
        # at once; cf. .
        beta = numpy.einsum('ij, ij->i', magnetic_potential.T, edge.T)

        return numpy.array([
            [edge_ce_ratio, -edge_ce_ratio * numpy.exp(1j * beta)],
            [-edge_ce_ratio * numpy.exp(-1j * beta), edge_ce_ratio]
            ])

vertices, cells = meshzoo.rectangle.create_mesh(0.0, 1.0, 0.0, 1.0, 31, 31)
mesh = pyfvm.mesh_tri.MeshTri(vertices, cells)

# Equivalently, one could have written
#
# from pyfvm.form_language import integrate, dV
# class GinzburgLandau(object):
#     def apply(self, psi):
#         return Energy() \
#             + integrate(lambda x: psi(x) * (V + g * abs(psi(x))**2), dV)
# f, _ = pyfvm.discretize(GinzburgLandau(), mesh)
#
# The Jacobian still has to be specified manually because of its special
# structure.

keo = pyfvm.get_fvm_matrix(mesh, [Energy()], [], [], [])

class D1(Subdomain):
    def is_inside(self, x): return x[1] < 0.5
    is_boundary_only = True


class Poisson(object):
    def apply(self, u):
        return integrate(lambda x: -n_dot_grad(u(x)), dS) \
                + integrate(lambda x: 3.0, dGamma) \
                - integrate(lambda x: 1.0, dV)

    def dirichlet(self, u):
        return [(u, D1())]

vertices, cells = meshzoo.rectangle.create_mesh(0.0, 1.0, 0.0, 1.0, 51, 51)
mesh = pyfvm.mesh_tri.MeshTri(vertices, cells)

matrix, rhs = pyfvm.discretize_linear(Poisson(), mesh)

u = linalg.spsolve(matrix, rhs)

mesh.write('out.vtu', point_data={'u': u})