diff --git a/test/freeflow/navierstokes/unstructured/problem.hh b/test/freeflow/navierstokes/unstructured/problem.hh
index d769b9e55ab762d29d791c9f13d6e85692b873a9..fa18189fd5c70358e5e1100958042d7a9ef3c63a 100644
--- a/test/freeflow/navierstokes/unstructured/problem.hh
+++ b/test/freeflow/navierstokes/unstructured/problem.hh
@@ -32,6 +32,7 @@
 #include <dumux/io/vtk/function.hh>
 #include <dumux/io/grid/griddata.hh>
 #include <dumux/discretization/evalsolution.hh>
+#include <dumux/discretization/evalgradients.hh>
 #include <dumux/discretization/elementsolution.hh>
 
 namespace Dumux {
@@ -210,7 +211,6 @@ public:
     {
         auto fvGeometry = localView(this->gridGeometry());
         auto elemVolVars = localView(gridVariables.curGridVolVars());
-        auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
 
         BoundaryFluxes forces(0.0);
         Scalar cylinderSurface = 0.0;
@@ -223,24 +223,17 @@ public:
                 {
                     if (scvf.boundary() && onCylinderBoundary_(scvf.ipGlobal()))
                     {
-                        elemVolVars.bind(element, fvGeometry, v);
-                        elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
-                        const auto& fluxVarCache = elemFluxVarsCache[scvf];
-
-                        Dune::FieldMatrix<Scalar, 2, 2> stress(0.0);
-                        for (const auto& scv : scvs(fvGeometry))
-                        {
-                            const auto velocity = elemVolVars[scv].velocity();
-                            for (int dir = 0; dir < dim; ++dir)
-                                stress[dir].axpy(velocity[dir], fluxVarCache.gradN(scv.indexInElement()));
-                        }
-
+                        const auto elemSol = elementSolution(element, elemVolVars, fvGeometry);
+                        auto stress = evalGradients(element, element.geometry(), fvGeometry.gridGeometry(), elemSol, scvf.ipGlobal());
                         stress *= viscosity_;
 
+                        BoundaryFluxes normalStress(0.0);
                         for (int dir = 0; dir < dim; ++dir)
+                        {
                             stress[dir][dir] -= this->pressure(element, fvGeometry, scvf);
+                            normalStress[dir] = stress[dir]*scvf.unitOuterNormal();
+                        }
 
-                        BoundaryFluxes normalStress = mv(stress, scvf.unitOuterNormal());
                         forces.axpy(scvf.area(), normalStress);
 
                         cylinderSurface += scvf.area();
@@ -373,54 +366,29 @@ public:
                     if (scvf.boundary() && onCylinderBoundary_(scvf.ipGlobal()))
                     {
                         elemVolVars.bind(element, fvGeometry, v);
-                        Dune::FieldMatrix<Scalar, 2, 2> stress(0.0);
-                        const auto& localBasis = fvGeometry.feLocalBasis();
-                        Scalar area = 0.0;
+                        Dune::FieldVector<GlobalPosition, 2> stress(GlobalPosition(0.0));
 
                         const auto geometry = element.geometry();
                         const auto isgeometry = fvGeometry.geometry(scvf);
                         const auto& quad = Dune::QuadratureRules<Scalar, std::decay_t<decltype(isgeometry)>::mydimension>::rule(isgeometry.type(), 5);
                         for (auto&& qp : quad)
                         {
-                            Dune::FieldMatrix<Scalar, 2, 2> stressLocal(0.0);
-                            using FeLocalBasis = typename GridGeometry::FeCache::FiniteElementType::Traits::LocalBasisType;
-                            using ShapeJacobian = typename FeLocalBasis::Traits::JacobianType;
-                            using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
-                            using JacobianInverseTransposed = typename Element::Geometry::JacobianInverseTransposed;
-                            std::vector<GlobalPosition> gradN;
-                            std::vector<ShapeJacobian> shapeJacobian;
-                            std::vector<ShapeValue> shapeValues;
                             const auto ipGlobal = isgeometry.global(qp.position());
-                            const auto ipLocal = geometry.local(ipGlobal);
-                            JacobianInverseTransposed jacInvT = geometry.jacobianInverseTransposed(ipLocal);
-                            localBasis.evaluateJacobian(ipLocal, shapeJacobian);
-                            localBasis.evaluateFunction(ipLocal, shapeValues);
-
-                            // compute the gradN at for every scv/dof
-                            gradN.resize(fvGeometry.numScv());
-                            for (const auto& scv: scvs(fvGeometry))
-                                jacInvT.mv(shapeJacobian[scv.localDofIndex()][0], gradN[scv.indexInElement()]);
-
-                            for (const auto& scv : scvs(fvGeometry))
-                            {
-                                const auto velocity = elemVolVars[scv].velocity();
-                                for (int dir = 0; dir < dim; ++dir)
-                                    stressLocal[dir].axpy(velocity[dir], gradN[scv.indexInElement()]);
-                            }
-                            area += qp.weight()*isgeometry.integrationElement(qp.position());
-                            stressLocal += getTransposed(stressLocal);
+                            const auto elemSol = elementSolution(element, elemVolVars, fvGeometry);
+                            auto stressLocal = evalGradients(element, geometry, fvGeometry.gridGeometry(), elemSol, ipGlobal);
                             stressLocal *= qp.weight()*isgeometry.integrationElement(qp.position());
-
-                            stress += stressLocal;
+                            stress = stress + stressLocal;
                         }
                         stress *= viscosity_;
                         stress /= scvf.area();
 
-                        // Add pressure contribuition
+                        BoundaryFluxes normalStress(0.0);
                         for (int dir = 0; dir < dim; ++dir)
+                        {
+                            // Add pressure contribuition
                             stress[dir][dir] -= this->pressure(element, fvGeometry, scvf);
-
-                        BoundaryFluxes normalStress = mv(stress, scvf.unitOuterNormal());
+                            normalStress[dir] = stress[dir]*scvf.unitOuterNormal();
+                        }
                         forces.axpy(scvf.area(), normalStress);
 
                         cylinderSurface += scvf.area();