[Chaoran]

I'm confused about the implementation of the DarcyForchheimer equation in Openfoam. As I found in openfoam wiki https://openfoamwiki.net/index.php/DarcyForchheimer
the source term was like

or

But when I dig into the source code in DarcyForchheimer.H. I found something like

 forAll(cells, i)
         {
             const label celli = cells[i];
             const label j = this->fieldIndex(i);
             const tensor Cd =
                 mu[celli]*dZones[j] + (rho[celli]*mag(U[celli]))*fZones[j];
 
             const scalar isoCd = tr(Cd);
 
             Udiag[celli] += V[celli]*isoCd;
             Usource[celli] -= V[celli]*((Cd - I*isoCd) & U[celli]);
         }

where user defined values of D and f from the fvOptions dictionary are stored in the diagonal elements of the tensor and 0.5*f in .
I just can't understand the relations between the source code and the above equation. Why use the trace of tensor and the term .

And I want to implement a source term like

Can anybody help? Thanks in advance.