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

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

Switch thermalStress(thermalProperties.lookup("thermalStress"));

volScalarField threeKalpha
(
    IOobject
    (
        "threeKalpha",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("0", dimensionSet(0, 2, -2 , -1, 0), 0.0)
);


volScalarField DT
(
    IOobject
    (
        "DT",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("0", dimensionSet(0, 2, -1 , 0, 0), 0.0)
);


if (thermalStress)
{

    autoPtr<volScalarField> CPtr;

    IOobject CIO
    (
        "C",
        runTime.timeName(0),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    );

    const dictionary& CDict(thermalProperties.subDict("C"));
    word CType(CDict.lookup("type"));
    if (CType == "uniform")
    {
        scalar CValue(readScalar(CDict.lookup("value")));

        CPtr.reset
        (
            new volScalarField
            (
                CIO,
                mesh,
                dimensionedScalar
                (
                    "C",
                    dimensionSet(0, 2, -2 , -1, 0),
                    CValue
                )
            )
        );

    }
    else if (CType == "field")
    {
        CIO.readOpt() = IOobject::MUST_READ;

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

    volScalarField& C = CPtr();

    autoPtr<volScalarField> rhoKPtr;

    IOobject rhoKIO
    (
        "k",
        runTime.timeName(0),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    );

    const dictionary& kDict(thermalProperties.subDict("k"));
    word kType(kDict.lookup("type"));
    if (kType == "uniform")
    {
        scalar rhoKValue(readScalar(kDict.lookup("value")));

        rhoKPtr.reset
        (
            new volScalarField
            (
                rhoKIO,
                mesh,
                dimensionedScalar
                (
                    "rhoK",
                    dimensionSet(1, 1, -3 , -1, 0),
                    rhoKValue
                )
            )
        );

    }
    else if (kType == "field")
    {
        rhoKIO.readOpt() = IOobject::MUST_READ;

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

    volScalarField& rhoK = rhoKPtr();

    autoPtr<volScalarField> alphaPtr;

    IOobject alphaIO
    (
        "alpha",
        runTime.timeName(0),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    );


    const dictionary& alphaDict(thermalProperties.subDict("alpha"));
    word alphaType(alphaDict.lookup("type"));

    if (alphaType == "uniform")
    {
        scalar alphaValue(readScalar(alphaDict.lookup("value")));
        alphaPtr.reset
        (
            new volScalarField
            (
                alphaIO,
                mesh,
                dimensionedScalar
                (
                    "alpha",
                    inv(dimTemperature),
                    alphaValue
                )
            )
        );
    }
    else if (alphaType == "field")
    {
        alphaIO.readOpt() = IOobject::MUST_READ;

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

    volScalarField& alpha = alphaPtr();

    Info<< "Normalising k : k/rho\n" << endl;
    volScalarField k(rhoK/rho);

    Info<< "Calculating thermal coefficients\n" << endl;

    threeKalpha = threeK*alpha;
    DT = k/C;
}