diff --git a/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
index 1ad7467efb81ae67a59a1f687f4721a5b1d0b65b..d64e47a13529f17c88dd40ad161d61d3b88aee97 100644
--- a/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
@@ -278,11 +278,9 @@ protected:
         // check if BJ-coeff is not zero to decide, if coupling condition
         // for the momentum balance (Dirichlet vor v.n) has to be set;
         // may be important at corners
-        const Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeff(this->element_(),
-                                                                              this->fvGeometry_(),
-                                                                              *isIt,
-                                                                              scvIdx,
-                                                                              boundaryFaceIdx);
+        const GlobalPosition &globalPos = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
+        Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeffAtPos(globalPos);
+
         // set velocity normal to the interface
         if (bcTypes.isCouplingInflow(momentumYIdx) && beaversJosephCoeff)
             this->residual_[scvIdx][momentumYIdx] =
@@ -325,11 +323,8 @@ protected:
     {
         const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
 
-        Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeff(this->element_(),
-                                                                        this->fvGeometry_(),
-                                                                        *isIt,
-                                                                        scvIdx,
-                                                                        boundaryFaceIdx);
+        const GlobalPosition &globalPos = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
+        Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeffAtPos(globalPos);
 
         // only enter here, if we are on a coupling boundary (see top)
         // and the BJ coefficient is not zero
diff --git a/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
index 48d51bcecb4e71486b3f8e09f3bebdde95495792..373546f19b10db8628ce38f8ba562902a051a7c1 100644
--- a/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
@@ -284,11 +284,9 @@ namespace Dumux
 			// check if BJ-coeff is not zero to decide, if coupling condition
 			// for the momentum balance (Dirichlet vor v.n) has to be set;
 			// may be important at corners
-			const Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeff(this->element_(),
-																				  this->fvGeometry_(),
-																				  *isIt,
-																				  scvIdx,
-																				  boundaryFaceIdx);
+      const GlobalPosition &globalPos = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
+      Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeffAtPos(globalPos);
+
 			// set velocity normal to the interface
 			if (bcTypes.isCouplingInflow(momentumYIdx) && beaversJosephCoeff)
 				this->residual_[scvIdx][momentumYIdx] =
@@ -336,12 +334,9 @@ namespace Dumux
 								const FluxVariables &boundaryVars) //TODO: required
 		{
 			const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
-			
-			Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeff(this->element_(),
-																			this->fvGeometry_(),
-																			*isIt,
-																			scvIdx,
-																			boundaryFaceIdx);
+
+      const GlobalPosition &globalPos = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
+      Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeffAtPos(globalPos);
 			
 			// only enter here, if we are on a coupling boundary (see top)
 			// and the BJ coefficient is not zero
diff --git a/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
index 327e0a70f23d0bee920bcb680528db864585205f..5f03262e77f5a05e7abff3b757707eff90e65400 100644
--- a/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
@@ -167,11 +167,11 @@ public:
         Scalar noDarcyX = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, NoDarcyX);
         Scalar xMin = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
 
-        bboxMin_[0] = std::max(xMin,noDarcyX);
-        bboxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+        bBoxMin_[0] = std::max(xMin,noDarcyX);
+        bBoxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
 
-        bboxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMin);
-        bboxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
+        bBoxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMin);
+        bBoxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
 
         runUpDistanceX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX); // first part of the interface without coupling
         initializationTime_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, InitTime);
@@ -234,14 +234,14 @@ public:
 
     /*!
      * \brief Specifies which kind of boundary condition should be
-     *        used for which equation on a given vertex
+     *        used for which equation on a given boundary segment
      *
      * \param values Stores the value of the boundary type
-     * \param vertex The vertex
-     */ 
-    void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+     * \param globalPos The global position
+     */
+    void boundaryTypesAtPos(BoundaryTypes &values,
+                            const GlobalPosition &globalPos) const
     {
-        const GlobalPosition globalPos = vertex.geometry().center();
         Scalar time = this->timeManager().time();
 
         values.setAllNeumann();
@@ -273,18 +273,15 @@ public:
     }
 
     /*!
-     * \brief Evaluate the boundary conditions for a dirichlet
-     *        control volume.
-     *
-     * \param values The Dirichlet values for the primary variables
-     * \param vertex The vertex representing the "half volume on the boundary"
+     * \brief Evaluate the boundary conditions for a Dirichlet
+     *        boundary segment
      *
-     * For this method, the \a values parameter stores primary variables.
+     * \param values Stores the Dirichlet values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variable} ] \f$
+     * \param globalPos The global position
      */
