diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index a80599495031fac69edfd231469b1783f94d863f..c8370bdfe610d418f2bf136240dd82a60cd45841 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -72,6 +72,7 @@ class FVProblem
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
     static constexpr bool isBox = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod == DiscretizationMethod::box;
+    static constexpr bool isStaggered = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod == DiscretizationMethod::staggered;
 
     using PointSourceMap = std::map<std::pair<std::size_t, std::size_t>,
                                     std::vector<PointSource> >;
@@ -206,9 +207,9 @@ public:
     PrimaryVariables dirichlet(const Element &element, const SubControlVolume &scv) const
     {
         // forward it to the method which only takes the global coordinate
-        if (!isBox)
+        if (!isBox && !isStaggered)
         {
-            DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for cell-centered method.");
+            DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for other than box or staggered method.");
         }
         else
             return asImp_().dirichletAtPos(scv.dofPosition());
diff --git a/dumux/common/staggeredfvproblem.hh b/dumux/common/staggeredfvproblem.hh
index ea82d38929e514487bdee09841d7210a628617d6..bed5fa94eb55d183ba88a7964b005021416e0f9a 100644
--- a/dumux/common/staggeredfvproblem.hh
+++ b/dumux/common/staggeredfvproblem.hh
@@ -130,20 +130,6 @@ public:
         return asImp_().neumannAtPos(scvf.ipGlobal());
     }
 
