def vis_splitting(Verts, splitting, output='vtk', fname='output.vtu'):
    """Coarse grid visualization for C/F splittings.

    Parameters
    ----------
    Verts : {array}
        coordinate array (N x D)
    splitting : {array}
        coarse(1)/fine(0) flags
    fname : {string, file object}
        file to be written, e.g. 'output.vtu'
    output : {string}
        'vtk' or 'matplotlib'

    Returns
    -------
        - Displays in screen or writes data to .vtu file for use in paraview
          (xml 0.1 format)

    Notes
    -----
    D :
        dimension of coordinate space
    N :
        # of vertices in the mesh represented in Verts
    Ndof :
        # of dof (= ldof * N)

        - simply color different points with different colors.  This works
          best with classical AMG.

        - writes a file (or opens a window) for each dof

        - for Ndof>1, they are assumed orderd [...dof1..., ...dof2..., etc]

    Examples
    --------
    >>> import numpy as np
    >>> from pyamg.vis.vis_coarse import vis_splitting
    >>> Verts = np.array([[0.0,0.0],
    ...                   [1.0,0.0],
    ...                   [0.0,1.0],
    ...                   [1.0,1.0]])
    >>> splitting = np.array([0,1,0,1,1,0,1,0])    # two variables
    >>> vis_splitting(Verts,splitting,output='vtk',fname='output.vtu')

    >>> from pyamg.classical import RS
    >>> from pyamg.vis.vis_coarse import vis_splitting
    >>> from pyamg.gallery import load_example
    >>> data = load_example('unit_square')
    >>> A = data['A'].tocsr()
    >>> V = data['vertices']
    >>> E2V = data['elements']
    >>> splitting = RS(A)
    >>> vis_splitting(Verts=V,splitting=splitting,output='vtk',
                      fname='output.vtu')

    """
    check_input(Verts, splitting)

    N = Verts.shape[0]
    Ndof = int(len(splitting) / N)

    E2V = np.arange(0, N, dtype=int)

    # adjust name in case of multiple variables
    a = fname.split('.')
    if len(a) < 2:
        fname1 = a[0]
        fname2 = '.vtu'
    elif len(a) >= 2:
        fname1 = "".join(a[:-1])
        fname2 = a[-1]
    else:
        raise ValueError('problem with fname')

    new_fname = fname
    for d in range(0, Ndof):
        # for each variables, write a file or open a figure

        if Ndof > 1:
            new_fname = fname1 + '_%d.' % (d+1) + fname2

        cdata = splitting[(d*N):((d+1)*N)]

        if output == 'vtk':
            write_basic_mesh(Verts=Verts, E2V=E2V, mesh_type='vertex',
                             cdata=cdata, fname=new_fname)
        elif output == 'matplotlib':
            from pylab import figure, show, plot, xlabel, ylabel, title, axis
            cdataF = np.where(cdata == 0)[0]
            cdataC = np.where(cdata == 1)[0]
            xC = Verts[cdataC, 0]
            yC = Verts[cdataC, 1]
            xF = Verts[cdataF, 0]
            yF = Verts[cdataF, 1]
            figure()
            plot(xC, yC, 'r.', xF, yF, 'b.', clip_on=True)
            title('C/F splitting (red=coarse, blue=fine)')
            xlabel('x')
            ylabel('y')
            axis('off')
            show()
        else:
            raise ValueError('problem with outputtype')