if (runTime.outputTime())
{
    volScalarField epsilonEq
    (
        IOobject
        (
            "epsilonEq",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        sqrt((2.0/3.0)*magSqr(dev(epsilon)))
    );

    Info<< "Max epsilonEq = " << max(epsilonEq).value() << endl;

    volScalarField sigmaEq
    (
        IOobject
        (
            "sigmaEq",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        sqrt((3.0/2.0)*magSqr(dev(sigma)))
    );

    Info<< "Max sigmaEq = " << max(sigmaEq).value() << endl;

    pointMesh pMesh(mesh);
    pointScalarField contactPointGap
    (
        IOobject
        (
            "contactPointGap",
            runTime.timeName(),
            mesh,
           IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        pMesh,
        dimensionedScalar("zero", dimless, 0.0)
    );

    forAll(mesh.boundary(), patchi)
    {
        if(DU.boundaryField()[patchi].type() == solidContactFvPatchVectorField::typeName)
        {
            const solidContactFvPatchVectorField& DUpatch =
                refCast<const solidContactFvPatchVectorField>
                    (DU.boundaryField()[patchi]);

            if(!DUpatch.master())
            {
                const labelList& meshPoints = mesh.boundaryMesh()[patchi].meshPoints();
                const scalarField gap = DUpatch.normalContactModelPtr()->slaveContactPointGap();
                forAll(meshPoints, pointi)
                {
                    contactPointGap[meshPoints[pointi]] = gap[pointi];
                }
            }
        }
    }


    //- total force
/*
    forAll(mesh.boundary(), patchi)
    {
        Info << "Patch " << mesh.boundary()[patchi].name() << endl;
        vectorField totalForce = mesh.Sf().boundaryField()[patchi] & sigma.boundaryField()[patchi];

        vector force = sum( totalForce );
        Info << "\ttotal force is " << force << " N" << endl;

        tensorField F = I + gradDU.boundaryField()[patchi];
        tensorField Finv = inv(F);
        //vectorField nCurrent = Finv & n.boundaryField()[patchi];
        //nCurrent /= mag(nCurrent);
        //scalar normalForce = sum( nCurrent & totalForce );
        scalar normalForce = sum( n.boundaryField()[patchi] & totalForce );
        Info << "\tnormal force is " << normalForce << " N" << endl;
        //scalar shearForce = mag(sum( (I - sqr(nCurrent)) & totalForce ));
        scalar shearForce = mag(sum( (I - sqr(n.boundaryField()[patchi])) & totalForce ));
        Info << "\tshear force is " << shearForce << " N" << endl;

        // if(mesh.boundary()[patchi].type() != "empty")
        //   {
    //     vector Sf0 = Sf.boundaryField()[patchi][0];
    //     symmTensor sigma0 = sigma.boundaryField()[patchi][0];
    //     Info << "sigmab[0] is " << sigma0 << nl
    //      << "Sfb is " << Sf0 << nl
    //      << "force is " << (Sf.boundaryField()[patchi][0]&sigma.boundaryField()[patchi][0]) << nl
    //      << "Sfx*sigmaxx " << (Sf0[vector::X]*sigma0[symmTensor::XX]) << nl
    //      << "Sfy*sigmaxy " << (Sf0[vector::Y]*sigma0[symmTensor::XY]) << nl
    //      << "Sfx*sigmayx " << (Sf0[vector::X]*sigma0[symmTensor::XY]) << nl
    //      << "Sfy*sigmayy " << (Sf0[vector::Y]*sigma0[symmTensor::YY]) << nl
    //      << endl;
        //vector SfTL(-0.000137451, 0.00383599, -4.76878e-20);
        // vector SfTL = Finv[0] & vector(0,0.004,0);
        // Info << "SfTLx*sigmaxx " << (SfTL[vector::X]*sigma0[symmTensor::XX]) << nl
        //      << "SfTLy*sigmaxy " << (SfTL[vector::Y]*sigma0[symmTensor::XY]) << nl
        //      << "SfTLx*sigmayx " << (SfTL[vector::X]*sigma0[symmTensor::XY]) << nl
        //      << "SfTLy*sigmayy " << (SfTL[vector::Y]*sigma0[symmTensor::YY]) << nl
        //      << endl;
    //      }
    }
*/

    runTime.write();
  }