-    /*!
-     * \brief Evaluate the initial value for a control volume.
-     *
-     * \param globalPos The global position
-     */
-    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
-    {
-        // Throw an exception (there is no reasonable default value
-        // for initial values)
-        DUNE_THROW(Dune::InvalidStateException,
-                   "The problem does not provide "
-                   "an initial() or an initialAtPos() method.");
-    }
-
     /*!
      * \brief Evaluate the initial value for an element (for cell-centered primary variables)
      * or face (for velocities)
diff --git a/dumux/freeflow/compositional/staggered/fluxvariables.hh b/dumux/freeflow/compositional/staggered/fluxvariables.hh
index a39d4890b3870c67e576bd1a5e6a1da0a445f874..6d54535a222b19caf6a16e3bf3cae899453740a9 100644
--- a/dumux/freeflow/compositional/staggered/fluxvariables.hh
+++ b/dumux/freeflow/compositional/staggered/fluxvariables.hh
@@ -85,7 +85,7 @@ public:
             bool isOutflow = false;
             if(scvf.boundary())
             {
-                const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+                const auto bcTypes = problem.boundaryTypes(element, scvf);
                     if(bcTypes.isOutflow(eqIdx))
                         isOutflow = true;
             }
diff --git a/dumux/freeflow/compositional/staggered/localresidual.hh b/dumux/freeflow/compositional/staggered/localresidual.hh
index c0288bebcac7d02cf8d96ff5c63ef8a5286d469b..e7df94c045184ce677f161fb4b5fa80acf54a569 100644
--- a/dumux/freeflow/compositional/staggered/localresidual.hh
+++ b/dumux/freeflow/compositional/staggered/localresidual.hh
@@ -46,6 +46,9 @@ class FreeflowNCResidualImpl<TypeTag, DiscretizationMethod::staggered>
 
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using GridView = typename FVGridGeometry::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
@@ -103,11 +106,12 @@ public:
     template<class ElementVolumeVariables, class BoundaryTypes>
     void setFixedCell(CellCenterResidual& residual,
                       const Problem& problem,
+                      const Element& element,
                       const SubControlVolume& insideScv,
                       const ElementVolumeVariables& elemVolVars,
                       const BoundaryTypes& bcTypes) const
     {
-        ParentType::setFixedCell(residual, problem, insideScv, elemVolVars, bcTypes);
+        ParentType::setFixedCell(residual, problem, element, insideScv, elemVolVars, bcTypes);
 
         for (int compIdx = 0; compIdx < numComponents; ++compIdx)
         {
@@ -119,7 +123,7 @@ public:
             {
                 const auto& insideVolVars = elemVolVars[insideScv];
                 const Scalar massOrMoleFraction = useMoles ? insideVolVars.moleFraction(compIdx) : insideVolVars.massFraction(compIdx);
-                residual[eqIdx - cellCenterOffset] = massOrMoleFraction - problem.dirichletAtPos(insideScv.center())[eqIdx];
+                residual[eqIdx - cellCenterOffset] = massOrMoleFraction - problem.dirichlet(element, insideScv)[eqIdx];
             }
         }
 
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index 1518b7bde8a8d87b0392b931bf3fcf8c34e7277f..1cba7502f6713f1d804bb27bd92ef1e0175d4f21 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -145,7 +145,7 @@ public:
 
         // Check if we are on an outflow boundary.
         const bool isOutflow = scvf.boundary()
-                               ? problem.boundaryTypesAtPos(scvf.center()).isOutflow(Indices::conti0EqIdx)
+                               ? problem.boundaryTypes(element, scvf).isOutflow(Indices::conti0EqIdx)
                                : false;
 
         // Call the generic flux function.
@@ -242,9 +242,10 @@ public:
 
         // The pressure term.
         // If specified, the pressure can be normalized using the initial value on the scfv of interest.
+        // The scvf is used to normalize by the same value from the left and right side.
         // Can potentially help to improve the condition number of the system matrix.
         const Scalar pressure = normalizePressure ?
-                                insideVolVars.pressure() - problem.initialAtPos(scvf.center())[Indices::pressureIdx]
+                                insideVolVars.pressure() - problem.initial(scvf)[Indices::pressureIdx]
                               : insideVolVars.pressure();
 
         // Account for the orientation of the staggered face's normal outer normal vector
@@ -254,7 +255,7 @@ public:
         // Treat outflow conditions.
         if(scvf.boundary())
         {
-            if(problem.boundaryTypesAtPos(scvf.center()).isOutflow(Indices::velocity(scvf.directionIndex())))
+            if(problem.boundaryTypes(element, scvf).isOutflow(Indices::velocity(scvf.directionIndex())))
             {
                 // Treat the staggered half-volume adjacent to the boundary as if it was on the opposite side of the boundary.
                 // The respective face's outer normal vector will point in the same direction as the scvf's one.
@@ -559,10 +560,10 @@ private:
             outflow += velocitySelf * velocitySelf * upVolVars.density();
         }
 
-        // Apply a pressure at the boudary.
+        // Apply a pressure at the boundary.
         const Scalar boundaryPressure = normalizePressure
                                         ? (problem.dirichlet(element, scvf)[Indices::pressureIdx] -
-                                           problem.initialAtPos(scvf.center())[Indices::pressureIdx])
+                                           problem.initial(scvf)[Indices::pressureIdx])
                                         : problem.dirichlet(element, scvf)[Indices::pressureIdx];
         outflow += boundaryPressure;
 
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index 569227760fd59c4003da1cc186fd7449ca5cb2e6..aa82eda41512f0b1e8e6bd68c8d8cc7b89157574 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -189,6 +189,7 @@ public:
     template<class BoundaryTypes>
     void setFixedCell(CellCenterResidual& residual,
                       const Problem& problem,
+                      const Element& element,
                       const SubControlVolume& insideScv,
                       const ElementVolumeVariables& elemVolVars,
                       const BoundaryTypes& bcTypes) const
@@ -197,7 +198,7 @@ public:
         if(bcTypes.isDirichletCell(Indices::conti0EqIdx))
         {
             const auto& insideVolVars = elemVolVars[insideScv];
-            residual[Indices::conti0EqIdx - cellCenterOffset] = insideVolVars.pressure() - problem.dirichletAtPos(insideScv.center())[Indices::pressureIdx];
+            residual[Indices::conti0EqIdx - cellCenterOffset] = insideVolVars.pressure() - problem.dirichlet(element, insideScv)[Indices::pressureIdx];
         }
     }
 
@@ -248,7 +249,7 @@ protected:
 
                 // if specified, set a fixed value at the center of a cell at the boundary
                 const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
-                asImp_().setFixedCell(residual, problem, scv, elemVolVars, bcTypes);
+                asImp_().setFixedCell(residual, problem, element, scv, elemVolVars, bcTypes);
             }
         }
     }
diff --git a/dumux/freeflow/nonisothermal/localresidual.hh b/dumux/freeflow/nonisothermal/localresidual.hh
index ca84910ccea765e3c8e4efb2032c5497cfa50ba3..030fb2a78a6768b2084fb636ab2c22a9dce50c77 100644
--- a/dumux/freeflow/nonisothermal/localresidual.hh
+++ b/dumux/freeflow/nonisothermal/localresidual.hh
@@ -111,7 +111,7 @@ public:
         bool isOutflow = false;
         if(scvf.boundary())
         {
-            const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+            const auto bcTypes = problem.boundaryTypes(element, scvf);
             if(bcTypes.isOutflow(Indices::energyBalanceIdx))
                 isOutflow = true;
         }
diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
index 032f791ab9bffcb2e3b54f42d29fc7f3804efdf9..287fdb2c9d25f4e16ef4f2bf39e92f7fd3984b36 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
@@ -93,7 +93,7 @@ public:
                                                                       elemVolVars, elemFaceVars, scvf, fluxVarsCache);
 
         // calculate advective flux
-        const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+        const auto bcTypes = problem.boundaryTypes(element, scvf);
         const bool isOutflowK = scvf.boundary() && bcTypes.isOutflow(turbulentKineticEnergyEqIdx);
         const bool isOutflowEpsilon = scvf.boundary() && bcTypes.isOutflow(dissipationEqIdx);
         auto upwindTermK = [](const auto& volVars)