def rectangular_plane_layout(mesh,corner, closed = False,I=1.):
    """
    corner: sorted list of four corners (x,y,z)

    2--3
    |  |
    1--4

    y
    |
    |--> x

    Output:
    Js

    """

    Jx = np.zeros(mesh.nEx)
    Jy = np.zeros(mesh.nEy)
    Jz = np.zeros(mesh.nEz)

    indy1 = np.logical_and( \
            np.logical_and( \
            np.logical_and(mesh.gridEy[:,0]>=corner[0,0],mesh.gridEy[:,0]<=corner[1,0]), \
            np.logical_and(mesh.gridEy[:,1] >=corner[0,1] , mesh.gridEy[:,1]<=corner[1,1] )),
            (mesh.gridEy[:,2] == corner[0,2]
            )
            )
    indx1 = np.logical_and( \
            np.logical_and( \
            np.logical_and(mesh.gridEx[:,0]>=corner[1,0],mesh.gridEx[:,0]<=corner[2,0]), \
            np.logical_and(mesh.gridEx[:,1] >=corner[1,1] , mesh.gridEx[:,1]<=corner[2,1] )),
            (mesh.gridEx[:,2] == corner[1,2]
            )
            )
    indy2 = np.logical_and( \
            np.logical_and( \
            np.logical_and(mesh.gridEy[:,0]>=corner[2,0],mesh.gridEy[:,0]<=corner[3,0]), \
            np.logical_and(mesh.gridEy[:,1] <=corner[2,1] , mesh.gridEy[:,1]>=corner[3,1] )),
            (mesh.gridEy[:,2] == corner[2,2]
            )
            )

    if closed:
        indx2 = np.logical_and( \
                np.logical_and( \
                np.logical_and(mesh.gridEx[:,0]>=corner[0,0],mesh.gridEx[:,0]<=corner[3,0]), \
                np.logical_and(mesh.gridEx[:,1] >=corner[0,1] , mesh.gridEx[:,1]<=corner[3,1] )),
                (mesh.gridEx[:,2] == corner[0,2]
                )
                )

    else:
        indx2 = []

    Jy[indy1] = -I
    Jx[indx1] = -I
    Jy[indy2] = I
    Jx[indx2] = I

    J = np.hstack((Jx,Jy,Jz))
    J = J*mesh.edge

    return J