function m = makeCellGrid( m, spacing )
%m = makeCellGrid( m, spacing )
%   Create a cellular layer in the form of a regular square grid.  The mesh
%   m is assumed to be flat or nearly so.

    xmin = min(m.nodes(:,1));
    xmax = max(m.nodes(:,1));
    xdiam = xmax-xmin;
    ymin = min(m.nodes(:,2));
    ymax = max(m.nodes(:,2));
    ydiam = ymax-ymin;
    nxcells = ceil( xdiam/spacing );
    nxpts = nxcells+1;
    xmargin = (nxcells*spacing - xdiam)/2;
    nycells = ceil( ydiam/spacing );
    nypts = nycells+1;
    ymargin = (nycells*spacing - ydiam)/2;
    npts = nxpts * nypts;
    vxpos = (0:nxcells)'*spacing + xmin - xmargin;
    vypos = (0:nycells)*spacing + ymin - ymargin;
    vpos = [ repmat( vxpos, nypts, 1 ), ...
             reshape( repmat( vypos, nxpts, 1 ), [], 1 ), ...
             zeros( npts, 1 ) ];
    ci = zeros( npts, 1 );
    bc = zeros( npts, 3 );
    bcerr = zeros( npts, 1 );
    abserr = zeros( npts, 1 );
    [ ci(1), bc(1,:), bcerr(1), abserr(1) ] = findFEMCell( m, vpos(1,:) );
    for i=2:size(vpos,1)
        [ ci(i), bc(i,:), bcerr(i), abserr(i) ] = findFEMCell( m, vpos(i,:), ci(i-1) );
        if bcerr(i) < 0
            [ ci(i), bc(i,:), bcerr(i), abserr(i) ] = findFEMCell( m, vpos(i,:) );
        end
      % fprintf( 1, 'FE %d bc [ %.2f %.2f %.2f ] bce %.2f abse %.2f\n', ...
      %     ci(i), bc(i,:), bcerr(i), abserr(i) );
    end
    a = (1:nxcells)';
    b = 0:nxpts:(npts-nxpts*2);
    topleft = reshape( repmat(a,1,nycells) + repmat(b,nxcells,1), [], 1 );
    topright = topleft+1;
    bottomleft = topleft+nxpts;
    bottomright = bottomleft+1;
    cells = [ topleft, topright, bottomright, bottomleft ];
    okpts = bcerr >= -0.01*spacing;
    okcells = all( okpts(cells), 2 );
    
    newToOldPt = find( okpts );
    nokpts = length(newToOldPt);
    oldToNewPt = zeros(npts,1);
    oldToNewPt(okpts) = 1:nokpts;
    cells = oldToNewPt( cells( okcells, : ) );
    nokcells = size(cells,1);
    vpos = vpos( okpts, : );
    bc = bc( okpts, : );
    ci = ci( okpts );

    m.secondlayer = newemptysecondlayer();
    for i=1:nokcells
        m.secondlayer.cells(i).vxs = cells(i,:);
    end
    m.secondlayer.vxFEMcell = ci;
    m.secondlayer.vxBaryCoords = bc;
    m.secondlayer.cell3dcoords = vpos;
    m.secondlayer.side = true( length(m.secondlayer.cells), 1 );
    m = completesecondlayer( m );
end
