calcViscousDiffusionTerm
{
    type            coded;
    libs            ("libutilityFunctionObjects.so");
    enabled         true;
    log             true;
    timeStart       $startTimeFieldAverage;
    timeEnd         $endTime;
    executeControl  adjustableRunTime;
    executeInterval $writeInterval;
    writeControl    adjustableRunTime;
    writeInterval   $writeInterval;

    codeOptions     #{ -I$(LIB_SRC)/OpenFOAM/lnInclude -I$(LIB_SRC)/finiteVolume/lnInclude #};

    codeInclude     #{ #include "wallFvPatch.H" #};

    // --- Inputs ---
    RMeanName           UPrime2Mean;
    transportProperties transportProperties;

    // --- Output ---
    OutputName          D_ij_RSTE;

    name                calcDij;

    codeWrite
    #{
        #include "fvCFD.H"

        const word RMeanNameIn      = this->dict().lookupOrDefault<word>("RMeanName", "UPrime2Mean");
        const word transportPropsNm = this->dict().lookupOrDefault<word>("transportProperties", "transportProperties");
        const word viscDiffOutNm    = this->dict().lookupOrDefault<word>("OutputName", "D_ij_RSTE");

        Info<< "Calculating Viscous Diffusion Term for FO: " << this->name() << " at t=" << mesh().time().timeName() << endl;

        // Get nu
        if (!mesh().foundObject<IOdictionary>(transportPropsNm))
        {
            FatalErrorInFunction << "Transport properties '" << transportPropsNm << "' not found." << exit(FatalError);
        }
        const dictionary& transportProps = mesh().lookupObject<IOdictionary>(transportPropsNm);
        dimensionedScalar nu("nu", dimViscosity, transportProps);

        // Lookup RMean
        if (!mesh().foundObject<volSymmTensorField>(RMeanNameIn))
        {
            FatalErrorInFunction << "Field " << RMeanNameIn << " not found." << exit(FatalError);
        }
        const volSymmTensorField& RMean = mesh().lookupObject<volSymmTensorField>(RMeanNameIn);
        Info<< "  Inputs: " << RMean.name() << ", nu=" << nu << endl;

        volSymmTensorField D_ij_RSTE (
            IOobject(viscDiffOutNm, mesh().time().timeName(), mesh(), IOobject::NO_READ, IOobject::AUTO_WRITE),
            nu * fvc::laplacian(RMean)
        );

        const vectorField& cellCentres = mesh().C();
        const vectorField& faceCentres = mesh().Cf();

        forAll(mesh().boundary(), patchI)
        {
            const fvPatch& patch = mesh().boundary()[patchI];

            if (isA<wallFvPatch>(patch))
            {
                const label startFaceI = patch.start();
                const labelList& faceCells = patch.faceCells();

                forAll(faceCells, faceI)
                {
                    const label c0 = faceCells[faceI];
                    const label meshFaceI = startFaceI + faceI;

                    const point faceCenter = faceCentres[meshFaceI];
                    const vector normal = (cellCentres[c0] - faceCenter) / mag(cellCentres[c0] - faceCenter);

                    const labelList& nbrs = mesh().cellCells()[c0];
                    //Info << "Cell " << c0 << " (" << patch.name() << ") has " << nbrs.size() << " neighbors." << endl;

                    label bestNbr = -1;
                    scalar maxProj = -GREAT;

                    for (const label nbr : nbrs)
                    {
                        const vector d = cellCentres[nbr] - faceCenter;
                        scalar proj = d & normal;

                        //Info << "    Neighbor " << nbr << " projection: " << proj << endl;

                        if (proj > maxProj)
                        {
                            maxProj = proj;
                            bestNbr = nbr;
                        }
                    }

                    if (bestNbr != -1)
                    {
                        D_ij_RSTE[c0] = D_ij_RSTE[bestNbr];
                        //Info << "  Overwriting cell " << c0
                        //     << " from neighbor " << bestNbr
                        //     << " with projection " << maxProj << endl;
                    }
                    else
                    {
                        //Info << "  WARNING: No valid neighbor found for cell " << c0
                        //     << " on patch " << patch.name() << endl;
                    }
                }
            }
        }

        D_ij_RSTE.correctBoundaryConditions();
        D_ij_RSTE.write();

        Info<< "  Written field: " << D_ij_RSTE.name() << endl;
        Info<< endl;
    #};
}