-    void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition globalPos = vertex.geometry().center();
-
         initial_(values, globalPos);
     }
 
@@ -294,21 +291,9 @@ public:
      *
      * \param values The Neumann values for the conservation equations in units of
      *                 \f$ [ \textnormal{unit of conserved quantity} / (m^{\textrm{dim}-1} \cdot s )] \f$
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param is The intersection between element and boundary
-     * \param scvIdx The local subcontrolvolume index
-     * \param boundaryFaceIdx The index of the boundary face
-     *
-     * For this method, the \a values parameter stores the mass flux
-     * in normal direction of each phase. Negative values mean influx.
+     * \param globalPos The global position
      */
-    void neumann(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const Intersection &is,
-                 const int scvIdx,
-                 const int boundaryFaceIdx) const
+    void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
         values = 0.;
     }
@@ -321,45 +306,27 @@ public:
     // \{
 
     /*!
-     * \brief Evaluate the source term for all phases within a given
-     *        sub-control-volume.
-     *
-     * \param values The source and sink values for the conservation equations in units of
-     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^\textrm{dim} \cdot s )] \f$
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param scvIdx The local subcontrolvolume index
+     * \brief Returns the source term
      *
-     * For this method, the \a values parameter stores the rate mass
-     * generated or annihilate per volume unit. Positive values mean
-     * that mass is created, negative ones mean that it vanishes.
+     * \param values Stores the source values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
+     * \param globalPos The global position
      */
-    void source(PrimaryVariables &values,
-                const Element &element,
-                const FVElementGeometry &fvGeometry,
-                const int scvIdx) const
+    void sourceAtPos(PrimaryVariables &values,
+                     const GlobalPosition &globalPos) const
     {
         values = Scalar(0);
     }
 
     /*!
-     * \brief Evaluate the initial value for a control volume.
+     * \brief Evaluates the initial values for a control volume
      *
-     * \param values The initial values for the primary variables
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param scvIdx The local subcontrolvolume index
-     *
-     * For this method, the \a values parameter stores primary
-     * variables.
+     * \param values Stores the initial values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variables} ] \f$
+     * \param globalPos The global position
      */
-    void initial(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const int scvIdx) const
+    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
-
         values = 0.;
 
         initial_(values, globalPos);
@@ -404,7 +371,7 @@ public:
                 storageChange /= (time - lastMassOutputTime_);
                 // 2d: interface length has to be accounted for
                 // in order to obtain kg/m²s
-                storageChange /= (bboxMax_[0]-bboxMin_[0]);
+                storageChange /= (bBoxMax_[0]-bBoxMin_[0]);
 
                 std::cout << "Time: " << time
                           << " TotalMass: " << storage[contiTotalMassIdx]
@@ -461,23 +428,23 @@ private:
     {
         //TODO: call density from fluidsystem
         values[pNIdx] = refPressure_
-                + 1000.*this->gravity()[1]*(globalPos[1]-bboxMax_[1]);
+                + 1000.*this->gravity()[1]*(globalPos[1]-bBoxMax_[1]);
 
         values[sWIdx] = initialSw_;
         values[temperatureIdx] = refTemperature_;
     }
 
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[0] < bboxMin_[0] + eps_; }
+    { return globalPos[0] < bBoxMin_[0] + eps_; }
 
     bool onRightBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[0] > bboxMax_[0] - eps_; }
