const scalarField&     pb  = p.boundaryField()[patchID];
const symmTensorField& tb  = turbulence->devReff()().boundaryField()[patchID];
const vectorField&     sfb =  mesh.Sf().boundaryField()[patchID];

vector force = gSum( (sfb*pb) + (sfb & tb));

Info << "Calculated Force " << force << endl;

forceFile << 2*force.component(vector::X)/rdim/Udim/Udim/Adim << "   " 
          << 2*force.component(vector::Y)/rdim/Udim/Udim/Adim << "   " 
          << 2*force.component(vector::Z)/rdim/Udim/Udim/Adim << "   " 
          << force.component(vector::X) << "   " 
          << force.component(vector::Y) << "   " 
          << force.component(vector::Z) << "   " << endl;

