function [m,ok] = onestep( m, useGrowthTensors )
%[m,ok] = onestep( m, useGrowthTensors )
%   Run one step of the iteration:
%   * Calculate the polarisation gradients.
%   * Perform a FEM calculation.
%   * Restore flatness/thickness.
%   * Perform dilution-by-growth.
%   * Generate the FEM cell normal vectors.
%   * Split large FEM cells.
%   * Recalculate FEM cell areas.
%   * Recalculate the second layer global positions and split large
%     second layer cells.
%   The result OK will be false if the user interrupts in the middle of a
%   simulation step or there is an error in the interaction function.

    global gOLD_THICKNESS gNEW_K_NOR
    
    ok = true;
    if nargin < 2, useGrowthTensors = 0; end
    m = calcPolGrad( m );
  % maxgrad1 = max(abs(m.gradpolgrowth(:,1)))
    m = makeCellFrames( m );
    m.globalDynamicProps.previousArea = m.globalDynamicProps.currentArea;
    if m.globalProps.flatten
        m = meshFlatStrain( m, m.globalProps.flattenratio );
    end
    if m.globalProps.growthEnabled
        transportspeed = cell( size( m.transportfield ) );
        for i=1:length(transportspeed)
            tv = transportvectors( m, i );
            if ~isempty( tv )
                transportspeed{i} = sqrt(sum(tv.^2,2));
            end
        end
        switch m.globalProps.thicknessMode
            case { 'direct', 'specified' }
                oldthicknesses = getThickness( m );
                KNOR_thicknesses = oldthicknesses .* ...
                    (1 + m.morphogens(:,m.mgenNameToIndex.KNOR) * m.globalProps.timestep);
        end
        if m.globalProps.plasticGrowth
            [m,u,G] = growthDisplacements( m );
            u = reshape( [ u, u ]', 3, [] )';
            m.effectiveGrowthTensor = G;
        else
            [m,u] = totalK( m, useGrowthTensors );
            ok = ~isempty(u);
            m.effectiveGrowthTensor = zeros( size(m.tricellvxs,1), ...
                                             size(m.effectiveGrowthTensor,2) );
        end
        % If the computation failed, u may have been returned as empty.
        if m.versioninfo.mgenversion==0
            THICKNESSMGEN = gOLD_THICKNESS;
        else
            THICKNESSMGEN = gNEW_K_NOR;
        end
        moved = (~isempty(u)) ...
                && (any(u(:)) ...
                    || (strcmp( m.globalProps.thicknessMode, 'direct' ) ...
                        && any(m.morphogens(THICKNESSMGEN,:))));
        if moved
            m.prismnodes = m.prismnodes + u;
            switch m.globalProps.thicknessMode
                case 'actual'
                    % Nothing.
                case { 'direct', 'specified' }
                    [m,displacements] = setThickness( m, KNOR_thicknesses );
                  % sd = sum(displacements,2);
                  % su = sum(u,2);
                  % r = zeros(size(sd));
                  % r(su~=0) = sd(su~=0)./su(su~=0);
                  % specerrmax = max(r)
                    u = u + displacements;
                case 'scaled'
                    [m,u] = restorethickness( m, u );
            end
            m.globalProps.trinodesvalid = 0;
            m = makeTRIvalid( m );
            if m.globalProps.alwaysFlat
                [m,u] = restoreflatness( m, u );
            else
                m = makeAreasAndNormals( m );
                if (~strcmp( m.globalProps.thicknessMode, 'actual' )) || m.globalProps.rectifyverticals
                    m = rectifyVerticals( m );
                end
                m = makebendangles( m );
            end
            if hasSecondLayer( m )
                m = calcCloneVxCoords( m );
            end
            % Rescale transport vectors.
            for i=1:length(transportspeed)
                if ~isempty( tv )
                    newtv = transportvectors( m, mi );
                    newtransportspeed = sqrt(sum(newtv.^2,2));
                    nonzeros = transportspeed{i} > 0;
                    scaling = ones(length(transportspeed{i}),1);
                    scaling(nonzeros) = newtransportspeed(nonzeros)./transportspeed{i}(nonzeros);
                    m.transportfield{i} = m.transportfield{i} .* repmat( scaling, 1, 3 );
                end
            end
            if ~m.globalProps.flatten
                m = dilateSubstances( m, u );
                m = trysplit( m );
                if m.globalProps.allowFlipEdges
                    m = flipedges( m );
                end
                if m.globalProps.allowElideEdges
                    m = tryElideEdge( m, m.globalProps.mincellangle );
                end
                if hasSecondLayer( m )
                    if m.globalProps.allowSplitBio
                        m = splitSecondLayer( m );
                    end
                    m = perturbSecondLayer( m, ...
                            m.secondlayer.jiggleAmount * m.secondlayer.splitThreshold );
                    m = calcBioACellAreas( m );
                    if strcmp( m.globalProps.biocolormode, 'area' )
                        m = setSecondLayerColorsByArea( m );
                    end
                end
            end
            m = calcmeshareas( m );        
        end
    end
    
    m = calcPolGrad( m );
  % maxgrad2 = max(abs(m.gradpolgrowth(:,1)))

    if m.globalProps.growthEnabled && (m.globalProps.boingNeeded == 1)
        m.globalProps.boingNeeded = 2;
    end
    m.saved = 0;
end
    