+    { return globalPos[0] > bBoxMax_[0] - eps_; }
 
     bool onLowerBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[1] < bboxMin_[1] + eps_; }
+    { return globalPos[1] < bBoxMin_[1] + eps_; }
 
     bool onUpperBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[1] > bboxMax_[1] - eps_; }
+    { return globalPos[1] > bBoxMax_[1] - eps_; }
 
     bool onBoundary_(const GlobalPosition &globalPos) const
     {
@@ -486,8 +453,8 @@ private:
     }
 
     static constexpr Scalar eps_ = 1e-8;
-    GlobalPosition bboxMin_;
-    GlobalPosition bboxMax_;
+    GlobalPosition bBoxMin_;
+    GlobalPosition bBoxMax_;
 
     int freqMassOutput_;
 
diff --git a/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
index 35818bc8880fc07f7795b16604a3d13da45ead48..1b698e0c8569373fc085266fa04c3028401aecd1 100644
--- a/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
@@ -155,10 +155,10 @@ public:
         : ParentType(timeManager, gridView),
           spatialParams_(gridView)
     {
-        bboxMin_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
-        bboxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
-        bboxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
-        bboxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMax);
+        bBoxMin_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+        bBoxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+        bBoxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
+        bBoxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMax);
         runUpDistanceX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX); // first part of the interface without coupling
 
         refVelocity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefVelocity);
@@ -202,14 +202,14 @@ public:
 
     /*!
      * \brief Specifies which kind of boundary condition should be
-     *        used for which equation on a given vertex
+     *        used for which equation on a given boundary segment
      *
      * \param values Stores the value of the boundary type
-     * \param vertex The vertex
+     * \param globalPos The global position
      */
-    void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+    void boundaryTypesAtPos(BoundaryTypes &values,
+                            const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos = vertex.geometry().center();
         const Scalar time = this->timeManager().time();
 
         values.setAllDirichlet();
@@ -294,16 +294,15 @@ public:
     }
 
     /*!
-     * \brief Evaluates the boundary conditions for a Dirichlet vertex
+     * \brief Evaluates the boundary conditions for a Dirichlet
+     *        boundary segment
      *
      * \param values Stores the Dirichlet values for the conservation equations in
      *               \f$ [ \textnormal{unit of primary variable} ] \f$
-     * \param vertex The vertex
+     * \param globalPos The global position
      */
-    void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition globalPos = vertex.geometry().center();
-
         values = 0.0;
 
         FluidState fluidState;
@@ -315,34 +314,21 @@ public:
         values[velocityXIdx] = xVelocity_(globalPos);
         values[velocityYIdx] = 0.0;
         values[pressureIdx] = refPressure()  +
-                density*this->gravity()[1]*(globalPos[1] - bboxMin_[1]);
+                density*this->gravity()[1]*(globalPos[1] - bBoxMin_[1]);
         values[massOrMoleFracIdx] = refMassfrac();
         values[temperatureIdx] = refTemperature();
     }
 
     /*!
-     * \brief Evaluates the boundary conditions for a Neumann intersection
-     *
-     * \param values Stores the Neumann values for the conservation equations in
-     *               \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
-     * \param element The finite element
-     * \param fvGeometry The finite volume geometry of the element
-     * \param is The intersection between element and boundary
-     * \param scvIdx The local index of the sub-control volume
-     * \param boundaryFaceIdx The index of the boundary face
+     * \brief Evaluate the boundary conditions for a Neumann
+     *        boundary segment.
      *
-     * Negative values indicate an inflow.
+     * \param values The Neumann values for the conservation equations in units of
+     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^{\textrm{dim}-1} \cdot s )] \f$
+     * \param globalPos The global position
      */
-    void neumann(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const Intersection &is,
-                 const int scvIdx,
-                 const int boundaryFaceIdx) const
+    void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos =
-                fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
-
         values = 0.;
 
         FluidState fluidState;
