Info<< "Reading mechanical properties\n" << endl;

    IOdictionary mechanicalProperties
    (
        IOobject
        (
            "mechanicalProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );

    const dictionary& rhoDict(mechanicalProperties.subDict("rho"));
    word rhoType(rhoDict.lookup("type"));

    autoPtr<volScalarField> rhoPtr;

    IOobject rhoIO
    (
        "rho",
        runTime.timeName(0),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    );

    if (rhoType == "uniform")
    {
        scalar rhoValue(readScalar(rhoDict.lookup("value")));

        rhoPtr.reset
        (
            new volScalarField
            (
                rhoIO,
                mesh,
                dimensionedScalar
                (
                    "rho",
                    dimMass/dimVolume,
                    rhoValue
                )
            )
        );
    }
    else if (rhoType == "field")
    {
        rhoIO.readOpt() = IOobject::MUST_READ;

        rhoPtr.reset
        (
            new volScalarField
            (
                rhoIO,
                mesh
            )
        );
    }
    else
    {
        FatalErrorInFunction
            << "Valid type entries are uniform or field for rho"
            << abort(FatalError);
    }

    volScalarField& rho = rhoPtr();

    const dictionary& EDict(mechanicalProperties.subDict("E"));
    word EType(EDict.lookup("type"));

    autoPtr<volScalarField> EPtr;

    IOobject EHeader
    (
        "E",
        runTime.timeName(0),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    );

    if (EType == "uniform")
    {
        scalar rhoEValue(readScalar(EDict.lookup("value")));

        EPtr.reset
        (
            new volScalarField
            (
                EHeader,
                mesh,
                dimensionedScalar
                (
                    "Erho",
                    dimMass/dimLength/sqr(dimTime),
                    rhoEValue
                )
            )
        );
    }
    else if (EType == "field")
    {
        EHeader.readOpt() = IOobject::MUST_READ;

        EPtr.reset
        (
            new volScalarField
            (
                EHeader,
                mesh
            )
        );
    }
    else
    {
        FatalErrorInFunction
            << "Valid type entries are uniform or field for E"
            << abort(FatalError);
    }

    volScalarField& rhoE = EPtr();

    autoPtr<volScalarField> nuPtr;

    IOobject nuIO
    (
        "nu",
        runTime.timeName(0),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    );

    const dictionary& nuDict(mechanicalProperties.subDict("nu"));
    word nuType(nuDict.lookup("type"));

    if (nuType == "uniform")
    {
        scalar nuValue(readScalar(nuDict.lookup("value")));
        nuPtr.reset
        (
            new volScalarField
            (
                nuIO,
                mesh,
                dimensionedScalar
                (
                    "nu",
                    dimless,
                    nuValue
                )
            )
        );
    }
    else if (nuType == "field")
    {
        nuIO.readOpt() = IOobject::MUST_READ;
        nuPtr.reset
        (
            new volScalarField
            (
                nuIO,
                mesh
            )
        );
    }
    else
    {
        FatalErrorInFunction
            << "Valid type entries are uniform or field for nu"
            << abort(FatalError);
    }

    volScalarField& nu = nuPtr();

    Info<< "Normalising E : E/rho\n" << endl;
    volScalarField E(rhoE/rho);

    Info<< "Calculating Lame's coefficients\n" << endl;

    volScalarField mu(E/(2.0*(1.0 + nu)));
    volScalarField lambda(nu*E/((1.0 + nu)*(1.0 - 2.0*nu)));
    volScalarField threeK(E/(1.0 - 2.0*nu));

    Switch planeStress(mechanicalProperties.lookup("planeStress"));

    if (planeStress)
    {
        Info<< "Plane Stress\n" << endl;

        lambda = nu*E/((1.0 + nu)*(1.0 - nu));
        threeK = E/(1.0 - nu);
    }
    else
    {
        Info<< "Plane Strain\n" << endl;
    }