function mesh = meshFlatStrain( mesh, scale )
%mesh = meshFlatStrain( mesh, scale )
%   Set the residual strain to be the strain the mesh must have, for its
%   relaxed state to be flat.  If scale is supplied, the strain is
%   multiplied by that amount

    if nargin < 2
        scale = 1;
    end

    mesh = rectifyVerticals( mesh );
    mesh = computeGNGlobal( mesh );
    numCells = size(mesh.tricellvxs,1);
    MAXFLATSTRAIN = 0.1;
    maxallrs = 0;
    for ci=1:numCells
        trivxs = mesh.tricellvxs(ci,:);
        prismvxs = [ trivxs*2-1, trivxs*2 ];
        cellvxCoords = mesh.prismnodes( prismvxs, : );
        d = cellFlatDisplacements( cellvxCoords, mesh.unitcellnormals(ci,:) );
        mesh.celldata(ci).eps0Gauss = zeros( 6, 6 );
        mesh.celldata(ci) = computeDisplacementStrains( mesh.celldata(ci), d );
        mesh.celldata(ci).residualStrain = mesh.celldata(ci).displacementStrain * scale;
        maxrs = max( abs(mesh.celldata(ci).residualStrain(:)) );
        if maxrs > 100000
          % x = 1; % Uncomment this command to allow a breakpoint to be set here.
        end
        if false && (maxrs > MAXFLATSTRAIN)
            strainScale = MAXFLATSTRAIN/maxrs;
            mesh.celldata(ci).residualStrain = strainScale * mesh.celldata(ci).residualStrain;
        end
        maxallrs = max( maxrs, maxallrs );
        mesh.celldata(ci).displacementStrain = zeros( 6, 6 );
        mesh.celldata(ci).residualStress = mesh.globalProps.D * mesh.celldata(ci).residualStrain;
    end
  % meshFlatStrain_maxallrs = maxallrs
    if true && (maxallrs > MAXFLATSTRAIN)
        strainScale = MAXFLATSTRAIN/maxallrs;
        for ci=1:numCells
            mesh.celldata(ci).residualStrain = strainScale * mesh.celldata(ci).residualStrain;
            mesh.celldata(ci).residualStress = strainScale * mesh.celldata(ci).residualStress;
        end
    end
end