@@ -355,7 +341,7 @@ public:
         const Scalar xVelocity = xVelocity_(globalPos);
 
         if (onLeftBoundary_(globalPos)
-                && globalPos[1] > bboxMin_[1] && globalPos[1] < bboxMax_[1])
+                && globalPos[1] > bBoxMin_[1] && globalPos[1] < bBoxMax_[1])
         {
             values[transportEqIdx] = -xVelocity*density*refMassfrac();
             values[energyEqIdx] = -xVelocity*density*enthalpy;
@@ -363,26 +349,14 @@ public:
     }
 
     /*!
-     * \brief Evaluate the Beavers-Joseph coefficient
-     *        at the center of a given intersection
+     * \brief Evaluate the Beavers-Joseph coefficient at given position
      *
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param is The intersection between element and boundary
-     * \param scvIdx The local subcontrolvolume index
-     * \param boundaryFaceIdx The index of the boundary face
+     * \param globalPos The global position
      *
      * \return Beavers-Joseph coefficient
      */
-    Scalar beaversJosephCoeff(const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const Intersection &is,
-                 const int scvIdx,
-                 const int boundaryFaceIdx) const
+    Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos =
-                fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
-
         if (onLowerBoundary_(globalPos))
             return alphaBJ_;
         else
@@ -408,23 +382,14 @@ public:
     // \}
 
     /*!
-     * \brief Evaluate the source term for all phases within a given
-     *        sub-control-volume.
-     *
-     * \param values The source and sink values for the conservation equations in units of
-     *               \f$ [ \textnormal{unit of conserved quantity} / (m^\textrm{dim} \cdot s )] \f$
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param scvIdx The local subcontrolvolume index
+     * \brief Returns the source term
      *
-     * For this method, the \a values parameter stores the rate mass
-     * generated or annihilate per volume unit. Positive values mean
-     * that mass is created, negative ones mean that it vanishes.
+     * \param values Stores the source values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
+     * \param globalPos The global position
      */
-    void source(PrimaryVariables &values,
-                const Element &element,
-                const FVElementGeometry &fvGeometry,
-                const int scvIdx) const
+    void sourceAtPos(PrimaryVariables &values,
+                     const GlobalPosition &globalPos) const
     {
         // ATTENTION: The source term of the mass balance has to be chosen as
         // div (q_momentum) in the problem file
@@ -434,22 +399,12 @@ public:
     /*!
      * \brief Evaluate the initial value for a control volume.
      *
-     * \param values The initial values for the primary variables
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param scvIdx The local subcontrolvolume index
-     *
-     * For this method, the \a values parameter stores primary
-     * variables.
+     * \param values Stores the initial values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variables} ] \f$
+     * \param globalPos The global position
      */
-    void initial(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const int scvIdx) const
+    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos
-            = element.geometry().corner(scvIdx);
-
         initial_(values, globalPos);
     }
     // \}
@@ -521,7 +476,7 @@ private:
         values[velocityYIdx] = 0.;
 
         values[pressureIdx] = refPressure()
-                + density*this->gravity()[1]*(globalPos[1] - bboxMin_[1]);
+                + density*this->gravity()[1]*(globalPos[1] - bBoxMin_[1]);
         values[massOrMoleFracIdx] = refMassfrac();
         values[temperatureIdx] = refTemperature();
     }
@@ -530,7 +485,7 @@ private:
     const Scalar xVelocity_(const GlobalPosition &globalPos) const
     {
         const Scalar vmax = refVelocity();
-        return  4*vmax*(globalPos[1] - bboxMin_[1])*(bboxMax_[1] - globalPos[1])
+        return  4*vmax*(globalPos[1] - bBoxMin_[1])*(bBoxMax_[1] - globalPos[1])
                 / (height_()*height_()) + 0.00134;
     }
 
@@ -559,16 +514,16 @@ private:
     { return sin(2*M_PI*this->timeManager().time()/period) * amplitude; }
 
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[0] < bboxMin_[0] + eps_; }
+    { return globalPos[0] < bBoxMin_[0] + eps_; }
 
     bool onRightBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[0] > bboxMax_[0] - eps_; }
