diff --git a/test/freeflow/navierstokes/unstructured/main.cc b/test/freeflow/navierstokes/unstructured/main.cc
index 495baf0fddff5bca8c63df5c25afe0534374427f..f8530d8f71b2b0f6c54f58efea945dd042abbae4 100644
--- a/test/freeflow/navierstokes/unstructured/main.cc
+++ b/test/freeflow/navierstokes/unstructured/main.cc
@@ -151,7 +151,7 @@ int main(int argc, char** argv)
static constexpr double cDrafReference = 5.57953523384;
static constexpr double cLiftReference = 0.010618948146;
static constexpr double pDiffReference = 0.11752016697;
- const auto [cDrag, cLift] = momentumProblem->evalDragAndLiftCoefficientWithIntegration(*momentumGridVariables, x[momentumIdx]);
+ const auto [cDrag, cLift] = momentumProblem->evalDragAndLiftCoefficient(*momentumGridVariables, x[momentumIdx]);
std::cout << "cDrag: " << cDrag
<< " (reference: " << cDrafReference << ")"
<< "\n"
diff --git a/test/freeflow/navierstokes/unstructured/problem.hh b/test/freeflow/navierstokes/unstructured/problem.hh
index fa18189fd5c70358e5e1100958042d7a9ef3c63a..e4babbf0f657393987895e3616fe71f590a33c17 100644
--- a/test/freeflow/navierstokes/unstructured/problem.hh
+++ b/test/freeflow/navierstokes/unstructured/problem.hh
@@ -212,148 +212,6 @@ public:
auto fvGeometry = localView(this->gridGeometry());
auto elemVolVars = localView(gridVariables.curGridVolVars());
- BoundaryFluxes forces(0.0);
- Scalar cylinderSurface = 0.0;
- for (const auto& element : elements(this->gridGeometry().gridView()))
- {
- fvGeometry.bind(element);
- if (fvGeometry.hasBoundaryScvf())
- {
- for (const auto& scvf : scvfs(fvGeometry))
- {
- if (scvf.boundary() && onCylinderBoundary_(scvf.ipGlobal()))
- {
- 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();
- }
-
- forces.axpy(scvf.area(), normalStress);
-
- cylinderSurface += scvf.area();
- }
- }
- }
- }
-
- const Scalar dragForce = -forces[0];
- const Scalar liftForce = -forces[1];
-
- const Scalar uMean = 0.2;
- const Scalar lChar = 0.1;
- const Scalar coefficientFactor = 2.0/(uMean*uMean*lChar);
-
- // sanity check
- std::cout << "Reynolds number: " << uMean*lChar/viscosity_ << std::endl;
- std::cout << "Cylinder surface: " << cylinderSurface
- << " (reference: 0.31415926535)"
- << std::endl;
-
- return std::make_pair( coefficientFactor*dragForce, coefficientFactor*liftForce );
- }
-
- //! Computes drag and lift coefficient benchmark indicator
- //! (only works for correct domain size 2.2 x 0.41 and boundary conditions, Re=20)
- template<class SolutionVector, class GridVariables, class P = ParentType, typename std::enable_if_t<P::isMomentumProblem(), int> = 0>
- auto evalDragAndLiftCoefficientWithoutBubble(const GridVariables& gridVariables, const SolutionVector& v) const
- {
- auto fvGeometry = localView(this->gridGeometry());
- auto elemVolVars = localView(gridVariables.curGridVolVars());
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
-
- BoundaryFluxes forces(0.0);
- Scalar cylinderSurface = 0.0;
- for (const auto& element : elements(this->gridGeometry().gridView()))
- {
- fvGeometry.bind(element);
- if (fvGeometry.hasBoundaryScvf())
- {
- for (const auto& scvf : scvfs(fvGeometry))
- {
- if (scvf.boundary() && onCylinderBoundary_(scvf.ipGlobal()))
- {
- using CoordScalar = typename Element::Geometry::GlobalCoordinate::value_type;
- using P1FiniteElement = Dune::LagrangeCubeLocalFiniteElement<CoordScalar, Scalar, dim, 1>;
- P1FiniteElement p1FiniteElement;
-
- elemVolVars.bind(element, fvGeometry, v);
- const auto& localBasis = p1FiniteElement.localBasis();
-
- const auto geometry = element.geometry();
- Dune::FieldMatrix<Scalar, 2, 2> stress(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 = scvf.ipGlobal();
- const auto ipLocal = geometry.local(ipGlobal);
- JacobianInverseTransposed jacInvT = geometry.jacobianInverseTransposed(ipLocal);
- localBasis.evaluateJacobian(ipLocal, shapeJacobian);
- localBasis.evaluateFunction(ipLocal, shapeValues);
-
- shapeJacobian.push_back(ShapeJacobian(0.0));
- shapeValues.push_back(ShapeValue(0.0));
-
- // 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)
- stress[dir].axpy(velocity[dir], gradN[scv.indexInElement()]);
- }
-
- stress += getTransposed(stress);
- stress *= viscosity_;
-
- for (int dir = 0; dir < dim; ++dir)
- stress[dir][dir] -= this->pressure(element, fvGeometry, scvf);
-
- BoundaryFluxes normalStress = mv(stress, scvf.unitOuterNormal());
- forces.axpy(scvf.area(), normalStress);
-
- cylinderSurface += scvf.area();
- }
- }
- }
- }
-
- const Scalar dragForce = -forces[0];
- const Scalar liftForce = -forces[1];
-
- const Scalar uMean = 0.2;
- const Scalar lChar = 0.1;
- const Scalar coefficientFactor = 2.0/(uMean*uMean*lChar);
-
- // sanity check
- std::cout << "Reynolds number: " << uMean*lChar/viscosity_ << std::endl;
- std::cout << "Cylinder surface: " << cylinderSurface
- << " (reference: 0.31415926535)"
- << std::endl;
-
- return std::make_pair( coefficientFactor*dragForce, coefficientFactor*liftForce );
- }
-
- //! Computes drag and lift coefficient benchmark indicator
- //! (only works for correct domain size 2.2 x 0.41 and boundary conditions, Re=20)
- template<class SolutionVector, class GridVariables, class P = ParentType, typename std::enable_if_t<P::isMomentumProblem(), int> = 0>
- auto evalDragAndLiftCoefficientWithIntegration(const GridVariables& gridVariables, const SolutionVector& v) const
- {
- auto fvGeometry = localView(this->gridGeometry());
- auto elemVolVars = localView(gridVariables.curGridVolVars());
-
BoundaryFluxes forces(0.0);
Scalar cylinderSurface = 0.0;
for (const auto& element : elements(this->gridGeometry().gridView()))