{
    scalar t = runTime.value();
    scalar tNext = t + runTime.deltaT().value();

    instantList Times = runTime.times();

    sigma = dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero);
    DSigmaCorr =
        dimensionedSymmTensor
        (
            "zero",
            dimForce/dimArea,
            symmTensor::zero
        );

    for (label i=1; i<Times.size(); i++)
    {
        runTime.setTime(Times[i], i);

        if(runTime.timeIndex() != i)
        {
            FatalErrorIn(args.executable())
                << "The strain increment DEpsilon must be stored for "
                << "each calculated time step. "
                << exit(FatalError);
        }

        IOobject DEpsilonHeader
        (
            "DEpsilon",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ
        );

        // Check DEpsilon exists
        if (DEpsilonHeader.headerOk())
        {
            volSymmTensorField DEpsilonOld(DEpsilonHeader, mesh);

            scalar tau = runTime.value() - m*runTime.deltaT().value();

            if(tau < 0)
            {
                sigma += 2.0*rheology.mu(t)*DEpsilonOld
                    + rheology.lambda(t)*(I*tr(DEpsilonOld));

                DSigmaCorr += 2.0*rheology.mu(tNext)*DEpsilonOld
                    + rheology.lambda(tNext)*(I*tr(DEpsilonOld));
            }
            else
            {
                sigma += 2.0*rheology.mu(t - tau)*DEpsilonOld
                    + rheology.lambda(t - tau)*(I*tr(DEpsilonOld));

                DSigmaCorr += 2.0*rheology.mu(tNext - tau)*DEpsilonOld
                    + rheology.lambda(tNext - tau)*(I*tr(DEpsilonOld));
            }
        }
        else
        {
            Info << "No DEpsilon" << endl;
        }
    }

    if(Times.size()>=2)
    {
        runTime++;
    }

    scalar tau = runTime.value() - m*runTime.deltaT().value();

    sigma += 2.0*rheology.mu(t - tau)*DEpsilon
        + rheology.lambda(t - tau)*(I*tr(DEpsilon));

    DSigmaCorr += 2.0*rheology.mu(tNext - tau)*DEpsilon
        + rheology.lambda(tNext - tau)*(I*tr(DEpsilon));

    DSigmaCorr -= sigma;
}