+    { return globalPos[0] > bBoxMax_[0] - eps_; }
 
     bool onLowerBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[1] < bboxMin_[1] + eps_; }
+    { return globalPos[1] < bBoxMin_[1] + eps_; }
 
     bool onUpperBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[1] > bboxMax_[1] - eps_; }
+    { return globalPos[1] > bBoxMax_[1] - eps_; }
 
     bool onBoundary_(const GlobalPosition &globalPos) const
     {
@@ -577,15 +532,15 @@ private:
     }
 
     const Scalar height_() const
-    { return bboxMax_[1] - bboxMin_[1]; }
+    { return bBoxMax_[1] - bBoxMin_[1]; }
 
     // spatial parameters
     SpatialParams spatialParams_;
 
     static constexpr Scalar eps_ = 1e-8;
 
-    GlobalPosition bboxMin_;
-    GlobalPosition bboxMax_;
+    GlobalPosition bBoxMin_;
+    GlobalPosition bBoxMax_;
 
     Scalar refVelocity_;
     Scalar refPressure_;
diff --git a/test/multidomain/2cstokes2p2c/2p2csubproblem.hh b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
index 3d262f62781c4a29c64e8198642909d33e1c0f26..a920e71fc1c64be297f30126a02e6980c6b8c415 100644
--- a/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
+++ b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
@@ -147,11 +147,11 @@ public:
         Scalar noDarcyX = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, NoDarcyX);
         Scalar xMin = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
 
-        bboxMin_[0] = std::max(xMin,noDarcyX);
-        bboxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+        bBoxMin_[0] = std::max(xMin,noDarcyX);
+        bBoxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
 
-        bboxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMin);
-        bboxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
+        bBoxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMin);
+        bBoxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
 
         runUpDistanceX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX); // first part of the interface without coupling
         initializationTime_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, InitTime);
@@ -223,14 +223,14 @@ public:
 
     /*!
      * \brief Specifies which kind of boundary condition should be
-     *        used for which equation on a given vertex
+     *        used for which equation on a given boundary segment
      *
      * \param values Stores the value of the boundary type
-     * \param vertex The vertex
+     * \param globalPos The global position
      */
-    void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+    void boundaryTypesAtPos(BoundaryTypes &values,
+                            const GlobalPosition &globalPos) const
     {
-        const GlobalPosition globalPos = vertex.geometry().center();
         Scalar time = this->timeManager().time();
 
         values.setAllNeumann();
@@ -253,18 +253,15 @@ public:
     }
 
     /*!
-     * \brief Evaluate the boundary conditions for a dirichlet
-     *        control volume.
-     *
-     * \param values The Dirichlet values for the primary variables
-     * \param vertex The vertex representing the "half volume on the boundary"
+     * \brief Evaluate the boundary conditions for a Dirichlet
+     *        boundary segment
      *
-     * For this method, the \a values parameter stores primary variables.
+     * \param values Stores the Dirichlet values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variable} ] \f$
+     * \param globalPos The global position
      */
-    void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition globalPos = vertex.geometry().center();
-
         initial_(values, globalPos);
     }
 
@@ -274,21 +271,9 @@ public:
      *
      * \param values The Neumann values for the conservation equations in units of
      *                 \f$ [ \textnormal{unit of conserved quantity} / (m^{\textrm{dim}-1} \cdot s )] \f$
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param is The intersection between element and boundary
-     * \param scvIdx The local subcontrolvolume index
-     * \param boundaryFaceIdx The index of the boundary face
-     *
-     * For this method, the \a values parameter stores the mass flux
-     * in normal direction of each phase. Negative values mean influx.
+     * \param globalPos The global position
      */
-    void neumann(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const Intersection &is,
-                 const int scvIdx,
-                 const int boundaryFaceIdx) const
+    void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
         values = Scalar(0);
     }
