Shear rate calculation with STL input

Hello users and developers,

I am simulating a simple parabolic poiseuille flow case, and trying to compute shear rates. I input my geometry as an STL file. I get the expected velocity field. However, when I compute the shear rate, I get weird values of shear rate in the envelope regions. I have attached a visualization here. Any suggestions how to avoid this? At the exact boundary of my domain, it would help to use a forward or backward finite difference scheme as opposed to the centered finite differences used by computeStrainRate

Here is my writeVTK function:

void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
              T dx, T dt, plint iter)
{
    auto vel = *computeVelocity(lattice);
    ParallelVtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), 3, dx);
    vtkOut.writeData<T>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt);
    vtkOut.writeData<3,T>(vel, "velocity", dx/dt);
    vtkOut.writeData<T>(*computeSymmetricTensorNorm(*computeStrainRate(vel,  lattice.getBoundingBox()), lattice.getBoundingBox()), "shear", 1/dt);

}

Instead of computing the shear rate from velocity, I tried computing it from shear stress using computeStrainRateFromStress, and that seems to give me satisfactory results. Why this discrepancy?