function m = renumberMesh( m, rnodes, redges, rcells, ...
                              rnodesmap, redgesmap, rcellsmap )
% m has had vertexes, edges, and cells merged or deleted.
% rnodes is an array whose length is the number of vertexes of the old
% version of m.  If vi is the index of such a vertex, then rnodes(vi) is
% the index of the same vertex in the current mesh.  If vi has been
% deleted, rnodes(vi) is zero.
% Similarly for redges and rcells.  The ordering of indexes is assumed to
% be unchanged.
% rnodesmap, redgesmap, rcellsmap are boolean maps of the same information
% in rnodes, redges, rcells.
% At least one of rnodes or rnodesmap must be supplied; if only one is, the
% other will be computed from it.  The same applies to redges and
% redgesmap, and rcells and rcellsmap.

    if nargin < 5
        rnodesmap = [];
    end
    if nargin < 6
        redgesmap = [];
    end
    if nargin < 7
        rcellsmap = [];
    end
    [rnodes,rnodesmap] = makeindexmap( rnodes, rnodesmap );
    [redges,redgesmap] = makeindexmap( redges, redgesmap );
    [rcells,rcellsmap] = makeindexmap( rcells, rcellsmap );

    % Use the maps to eliminate items; use the renumbering arrays to
    % transform values.

    % Per-node data.
    if ~isempty(rnodesmap)
        m.nodes = m.nodes(rnodesmap,:);
        m.morphogens = m.morphogens(rnodesmap,:);
        m.mgen_production = m.mgen_production(rnodesmap,:);
        m.morphogenclamp = m.morphogenclamp(rnodesmap,:);
        m.fixedDFmap = m.fixedDFmap(rnodesmap,:);
        m.nodecelledges = { m.nodecelledges{rnodesmap} };
        if ~isempty(m.growthanglepervertex)
            m.growthanglepervertex = m.growthanglepervertex(rnodesmap);
        end
      % if false && (length(m.globalProps.stitchDFs) > 0)
      %     for i=1:length(m.globalProps.stitchDFs)
      %         dfs = m.globalProps.stitchDFs{i};
      %     end
      % end

        if isfield( m, 'growthTensorPerVertex' )
            m.growthTensorPerVertex = ...
                m.growthTensorPerVertex( rnodesmap, : );
        end
        % Delete the locator node if necessary.
        if (m.globalDynamicProps.locatenode > 0) && rnodesmap(m.globalDynamicProps.locatenode)
            m.globalDynamicProps.locatenode = 0;
            m.globalDynamicProps.locateDFs = [0 0 0];
        end

        % Per-prism-node data.
        pnodesmap = reshape( [ rnodesmap(:), rnodesmap(:) ]', [], 1 );
        m.prismnodes = m.prismnodes(pnodesmap,:);
        if ~isempty(m.displacements)
            m.displacements = m.displacements(pnodesmap,:);
        end
    end
    if (~isempty(redges)) || (~isempty(rcells))
        for i=1:length(m.nodecelledges)
            nce = m.nodecelledges{i};
            if isempty(redges)
                e = nce(1,:);
            else
                e = redges(nce(1,:));
                e = e(e~=0);
            end
            c = nce(2,e~=0);
            if ~isempty(rcells)
                c(c ~= 0) = rcells(c(c ~= 0));
            end
            zi = find(c==0,1);
            % Need to rotate to put any zero cell at the end.  Perhaps complain
            % if more than one zero.
            if ~isempty(zi)
                recycle = [(zi+1):length(e), 1:zi];
                e = e(recycle);
                c = c(recycle);
            end
            m.nodecelledges{i} = [ e; c ];
        end
    end
    
    % Per-edge data.
        
    if ~isempty(redgesmap)
        m.edgecells = m.edgecells(redgesmap,:);
    end
    if ~isempty(rcells)
        m.edgecells(m.edgecells~=0) = rcells(m.edgecells(m.edgecells~=0));
    end
    badedges = m.edgecells(:,1)==0;
    m.edgecells(badedges,:) = m.edgecells(badedges,[2 1]);
    if ~isempty(redgesmap)
        m.currentbendangle = m.currentbendangle(redgesmap);
        m.initialbendangle = m.initialbendangle(redgesmap);
        m.seams = m.seams(redgesmap);
    end
    if ~isempty(redgesmap)
        m.edgeends = m.edgeends(redgesmap,:);
    end
    if ~isempty(rnodes)
        m.edgeends = rnodes(m.edgeends);
    end
    
    % Per-cell data.
    if ~isempty(rcellsmap)
        m.tricellvxs = m.tricellvxs(rcellsmap,:);
        m.celledges = m.celledges(rcellsmap,:);
        m.celldata = m.celldata(rcellsmap);
        if ~isempty( m.cellFrames )
            m.cellFrames = m.cellFrames(:,:,rcellsmap);
        end
        m.cellbulkmodulus = m.cellbulkmodulus(rcellsmap);
        m.cellpoisson = m.cellpoisson(rcellsmap);
        m.cellstiffness = m.cellstiffness(:,:,rcellsmap);
        m.effectiveGrowthTensor = m.effectiveGrowthTensor(rcellsmap,:);
        m.unitcellnormals = m.unitcellnormals(rcellsmap,:);
        m.gradpolgrowth = m.gradpolgrowth(rcellsmap,:);
        m.polfreeze = m.polfreeze(rcellsmap,:);
        m.polfreezebc = m.polfreezebc(rcellsmap,:);
        m.polfrozen = m.polfrozen(rcellsmap);
        m.cellareas = m.cellareas(rcellsmap);
        if ~isempty(m.decorBCs)
            m.decorBCs = m.decorBCs( rcellsmap( m.decorFEs ), : );
        end
        if ~isempty(m.decorFEs)
            m.decorFEs = rcells( m.decorFEs( rcellsmap( m.decorFEs ) ) );
        end
        if ~isempty(m.growthangleperFE)
            m.growthangleperFE = m.growthangleperFE(rcellsmap);
        end
        if isfield( m, 'celllabel' )
            m.celllabel = m.celllabel(rcellsmap);
        end
    end
    if ~isempty(rnodes)
        m.tricellvxs = rnodes(m.tricellvxs);
    end
    if ~isempty(redges)
        m.celledges = redges(m.celledges);
    end

    if hasSecondLayer( m ) && ~isempty(rcells)
        m.secondlayer.vxFEMcell = rcells(m.secondlayer.vxFEMcell);
    end
end

function [reindex,bitmap] = makeindexmap( reindex, bitmap )
    if isempty(bitmap)
        bitmap = reindex ~= 0;
    elseif isempty(reindex)
        reindex = makereindex(bitmap);
    end
end

function ri = makereindex( bitmap )
    x = find(bitmap);
    ri = zeros(1,length(bitmap),'int32');
    ri(x) = 1:length(x);
end