@@ -301,45 +286,27 @@ public:
     // \{
 
     /*!
-     * \brief Evaluate the source term for all phases within a given
-     *        sub-control-volume.
-     *
-     * \param values The source and sink values for the conservation equations in units of
-     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^\textrm{dim} \cdot s )] \f$
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param scvIdx The local subcontrolvolume index
+     * \brief Returns the source term
      *
-     * For this method, the \a values parameter stores the rate mass
-     * generated or annihilate per volume unit. Positive values mean
-     * that mass is created, negative ones mean that it vanishes.
+     * \param values Stores the source values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
+     * \param globalPos The global position
      */
-    void source(PrimaryVariables &values,
-                const Element &element,
-                const FVElementGeometry &fvGeometry,
-                const int scvIdx) const
+    void sourceAtPos(PrimaryVariables &values,
+                     const GlobalPosition &globalPos) const
     {
         values = Scalar(0);
     }
 
     /*!
-     * \brief Evaluate the initial value for a control volume.
+     * \brief Evaluates the initial values for a control volume
      *
-     * \param values The initial values for the primary variables
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param scvIdx The local subcontrolvolume index
-     *
-     * For this method, the \a values parameter stores primary
-     * variables.
+     * \param values Stores the initial values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variables} ] \f$
+     * \param globalPos The global position
      */
-    void initial(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const int scvIdx) const
+    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
-
         values = Scalar(0);
 
         initial_(values, globalPos);
@@ -384,7 +351,7 @@ public:
                 storageChange /= (time - lastMassOutputTime_);
                 // 2d: interface length has to be accounted for
                 // in order to obtain kg/m²s
-                storageChange /= (bboxMax_[0]-bboxMin_[0]);
+                storageChange /= (bBoxMax_[0]-bBoxMin_[0]);
 
                 std::cout << "Time: " << time
                           << " TotalMass: " << storage[contiTotalMassIdx]
@@ -437,21 +404,21 @@ private:
                   const GlobalPosition &globalPos) const
     {
         values[pressureIdx] = refPressure_
-                + 1000.*this->gravity()[1]*(globalPos[1]-bboxMax_[1]);
+                + 1000.*this->gravity()[1]*(globalPos[1]-bBoxMax_[1]);
         values[switchIdx] = initialSw_;
     }
 
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[0] < bboxMin_[0] + eps_; }
+    { return globalPos[0] < bBoxMin_[0] + eps_; }
 
     bool onRightBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[0] > bboxMax_[0] - eps_; }
+    { return globalPos[0] > bBoxMax_[0] - eps_; }
 
     bool onLowerBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[1] < bboxMin_[1] + eps_; }
+    { return globalPos[1] < bBoxMin_[1] + eps_; }
 
     bool onUpperBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[1] > bboxMax_[1] - eps_; }
+    { return globalPos[1] > bBoxMax_[1] - eps_; }
 
     bool onBoundary_(const GlobalPosition &globalPos) const
     {
@@ -460,8 +427,8 @@ private:
     }
 
     static constexpr Scalar eps_ = 1e-8;
-    GlobalPosition bboxMin_;
-    GlobalPosition bboxMax_;
+    GlobalPosition bBoxMin_;
+    GlobalPosition bBoxMax_;
 
     int freqMassOutput_;
 
diff --git a/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
index 648e4da76f1b681dcf197f9176532b089ad50fc8..2414382d0a09f3c17050cfa089f11343c57001b7 100644
--- a/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
+++ b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
@@ -151,10 +151,10 @@ public:
         : ParentType(timeManager, gridView),
           spatialParams_(gridView)
     {
-        bboxMin_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
-        bboxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
-        bboxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
-        bboxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMax);
+        bBoxMin_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+        bBoxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+        bBoxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
+        bBoxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMax);
         runUpDistanceX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX); // first part of the interface without coupling
 
         refVelocity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefVelocity);
@@ -208,15 +208,14 @@ public:
 
     /*!
      * \brief Specifies which kind of boundary condition should be
-     *        used for which equation on a given vertex
+     *        used for which equation on a given boundary segment
      *
      * \param values Stores the value of the boundary type
-     * \param vertex The vertex
+     * \param globalPos The global position
      */
-    void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+    void boundaryTypesAtPos(BoundaryTypes &values,
+                            const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos = 
-                vertex.geometry().center();
         const Scalar time = this->timeManager().time();
 
         values.setAllDirichlet();
@@ -268,17 +267,15 @@ public:
     }
 
     /*!
-     * \brief Evaluates the boundary conditions for a Dirichlet vertex
+     * \brief Evaluates the boundary conditions for a Dirichlet
+     *        boundary segment
      *
      * \param values Stores the Dirichlet values for the conservation equations in
      *               \f$ [ \textnormal{unit of primary variable} ] \f$
-     * \param vertex The vertex
+     * \param globalPos The global position
      */
-    void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition globalPos = vertex.geometry().center();
-//        initial_(values, globalPos);
-
         FluidState fluidState;
         updateFluidStateForBC_(fluidState);
 
@@ -288,33 +285,20 @@ public:
         values[velocityXIdx] = xVelocity_(globalPos);
         values[velocityYIdx] = 0.0;
         values[pressureIdx] = refPressure()
-                + density*this->gravity()[1]*(globalPos[1] - bboxMin_[1]);
+                + density*this->gravity()[1]*(globalPos[1] - bBoxMin_[1]);
         values[massOrMoleFracIdx] = refMassfrac();
     }
 
     /*!
-     * \brief Evaluates the boundary conditions for a Neumann intersection
-     *
-     * \param values Stores the Neumann values for the conservation equations in
-     *               \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
-     * \param element The finite element
-     * \param fvGeometry The finite volume geometry of the element
-     * \param is The intersection between element and boundary
-     * \param scvIdx The local index of the sub-control volume
-     * \param boundaryFaceIdx The index of the boundary face
+     * \brief Evaluate the boundary conditions for a Neumann
+     *        boundary segment.
      *
-     * Negative values indicate an inflow.
+     * \param values The Neumann values for the conservation equations in units of
+     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^{\textrm{dim}-1} \cdot s )] \f$
+     * \param globalPos The global position
      */
-    void neumann(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const Intersection &is,
-                 const int scvIdx,
-                 const int boundaryFaceIdx) const
+    void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos =
-                fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
-
         values = 0.;
 
         FluidState fluidState;
@@ -325,7 +309,7 @@ public:
         const Scalar xVelocity = xVelocity_(globalPos);
 
         if (onLeftBoundary_(globalPos)
-                && globalPos[1] > bboxMin_[1] && globalPos[1] < bboxMax_[1])
+                && globalPos[1] > bBoxMin_[1] && globalPos[1] < bBoxMax_[1])
         {
         	// rho*v*X at inflow
         	values[transportEqIdx] = -xVelocity * density * refMassfrac();
@@ -333,26 +317,14 @@ public:
     }
 
     /*!
-     * \brief Evaluate the Beavers-Joseph coefficient
-     *        at the center of a given intersection
+     * \brief Evaluate the Beavers-Joseph coefficient at given position
      *
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param is The intersection between element and boundary
-     * \param scvIdx The local subcontrol volume index
-     * \param boundaryFaceIdx The index of the boundary face
+     * \param globalPos The global position
      *
      * \return Beavers-Joseph coefficient
      */
-    Scalar beaversJosephCoeff(const Element &element,
-                              const FVElementGeometry &fvGeometry,
-                              const Intersection &is,
-                              const int scvIdx,
-                              const int boundaryFaceIdx) const
+    Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos =
-                fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
-
         if (onLowerBoundary_(globalPos))
             return alphaBJ_;
         else
@@ -382,23 +354,14 @@ public:
     // \{
 
     /*!
-     * \brief Evaluate the source term for all phases within a given
-     *        sub-control-volume.
+     * \brief Returns the source term
      *
-     * \param values The source and sink values for the conservation equations in units of
-     *               \f$ [ \textnormal{unit of conserved quantity} / (m^\textrm{dim} \cdot s )] \f$
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param scvIdx The local subcontrolvolume index
-     *
-     * For this method, the \a values parameter stores the rate mass
-     * generated or annihilate per volume unit. Positive values mean
-     * that mass is created, negative ones mean that it vanishes.
+     * \param values Stores the source values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
+     * \param globalPos The global position
      */
-    void source(PrimaryVariables &values,
-                const Element &element,
-                const FVElementGeometry &fvGeometry,
-                const int scvIdx) const
+    void sourceAtPos(PrimaryVariables &values,
+                     const GlobalPosition &globalPos) const
     {
         // ATTENTION: The source term of the mass balance has to be chosen as
         // div (q_momentum) in the problem file
@@ -408,22 +371,12 @@ public:
     /*!
      * \brief Evaluate the initial value for a control volume.
      *
-     * \param values The initial values for the primary variables
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param scvIdx The local subcontrolvolume index
-     *
-     * For this method, the \a values parameter stores primary
-     * variables.
+     * \param values Stores the initial values for the conservation equations in
+     *               \f$ [ \textnormal{unit of primary variables} ] \f$
+     * \param globalPos The global position
      */
-    void initial(PrimaryVariables &values,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const int scvIdx) const
+    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
     {
-        const GlobalPosition &globalPos
-            = element.geometry().corner(scvIdx);
-
         initial_(values, globalPos);
     }
     // \}
@@ -495,7 +448,7 @@ private:
         values[velocityYIdx] = 0.;
 
         values[pressureIdx] = refPressure()
-                + density*this->gravity()[1]*(globalPos[1] - bboxMin_[1]);
+                + density*this->gravity()[1]*(globalPos[1] - bBoxMin_[1]);
         values[massOrMoleFracIdx] = refMassfrac();
     }
 
@@ -503,7 +456,7 @@ private:
     const Scalar xVelocity_(const GlobalPosition &globalPos) const
     {
         const Scalar vmax = refVelocity();
-        return  4*vmax*(globalPos[1] - bboxMin_[1])*(bboxMax_[1] - globalPos[1])
+        return  4*vmax*(globalPos[1] - bBoxMin_[1])*(bBoxMax_[1] - globalPos[1])
                 / (height_()*height_()) + 0.00134;
     }
 
@@ -532,16 +485,16 @@ private:
     { return sin(2*M_PI*this->timeManager().time()/period) * amplitude; }
 
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[0] < bboxMin_[0] + eps_; }
+    { return globalPos[0] < bBoxMin_[0] + eps_; }
 
     bool onRightBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[0] > bboxMax_[0] - eps_; }
+    { return globalPos[0] > bBoxMax_[0] - eps_; }
 
     bool onLowerBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[1] < bboxMin_[1] + eps_; }
+    { return globalPos[1] < bBoxMin_[1] + eps_; }
 
     bool onUpperBoundary_(const GlobalPosition &globalPos) const
-    { return globalPos[1] > bboxMax_[1] - eps_; }
+    { return globalPos[1] > bBoxMax_[1] - eps_; }
 
     bool onBoundary_(const GlobalPosition &globalPos) const
     {
@@ -551,14 +504,14 @@ private:
 
     // the height of the free-flow domain
     const Scalar height_() const
-    { return bboxMax_[1] - bboxMin_[1]; }
+    { return bBoxMax_[1] - bBoxMin_[1]; }
 
     // spatial parameters
     SpatialParams spatialParams_;
 
     static constexpr Scalar eps_ = 1e-8;
-    GlobalPosition bboxMin_;
-    GlobalPosition bboxMax_;
+    GlobalPosition bBoxMin_;
+    GlobalPosition bBoxMax_;
 
     Scalar refVelocity_;
     Scalar refPressure_;