diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt
index 8901582125863d1de578ecad4b61823c55b4eacb..a97251962f8e61a825ae7d5b60ad11374bf59479 100644
--- a/dumux/CMakeLists.txt
+++ b/dumux/CMakeLists.txt
@@ -1,4 +1,5 @@
 add_subdirectory("common")
+add_subdirectory("discretization")
 add_subdirectory("freeflow")
 add_subdirectory("geomechanics")
 add_subdirectory("implicit")
diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..edb1a46d26076a164a4fa958deab18facd9dbd0b
--- /dev/null
+++ b/dumux/discretization/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_subdirectory("box")
+add_subdirectory("cellcentered")
diff --git a/dumux/discretization/box/CMakeLists.txt b/dumux/discretization/box/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b0883c5d5af56bdcb5c01cf7e5756cd95f338cff
--- /dev/null
+++ b/dumux/discretization/box/CMakeLists.txt
@@ -0,0 +1,8 @@
+install(FILES
+fluxvariablescachevector.hh
+fvelementgeometryvector.hh
+stencils.hh
+subcontrolvolumeface.hh
+subcontrolvolume.hh
+volumevariablesvector.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/box)
diff --git a/dumux/discretization/box/darcyslaw.hh b/dumux/discretization/box/darcyslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d36eb7f0aff9cccf4dd9b2defa29483fad268126
--- /dev/null
+++ b/dumux/discretization/box/darcyslaw.hh
@@ -0,0 +1,163 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ *        volume and mass fluxes of fluid phases over a face of a finite volume by means
+ *        of the Darcy approximation. Specializations are provided for the different discretization methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_DARCYS_LAW_HH
+#define DUMUX_DISCRETIZATION_BOX_DARCYS_LAW_HH
+
+#include <memory>
+
+#include <dune/common/float_cmp.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/implicit/properties.hh>
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration of properties
+NEW_PROP_TAG(ProblemEnableGravity);
+}
+
+/*!
+ * \ingroup DarcysLaw
+ * \brief Specialization of Darcy's Law for the box method.
+ */
+template <class TypeTag>
+class DarcysLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::Box>::type >
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
+    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
+    typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesCache) FluxVariablesCache;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::IndexSet::IndexType IndexType;
+    typedef typename GridView::ctype CoordScalar;
+    typedef std::vector<IndexType> Stencil;
+
+    enum { dim = GridView::dimension} ;
+    enum { dimWorld = GridView::dimensionworld} ;
+
+    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
+    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
+
+    typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> FeCache;
+    typedef typename FeCache::FiniteElementType::Traits::LocalBasisType FeLocalBasis;
+    typedef typename FluxVariablesCache::FaceData FaceData;
+
+public:
+
+    static Scalar flux(const Problem& problem,
+                       const Element& element,
+                       const SubControlVolumeFace& scvFace,
+                       const IndexType phaseIdx)
+    {
+        // get the precalculated local jacobian and shape values at the integration point
+        const auto& faceData = problem.model().fluxVarsCache()[scvFace].faceData();
+
+        const auto& fvGeometry = problem.model().fvGeometries(element);
+        const auto& insideScv = problem.model().fvGeometries().subControlVolume(scvFace.insideScvIdx());
+        const auto extrusionFactor = problem.model().curVolVars(insideScv).extrusionFactor();
+        const auto K = problem.spatialParams().intrinsicPermeability(insideScv);
+
+        // evaluate gradP - rho*g at integration point
+        DimVector gradP(0.0);
+        for (const auto& scv : fvGeometry.scvs())
+        {
+            // the global shape function gradient
+            DimVector gradI;
+            faceData.jacInvT.mv(faceData.localJacobian[scv.indexInElement()][0], gradI);
+
+            gradI *= problem.model().curVolVars(scv).pressure(phaseIdx);
+            gradP += gradI;
+        }
+        if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
+        {
+            // gravitational acceleration
+            DimVector g(problem.gravityAtPos(scvFace.center()));
+
+            // interpolate the density at the IP
+            const auto& insideVolVars = problem.model().curVolVars(insideScv);
+            const auto& outsideScv = problem.model().fvGeometries().subControlVolume(scvFace.outsideScvIdx());
+            const auto& outsideVolVars = problem.model().curVolVars(outsideScv);
+            Scalar rho = 0.5*(insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx));
+
+            // turn gravity into a force
+            g *= rho;
+
+            // subtract from pressure gradient
+            gradP -= g;
+        }
+
+        // apply the permeability and return the flux
+        auto KGradP = applyPermeability(K, gradP);
+        return -1.0*(KGradP*scvFace.unitOuterNormal())*scvFace.area()*extrusionFactor;
+    }
+
+    static DimVector applyPermeability(const DimWorldMatrix& K, const DimVector& gradI)
+    {
+        DimVector result(0.0);
+        K.mv(gradI, result);
+
+        return result;
+    }
+
+    static DimVector applyPermeability(const Scalar k, const DimVector& gradI)
+    {
+        DimVector result(gradI);
+        result *= k;
+        return result;
+    }
+
+    // This is for compatibility with the cc methods. The flux stencil info is obsolete for the box method.
+    static Stencil stencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
+    {
+        return Stencil(0);
+    }
+
+    static FaceData calculateFaceData(const Problem& problem, const Element& element, const typename Element::Geometry& geometry, const FeLocalBasis& localBasis, const SubControlVolumeFace& scvFace)
+    {
+        FaceData faceData;
+
+        // evaluate shape functions and gradients at the integration point
+        const auto ipLocal = geometry.local(scvFace.center());
+        faceData.jacInvT = geometry.jacobianInverseTransposed(ipLocal);
+        localBasis.evaluateJacobian(ipLocal, faceData.localJacobian);
+        localBasis.evaluateFunction(ipLocal, faceData.shapeValues);
+
+        return std::move(faceData);
+    }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/box/fickslaw.hh b/dumux/discretization/box/fickslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8406f61411ba489d18549e381b5c1b77924e6728
--- /dev/null
+++ b/dumux/discretization/box/fickslaw.hh
@@ -0,0 +1,179 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ *        diffusive mass fluxes due to molecular diffusion with Fick's law.
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_FICKS_LAW_HH
+#define DUMUX_DISCRETIZATION_BOX_FICKS_LAW_HH
+
+#include <dune/common/float_cmp.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/implicit/properties.hh>
+
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration of properties
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(FluidSystem);
+NEW_PROP_TAG(EffectiveDiffusivityModel);
+}
+
+/*!
+ * \ingroup CCTpfaFicksLaw
+ * \brief Specialization of Fick's Law for the box method.
+ */
+template <class TypeTag>
+class FicksLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::Box>::type >
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffDiffModel;
+    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
+    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::IndexSet::IndexType IndexType;
+    typedef typename std::vector<IndexType> Stencil;
+
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    enum { dim = GridView::dimension} ;
+    enum { dimWorld = GridView::dimensionworld} ;
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
+
+    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+
+    static Scalar flux(const Problem& problem,
+                       const Element& element,
+                       const SubControlVolumeFace& scvFace,
+                       const int phaseIdx,
+                       const int compIdx)
+    {
+        // diffusion tensors are always solution dependent
+        Scalar tij = calculateTransmissibility_(problem, scvFace, phaseIdx, compIdx);
+
+        // Get the inside volume variables
+        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+        const auto& insideVolVars = problem.model().curVolVars(insideScv);
+
+        // and the outside volume variables
+        const auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx());
+
+        // compute the diffusive flux
+        const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
+        const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
+        const auto rho = 0.5*(insideVolVars.molarDensity(phaseIdx) + outsideVolVars.molarDensity(phaseIdx));
+
+        return rho*tij*(xInside - xOutside);
+    }
+
+    static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace)
+    {
+        std::vector<IndexType> stencil;
+        stencil.clear();
+        if (!scvFace.boundary())
+        {
+            stencil.push_back(scvFace.insideScvIdx());
+            stencil.push_back(scvFace.outsideScvIdx());
+        }
+        else
+            stencil.push_back(scvFace.insideScvIdx());
+
+        return stencil;
+    }
+
+private:
+
+
+    static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace, const int phaseIdx, const int compIdx)
+    {
+        Scalar tij;
+
+        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+        const auto& insideVolVars = problem.model().curVolVars(insideScvIdx);
+
+        auto insideD = insideVolVars.diffusionCoefficient(phaseIdx, compIdx);
+        insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), insideD);
+        Scalar ti = calculateOmega_(problem, scvFace, insideD, insideScv);
+
+        if (!scvFace.boundary())
+        {
+            const auto outsideScvIdx = scvFace.outsideScvIdx();
+            const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
+            const auto& outsideVolVars = problem.model().curVolVars(outsideScvIdx);
+
+            auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx, compIdx);
+            outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), outsideD);
+            Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideD, outsideScv);
+
+            tij = scvFace.area()*(ti * tj)/(ti + tj);
+        }
+        else
+        {
+            tij = scvFace.area()*ti;
+        }
+
+        return tij;
+    }
+
+    static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &D, const SubControlVolume &scv)
+    {
+        GlobalPosition Dnormal;
+        D.mv(scvFace.unitOuterNormal(), Dnormal);
+
+        auto distanceVector = scvFace.center();
+        distanceVector -= scv.center();
+        distanceVector /= distanceVector.two_norm2();
+
+        Scalar omega = Dnormal * distanceVector;
+        omega *= problem.model().curVolVars(scv).extrusionFactor();
+
+        return omega;
+    }
+
+    static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, Scalar D, const SubControlVolume &scv)
+    {
+        auto distanceVector = scvFace.center();
+        distanceVector -= scv.center();
+        distanceVector /= distanceVector.two_norm2();
+
+        Scalar omega = D * (distanceVector * scvFace.unitOuterNormal());
+        omega *= problem.model().curVolVars(scv).extrusionFactor();
+
+        return omega;
+    }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/box/fluxvariablescachevector.hh b/dumux/discretization/box/fluxvariablescachevector.hh
similarity index 98%
rename from dumux/implicit/box/fluxvariablescachevector.hh
rename to dumux/discretization/box/fluxvariablescachevector.hh
index 1208662561943430a8fe4ce560cc5fc56bab3a4b..83e69beb5d4cc284941573a21312bfeaff6e1320 100644
--- a/dumux/implicit/box/fluxvariablescachevector.hh
+++ b/dumux/discretization/box/fluxvariablescachevector.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief Base class for the volume variables vector
  */
-#ifndef DUMUX_IMPLICIT_BOX_FLUXVARSCACHEVECTOR_HH
-#define DUMUX_IMPLICIT_BOX_FLUXVARSCACHEVECTOR_HH
+#ifndef DUMUX_DISCRETIZATION_BOX_FLUXVARSCACHEVECTOR_HH
+#define DUMUX_DISCRETIZATION_BOX_FLUXVARSCACHEVECTOR_HH
 
 #include <dumux/implicit/properties.hh>
 
diff --git a/dumux/implicit/box/fvelementgeometryvector.hh b/dumux/discretization/box/fvelementgeometryvector.hh
similarity index 99%
rename from dumux/implicit/box/fvelementgeometryvector.hh
rename to dumux/discretization/box/fvelementgeometryvector.hh
index b2da799aca34b53a0088aabf8c54deddaccbbfd5..48c11d229f87c976f351d25437dc735b2ae6de00 100644
--- a/dumux/implicit/box/fvelementgeometryvector.hh
+++ b/dumux/discretization/box/fvelementgeometryvector.hh
@@ -22,17 +22,14 @@
  *        This builds up the sub control volumes and sub control volume faces
  *        for each element.
  */
-#ifndef DUMUX_IMPLICIT_BOX_FV_GEOMETRY_VECTOR_HH
-#define DUMUX_IMPLICIT_BOX_FV_GEOMETRY_VECTOR_HH
+#ifndef DUMUX_DISCRETIZATION_BOX_FV_GEOMETRY_VECTOR_HH
+#define DUMUX_DISCRETIZATION_BOX_FV_GEOMETRY_VECTOR_HH
 
 #include <dune/geometry/multilineargeometry.hh>
 #include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
 #include <dumux/implicit/box/properties.hh>
-#include <dumux/implicit/subcontrolvolume.hh>
-#include <dumux/implicit/subcontrolvolumeface.hh>
-#include <dumux/implicit/fvelementgeometry.hh>
 #include <dumux/common/elementmap.hh>
 #include <dumux/common/math.hh>
 
diff --git a/dumux/implicit/box/stencils.hh b/dumux/discretization/box/stencils.hh
similarity index 98%
rename from dumux/implicit/box/stencils.hh
rename to dumux/discretization/box/stencils.hh
index 3185f8a55f83af0bcf5f39348ef3064057f13f6e..8b67b95636d501f453ae3f5a4be2883b1af2ae3a 100644
--- a/dumux/implicit/box/stencils.hh
+++ b/dumux/discretization/box/stencils.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief Implements the notion of stencils for vertex-centered models
  */
-#ifndef DUMUX_BOX_STENCILS_HH
-#define DUMUX_BOX_STENCILS_HH
+#ifndef DUMUX_DISCRETIZATION_BOX_STENCILS_HH
+#define DUMUX_DISCRETIZATION_BOX_STENCILS_HH
 
 #include <set>
 #include <dumux/implicit/box/properties.hh>
diff --git a/dumux/discretization/box/subcontrolvolume.hh b/dumux/discretization/box/subcontrolvolume.hh
new file mode 100644
index 0000000000000000000000000000000000000000..df886d74d351a46a4a3e58b029edd245d67efbda
--- /dev/null
+++ b/dumux/discretization/box/subcontrolvolume.hh
@@ -0,0 +1,84 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for a sub control volume
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_SUBCONTROLVOLUME_HH
+#define DUMUX_DISCRETIZATION_BOX_SUBCONTROLVOLUME_HH
+
+#include <dune/common/fvector.hh>
+#include <dumux/discretization/subcontrolvolumebase.hh>
+
+namespace Dumux
+{
+template<class G, typename I>
+class BoxSubControlVolume : public SubControlVolumeBase<G, I>
+{
+public:
+    // exported types
+    using Geometry = G;
+    using IndexType = I;
+
+private:
+    using Scalar = typename Geometry::ctype;
+    enum { dimworld = Geometry::coorddimension };
+    using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+public:
+    // the contructor in the box case
+    BoxSubControlVolume(Geometry geometry,
+                     IndexType scvIdx,
+                     IndexType elementIndex,
+                     IndexType indexInElement,
+                     IndexType dofIndex)
+    : SubControlVolumeBase<G, I>(std::move(geometry), std::move(elementIndex)),
+      scvIdx_(scvIdx), indexInElement_(std::move(indexInElement)),
+      dofIndex_(std::move(dofIndex)) {}
+
+    IndexType indexInElement() const
+    {
+        return indexInElement_;
+    }
+
+    //! The global index of this scv
+    IndexType index() const
+    {
+        return scvIdx_;
+    }
+
+    IndexType dofIndex() const
+    {
+        return dofIndex_;
+    }
+
+    GlobalPosition dofPosition() const
+    {
+        return this->geometry().corner(indexInElement_);
+    }
+
+private:
+    IndexType scvIdx_;
+    IndexType indexInElement_;
+    IndexType dofIndex_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/box/subcontrolvolumeface.hh b/dumux/discretization/box/subcontrolvolumeface.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b96d5085c58f0d20dbcac3e393b966cbb64aeb77
--- /dev/null
+++ b/dumux/discretization/box/subcontrolvolumeface.hh
@@ -0,0 +1,61 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for a sub control volume face
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_SUBCONTROLVOLUMEFACE_HH
+#define DUMUX_DISCRETIZATION_BOX_SUBCONTROLVOLUMEFACE_HH
+
+#include <utility>
+#include <dune/common/fvector.hh>
+#include <dumux/discretization/subcontrolvolumefacebase.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Discretization
+ * \brief Class for a sub control volume face in the box method, i.e a part of the boundary
+ *        of a sub control volume we compute fluxes on. We simply use the base class here.
+ */
+template<class Geometry, typename IndexType>
+class BoxSubControlVolumeFace : public SubControlVolumeFaceBase<Geometry, IndexType>
+{
+    using Scalar = typename Geometry::ctype;
+    static const int dim = Geometry::mydimension;
+    static const int dimworld = Geometry::coorddimension;
+
+    using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+    using LocalPosition = Dune::FieldVector<Scalar, dim>;
+
+public:
+    BoxSubControlVolumeFace(const Geometry& geometry,
+                             const GlobalPosition& ipGlobal,
+                             const GlobalPosition& unitOuterNormal,
+                             IndexType scvfIndex,
+                             IndexType indexInElement,
+                             const std::vector<IndexType>& scvIndices,
+                             bool boundary = false)
+    : SubControlVolumeFaceBase<Geometry, IndexType>(geometry, ipGlobal, unitOuterNormal, scvfIndex, indexInElement, scvIndices, boundary) {}
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/box/volumevariablesvector.hh b/dumux/discretization/box/volumevariablesvector.hh
similarity index 99%
rename from dumux/implicit/box/volumevariablesvector.hh
rename to dumux/discretization/box/volumevariablesvector.hh
index 0b97fdd5c6fe8ec86c8733bf31349937154383df..4d5e15a6867971c357cbb7da75f3cfb4975239a2 100644
--- a/dumux/implicit/box/volumevariablesvector.hh
+++ b/dumux/discretization/box/volumevariablesvector.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief Base class for the volume variables vector
  */
-#ifndef DUMUX_IMPLICIT_BOX_VOLVARSVECTOR_HH
-#define DUMUX_IMPLICIT_BOX_VOLVARSVECTOR_HH
+#ifndef DUMUX_DISCRETIZATION_BOX_VOLVARSVECTOR_HH
+#define DUMUX_DISCRETIZATION_BOX_VOLVARSVECTOR_HH
 
 #include <dumux/implicit/properties.hh>
 
diff --git a/dumux/discretization/cellcentered/CMakeLists.txt b/dumux/discretization/cellcentered/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..313cc2adfab0b789e000d868eec4c8c31c0f38f8
--- /dev/null
+++ b/dumux/discretization/cellcentered/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_subdirectory("tpfa")
+
+install(FILES
+fluxvariablescachevector.hh
+stencils.hh
+volumevariablesvector.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/cellcentered)
\ No newline at end of file
diff --git a/dumux/implicit/cellcentered/fluxvariablescachevector.hh b/dumux/discretization/cellcentered/fluxvariablescachevector.hh
similarity index 98%
rename from dumux/implicit/cellcentered/fluxvariablescachevector.hh
rename to dumux/discretization/cellcentered/fluxvariablescachevector.hh
index ffae3dbf9653ab3b948c24370c2aa64a1dfd71a2..13026e45e9b5c2c99839ac7c00375669866e7f2c 100644
--- a/dumux/implicit/cellcentered/fluxvariablescachevector.hh
+++ b/dumux/discretization/cellcentered/fluxvariablescachevector.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief Base class for the volume variables vector
  */
-#ifndef DUMUX_IMPLICIT_CC_FLUXVARSCACHEVECTOR_HH
-#define DUMUX_IMPLICIT_CC_FLUXVARSCACHEVECTOR_HH
+#ifndef DUMUX_DISCRETIZATION_CC_FLUXVARSCACHEVECTOR_HH
+#define DUMUX_DISCRETIZATION_CC_FLUXVARSCACHEVECTOR_HH
 
 #include <dumux/implicit/properties.hh>
 
diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/discretization/cellcentered/stencils.hh
similarity index 98%
rename from dumux/implicit/cellcentered/stencils.hh
rename to dumux/discretization/cellcentered/stencils.hh
index 34e156c31ff307455584f523512a63df77607421..b5de5588ef81b8b5720ac4b4b2f46db1852192c9 100644
--- a/dumux/implicit/cellcentered/stencils.hh
+++ b/dumux/discretization/cellcentered/stencils.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief Implements the notion of stencils for cell-centered models
  */
-#ifndef DUMUX_CC_STENCILS_HH
-#define DUMUX_CC_STENCILS_HH
+#ifndef DUMUX_DISCRETIZATION_CC_STENCILS_HH
+#define DUMUX_DISCRETIZATION_CC_STENCILS_HH
 
 #include <set>
 #include <dumux/implicit/cellcentered/properties.hh>
diff --git a/dumux/discretization/cellcentered/tpfa/CMakeLists.txt b/dumux/discretization/cellcentered/tpfa/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ff7ebeb2ba8ab4c57db241f4588d5f05ef4d7663
--- /dev/null
+++ b/dumux/discretization/cellcentered/tpfa/CMakeLists.txt
@@ -0,0 +1,5 @@
+install(FILES
+fvelementgeometryvector.hh
+subcontrolvolumeface.hh
+subcontrolvolume.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/cellcentered/tpfa)
\ No newline at end of file
diff --git a/dumux/porousmediumflow/constitutivelaws/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
similarity index 60%
rename from dumux/porousmediumflow/constitutivelaws/darcyslaw.hh
rename to dumux/discretization/cellcentered/tpfa/darcyslaw.hh
index f5f6fb07fd441175f3235c606e40b444f3ea1576..c24d224618ad2a754c52a16ac3c651c5ae1d71c0 100644
--- a/dumux/porousmediumflow/constitutivelaws/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
@@ -22,8 +22,8 @@
  *        volume and mass fluxes of fluid phases over a face of a finite volume by means
  *        of the Darcy approximation. Specializations are provided for the different discretization methods.
  */
-#ifndef DUMUX_POROUSMEDIUMFLOW_DARCYS_LAW_HH
-#define DUMUX_POROUSMEDIUMFLOW_DARCYS_LAW_HH
+#ifndef DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH
+#define DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH
 
 #include <memory>
 
@@ -47,15 +47,8 @@ NEW_PROP_TAG(ProblemEnableGravity);
 
 /*!
  * \ingroup DarcysLaw
- * \brief Evaluates the normal component of the Darcy velocity
- * on a (sub)control volume face. Specializations are provided
- * for the different discretization methods.
+ * \brief Specialization of Darcy's Law for the CCTpfa method.
  */
-template <class TypeTag, typename DiscretizationMethod = void>
-class DarcysLaw
-{};
-
-// Specialization for the CC-Tpfa method
 template <class TypeTag>
 class DarcysLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::CCTpfa>::type >
 {
@@ -198,116 +191,6 @@ private:
     }
 };
 
-// Specialization for the Box Method
-template <class TypeTag>
-class DarcysLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::Box>::type >
-{
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
-    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
-    typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesCache) FluxVariablesCache;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-
-    typedef typename GridView::template Codim<0>::Entity Element;
-    typedef typename GridView::IndexSet::IndexType IndexType;
-    typedef typename GridView::ctype CoordScalar;
-    typedef std::vector<IndexType> Stencil;
-
-    enum { dim = GridView::dimension} ;
-    enum { dimWorld = GridView::dimensionworld} ;
-
-    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
-    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
-
-    typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> FeCache;
-    typedef typename FeCache::FiniteElementType::Traits::LocalBasisType FeLocalBasis;
-    typedef typename FluxVariablesCache::FaceData FaceData;
-
-public:
-
-    static Scalar flux(const Problem& problem,
-                       const Element& element,
-                       const SubControlVolumeFace& scvFace,
-                       const IndexType phaseIdx)
-    {
-        // get the precalculated local jacobian and shape values at the integration point
-        const auto& faceData = problem.model().fluxVarsCache()[scvFace].faceData();
-
-        const auto& fvGeometry = problem.model().fvGeometries(element);
-        const auto& insideScv = problem.model().fvGeometries().subControlVolume(scvFace.insideScvIdx());
-        const auto extrusionFactor = problem.model().curVolVars(insideScv).extrusionFactor();
-        const auto K = problem.spatialParams().intrinsicPermeability(insideScv);
-
-        // evaluate gradP - rho*g at integration point
-        DimVector gradP(0.0);
-        for (const auto& scv : fvGeometry.scvs())
-        {
-            // the global shape function gradient
-            DimVector gradI;
-            faceData.jacInvT.mv(faceData.localJacobian[scv.indexInElement()][0], gradI);
-
-            gradI *= problem.model().curVolVars(scv).pressure(phaseIdx);
-            gradP += gradI;
-        }
-        if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
-        {
-            // gravitational acceleration
-            DimVector g(problem.gravityAtPos(scvFace.center()));
-
-            // interpolate the density at the IP
-            const auto& insideVolVars = problem.model().curVolVars(insideScv);
-            const auto& outsideScv = problem.model().fvGeometries().subControlVolume(scvFace.outsideScvIdx());
-            const auto& outsideVolVars = problem.model().curVolVars(outsideScv);
-            Scalar rho = 0.5*(insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx));
-
-            // turn gravity into a force
-            g *= rho;
-
-            // subtract from pressure gradient
-            gradP -= g;
-        }
-
-        // apply the permeability and return the flux
-        auto KGradP = applyPermeability(K, gradP);
-        return -1.0*(KGradP*scvFace.unitOuterNormal())*scvFace.area()*extrusionFactor;
-    }
-
-    static DimVector applyPermeability(const DimWorldMatrix& K, const DimVector& gradI)
-    {
-        DimVector result(0.0);
-        K.mv(gradI, result);
-
-        return result;
-    }
-
-    static DimVector applyPermeability(const Scalar k, const DimVector& gradI)
-    {
-        DimVector result(gradI);
-        result *= k;
-        return result;
-    }
-
-    // This is for compatibility with the cc methods. The flux stencil info is obsolete for the box method.
-    static Stencil stencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
-    {
-        return Stencil(0);
-    }
-
-    static FaceData calculateFaceData(const Problem& problem, const Element& element, const typename Element::Geometry& geometry, const FeLocalBasis& localBasis, const SubControlVolumeFace& scvFace)
-    {
-        FaceData faceData;
-
-        // evaluate shape functions and gradients at the integration point
-        const auto ipLocal = geometry.local(scvFace.center());
-        faceData.jacInvT = geometry.jacobianInverseTransposed(ipLocal);
-        localBasis.evaluateJacobian(ipLocal, faceData.localJacobian);
-        localBasis.evaluateFunction(ipLocal, faceData.shapeValues);
-
-        return std::move(faceData);
-    }
-};
-
 } // end namespace
 
 #endif
diff --git a/dumux/discretization/cellcentered/tpfa/fickslaw.hh b/dumux/discretization/cellcentered/tpfa/fickslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cddca62731be7c0ca1fb190a13b38dd7a785f4ed
--- /dev/null
+++ b/dumux/discretization/cellcentered/tpfa/fickslaw.hh
@@ -0,0 +1,178 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ *        diffusive mass fluxes due to molecular diffusion with Fick's law.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_TPFA_FICKS_LAW_HH
+#define DUMUX_DISCRETIZATION_CC_TPFA_FICKS_LAW_HH
+
+#include <dune/common/float_cmp.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/implicit/properties.hh>
+
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration of properties
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(FluidSystem);
+NEW_PROP_TAG(EffectiveDiffusivityModel);
+}
+
+/*!
+ * \ingroup CCTpfaFicksLaw
+ * \brief Specialization of Fick's Law for the CCTpfa method.
+ */
+template <class TypeTag>
+class FicksLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::CCTpfa>::type >
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffDiffModel;
+    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
+    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::IndexSet::IndexType IndexType;
+    typedef typename std::vector<IndexType> Stencil;
+
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    enum { dim = GridView::dimension} ;
+    enum { dimWorld = GridView::dimensionworld} ;
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
+
+    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+
+    static Scalar flux(const Problem& problem,
+                       const Element& element,
+                       const SubControlVolumeFace& scvFace,
+                       const int phaseIdx,
+                       const int compIdx)
+    {
+        // diffusion tensors are always solution dependent
+        Scalar tij = calculateTransmissibility_(problem, element, scvFace, phaseIdx, compIdx);
+
+        // Get the inside volume variables
+        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+        const auto& insideVolVars = problem.model().curVolVars(insideScv);
+
+        // and the outside volume variables
+        const auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx());
+
+        // compute the diffusive flux
+        const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
+        const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
+        const auto rho = 0.5*(insideVolVars.molarDensity(phaseIdx) + outsideVolVars.molarDensity(phaseIdx));
+
+        return rho*tij*(xInside - xOutside);
+    }
+
+    static Stencil stencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
+    {
+        std::vector<IndexType> stencil;
+        stencil.clear();
+        if (!scvFace.boundary())
+        {
+            stencil.push_back(scvFace.insideScvIdx());
+            stencil.push_back(scvFace.outsideScvIdx());
+        }
+        else
+            stencil.push_back(scvFace.insideScvIdx());
+
+        return stencil;
+    }
+
+private:
+
+
+    static Scalar calculateTransmissibility_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace, const int phaseIdx, const int compIdx)
+    {
+        Scalar tij;
+
+        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+        const auto& insideVolVars = problem.model().curVolVars(insideScvIdx);
+
+        auto insideD = insideVolVars.diffusionCoefficient(phaseIdx, compIdx);
+        insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), insideD);
+        Scalar ti = calculateOmega_(problem, element, scvFace, insideD, insideScv);
+
+        if (!scvFace.boundary())
+        {
+            const auto outsideScvIdx = scvFace.outsideScvIdx();
+            const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
+            const auto& outsideVolVars = problem.model().curVolVars(outsideScvIdx);
+
+            auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx, compIdx);
+            outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), outsideD);
+            Scalar tj = -1.0*calculateOmega_(problem, element, scvFace, outsideD, outsideScv);
+
+            tij = scvFace.area()*(ti * tj)/(ti + tj);
+        }
+        else
+        {
+            tij = scvFace.area()*ti;
+        }
+
+        return tij;
+    }
+
+    static Scalar calculateOmega_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace, const DimWorldMatrix &D, const SubControlVolume &scv)
+    {
+        GlobalPosition Dnormal;
+        D.mv(scvFace.unitOuterNormal(), Dnormal);
+
+        auto distanceVector = scvFace.center();
+        distanceVector -= scv.center();
+        distanceVector /= distanceVector.two_norm2();
+
+        Scalar omega = Dnormal * distanceVector;
+        omega *= problem.boxExtrusionFactor(element, scv);
+
+        return omega;
+    }
+
+    static Scalar calculateOmega_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace, Scalar D, const SubControlVolume &scv)
+    {
+        auto distanceVector = scvFace.center();
+        distanceVector -= scv.center();
+        distanceVector /= distanceVector.two_norm2();
+
+        Scalar omega = D * (distanceVector * scvFace.unitOuterNormal());
+        omega *= problem.boxExtrusionFactor(element, scv);
+
+        return omega;
+    }
+};
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh
similarity index 93%
rename from dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh
rename to dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh
index 13201044714489b151159f40f3ffb8249d336b65..1a239b4e38eee895288e0a3940b2ba40aac6250f 100644
--- a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh
@@ -22,12 +22,10 @@
  *        This builds up the sub control volumes and sub control volume faces
  *        for each element.
  */
-#ifndef DUMUX_IMPLICIT_TPFA_FV_GEOMETRY_VECTOR_HH
-#define DUMUX_IMPLICIT_TPFA_FV_GEOMETRY_VECTOR_HH
+#ifndef DUMUX_DISCRETIZATION_CCTPFA_FV_GEOMETRY_VECTOR_HH
+#define DUMUX_DISCRETIZATION_CCTPFA_FV_GEOMETRY_VECTOR_HH
 
-#include <dumux/implicit/subcontrolvolume.hh>
-#include <dumux/implicit/subcontrolvolumeface.hh>
-#include <dumux/implicit/fvelementgeometry.hh>
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
 #include <dumux/common/elementmap.hh>
 
 namespace Dumux
@@ -286,21 +284,21 @@ public:
             IndexType numLocalFaces = element.subEntities(1);
             std::vector<IndexType> scvfsIndexSet(numLocalFaces);
             std::vector<IndexType> neighborVolVarIndexSet(numLocalFaces);
+            IndexType localFaceIdx = 0;
             for (const auto& intersection : intersections(gridView_, element))
             {
-                IndexType localFaceIdx = intersection.indexInInside();
                 // inner sub control volume faces
                 if (intersection.neighbor())
                 {
                     scvfsIndexSet[localFaceIdx] = numScvf_++;
                     auto nIdx = problem.elementMapper().index(intersection.outside());
-                    neighborVolVarIndexSet[localFaceIdx] = nIdx;
+                    neighborVolVarIndexSet[localFaceIdx++] = nIdx;
                 }
                 // boundary sub control volume faces
                 else if (intersection.boundary())
                 {
                     scvfsIndexSet[localFaceIdx] = numScvf_++;
-                    neighborVolVarIndexSet[localFaceIdx] = numScvs_ + numBoundaryScvf_++;
+                    neighborVolVarIndexSet[localFaceIdx++] = numScvs_ + numBoundaryScvf_++;
                 }
             }
 
@@ -361,12 +359,12 @@ private:
             if (intersection.neighbor() || intersection.boundary())
             {
                 localScvfs_.push_back(SubControlVolumeFace(intersection.geometry(),
-                                                             intersection.geometry().center(),
-                                                             intersection.centerUnitOuterNormal(),
-                                                             scvFaceIndices[scvfCounter],
-                                                             intersection.indexInInside(),
-                                                             std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]}),
-                                                             intersection.boundary()));
+                                                           intersection.geometry().center(),
+                                                           intersection.centerUnitOuterNormal(),
+                                                           scvFaceIndices[scvfCounter],
+                                                           intersection.indexInInside(),
+                                                           std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]}),
+                                                           intersection.boundary()));
 
                 localScvfIndices_.push_back(scvFaceIndices[scvfCounter]);
                 scvfCounter++;
diff --git a/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh b/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh
new file mode 100644
index 0000000000000000000000000000000000000000..aaba6ee555717cba8bf670ebffc7a45bb9bc9222
--- /dev/null
+++ b/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh
@@ -0,0 +1,73 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for a sub control volume
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_TPFA_SUBCONTROLVOLUME_HH
+#define DUMUX_DISCRETIZATION_CC_TPFA_SUBCONTROLVOLUME_HH
+
+#include <dune/common/fvector.hh>
+#include <dumux/discretization/subcontrolvolumebase.hh>
+
+namespace Dumux
+{
+template<class G, typename I>
+class CCTpfaSubControlVolume : public SubControlVolumeBase<G, I>
+{
+public:
+    // exported types
+    using Geometry = G;
+    using IndexType = I;
+
+private:
+    using Scalar = typename Geometry::ctype;
+    enum { dimworld = Geometry::coorddimension };
+    using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+public:
+    // the contructor in the cc case
+    CCTpfaSubControlVolume(Geometry geometry,
+                     IndexType elementIndex)
+    : SubControlVolumeBase<G, I>(std::move(geometry), std::move(elementIndex)) {}
+
+    IndexType indexInElement() const
+    {
+        return IndexType(0);
+    }
+
+    //! The global index of this scv
+    IndexType index() const
+    {
+        return this->elementIndex();
+    }
+
+    IndexType dofIndex() const
+    {
+        return this->elementIndex();
+    }
+
+    GlobalPosition dofPosition() const
+    {
+        return this->geometry().center();
+    }
+};
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9a3d763c5c3d3cca5de94071ee97353608e2d30d
--- /dev/null
+++ b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
@@ -0,0 +1,62 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for a sub control volume face
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_TPFA_SUBCONTROLVOLUMEFACE_HH
+#define DUMUX_DISCRETIZATION_CC_TPFA_SUBCONTROLVOLUMEFACE_HH
+
+#include <utility>
+#include <dune/common/fvector.hh>
+#include <dumux/discretization/subcontrolvolumefacebase.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Discretization
+ * \brief Class for a sub control volume face in the box method, i.e a part of the boundary
+ *        of a sub control volume we compute fluxes on. We simply use the base class here.
+ */
+template<class Geometry, typename IndexType>
+class CCTpfaSubControlVolumeFace : public SubControlVolumeFaceBase<Geometry, IndexType>
+{
+    using Scalar = typename Geometry::ctype;
+    static const int dim = Geometry::mydimension;
+    static const int dimworld = Geometry::coorddimension;
+
+    using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+    using LocalPosition = Dune::FieldVector<Scalar, dim>;
+
+public:
+    CCTpfaSubControlVolumeFace(const Geometry& geometry,
+                               const GlobalPosition& ipGlobal,
+                               const GlobalPosition& unitOuterNormal,
+                               IndexType scvfIndex,
+                               IndexType indexInElement,
+                               const std::vector<IndexType>& scvIndices,
+                               bool boundary = false)
+    : SubControlVolumeFaceBase<Geometry, IndexType>(geometry, ipGlobal, unitOuterNormal, scvfIndex, indexInElement, scvIndices, boundary)
+    {}
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/cellcentered/volumevariablesvector.hh b/dumux/discretization/cellcentered/volumevariablesvector.hh
similarity index 99%
rename from dumux/implicit/cellcentered/volumevariablesvector.hh
rename to dumux/discretization/cellcentered/volumevariablesvector.hh
index bac9996b575c165da7e132541fd5f9e0b5eb84cf..85534726f162258f249310c5aa43a1a34bb03191 100644
--- a/dumux/implicit/cellcentered/volumevariablesvector.hh
+++ b/dumux/discretization/cellcentered/volumevariablesvector.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief Base class for the volume variables vector
  */
-#ifndef DUMUX_IMPLICIT_CC_VOLVARSVECTOR_HH
-#define DUMUX_IMPLICIT_CC_VOLVARSVECTOR_HH
+#ifndef DUMUX_DISCRETIZATION_CC_VOLVARSVECTOR_HH
+#define DUMUX_DISCRETIZATION_CC_VOLVARSVECTOR_HH
 
 #include <dumux/implicit/properties.hh>
 
diff --git a/dumux/discretization/darcyslaw.hh b/dumux/discretization/darcyslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..14e35ee628bfa57ce662ffe68a6b95cd0ed55137
--- /dev/null
+++ b/dumux/discretization/darcyslaw.hh
@@ -0,0 +1,47 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ *        volume and mass fluxes of fluid phases over a face of a finite volume by means
+ *        of the Darcy approximation. Specializations are provided for the different discretization methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_DARCYS_LAW_HH
+#define DUMUX_DISCRETIZATION_DARCYS_LAW_HH
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup DarcysLaw
+ * \brief Evaluates the normal component of the Darcy velocity
+ * on a (sub)control volume face. Specializations are provided
+ * for the different discretization methods. These specializations
+ * are found in the headers included below.
+ */
+template <class TypeTag, typename DiscretizationMethod = void>
+class DarcysLaw
+{};
+
+} // end namespace
+
+#include <dumux/discretization/box/darcyslaw.hh>
+#include <dumux/discretization/cellcentered/tpfa/darcyslaw.hh>
+
+#endif
diff --git a/dumux/discretization/fickslaw.hh b/dumux/discretization/fickslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a1de3ebc2dbdf1fa18255e7ce2906c17bcf80433
--- /dev/null
+++ b/dumux/discretization/fickslaw.hh
@@ -0,0 +1,59 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ *        diffusive mass fluxes due to molecular diffusion with Fick's law.
+ */
+#ifndef DUMUX_DISCRETIZATION_FICKS_LAW_HH
+#define DUMUX_DISCRETIZATION_FICKS_LAW_HH
+
+#include <dune/common/float_cmp.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/implicit/properties.hh>
+
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration of properties
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(FluidSystem);
+NEW_PROP_TAG(EffectiveDiffusivityModel);
+}
+
+/*!
+ * \ingroup CCTpfaFicksLaw
+ * \brief Evaluates the diffusive mass flux according to Fick's law
+ */
+template <class TypeTag, typename DiscretizationMethod = void>
+class FicksLaw
+{};
+
+} // end namespace
+
+#include <dumux/discretization/cellcentered/tpfa/fickslaw.hh>
+#include <dumux/discretization/box/fickslaw.hh>
+
+#endif
diff --git a/dumux/implicit/fluxvariablesbase.hh b/dumux/discretization/fluxvariablesbase.hh
similarity index 97%
rename from dumux/implicit/fluxvariablesbase.hh
rename to dumux/discretization/fluxvariablesbase.hh
index 6fe0a78696277116e378ed6eec9c1a83d1029a83..de8061254d48ce12d473e717c04a6ee5430b64b1 100644
--- a/dumux/implicit/fluxvariablesbase.hh
+++ b/dumux/discretization/fluxvariablesbase.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief Base class for the flux variables
  */
-#ifndef DUMUX_IMPLICIT_FLUXVARIABLESBASE_HH
-#define DUMUX_IMPLICIT_FLUXVARIABLESBASE_HH
+#ifndef DUMUX_DISCRETIZATION_FLUXVARIABLESBASE_HH
+#define DUMUX_DISCRETIZATION_FLUXVARIABLESBASE_HH
 
 #include <dumux/implicit/properties.hh>
 
@@ -29,7 +29,7 @@ namespace Dumux
 {
 
 /*!
- * \ingroup ImplicitModel
+ * \ingroup Discretization
  * \brief Base class for the flux variables
  *        Actual flux variables inherit from this class
  */
diff --git a/dumux/implicit/fvelementgeometry.hh b/dumux/discretization/fvelementgeometry.hh
similarity index 97%
rename from dumux/implicit/fvelementgeometry.hh
rename to dumux/discretization/fvelementgeometry.hh
index df59ef883b83a6e500551eb469ea18e64fd1d7c7..b2b5cf8e9b88792417a04edda26e110fba2640e1 100644
--- a/dumux/implicit/fvelementgeometry.hh
+++ b/dumux/discretization/fvelementgeometry.hh
@@ -21,15 +21,12 @@
  * \brief Class providing iterators over sub control volumes and sub control
  *        volume faces of an element.
  */
-#ifndef DUMUX_FV_ELEMENTGEOMETRY_HH
-#define DUMUX_FV_ELEMENTGEOMETRY_HH
+#ifndef DUMUX_DISCRETIZATION_FV_ELEMENTGEOMETRY_HH
+#define DUMUX_DISCRETIZATION_FV_ELEMENTGEOMETRY_HH
 
 #include <dune/common/iteratorrange.hh>
 #include <dune/common/iteratorfacades.hh>
 
-#include <dumux/implicit/subcontrolvolume.hh>
-#include <dumux/implicit/subcontrolvolumeface.hh>
-
 namespace Dumux
 {
 namespace Properties
@@ -41,7 +38,7 @@ NEW_PROP_TAG(FVElementGeometryVector);
 }
 
 /*!
- * \ingroup ImplcititModel
+ * \ingroup Discretization
  * \brief An iterator over sub control volumes
  */
 template<class SubControlVolume, class Vector, class FVElementGeometryVector>
diff --git a/dumux/implicit/subcontrolvolume.hh b/dumux/discretization/subcontrolvolumebase.hh
similarity index 100%
rename from dumux/implicit/subcontrolvolume.hh
rename to dumux/discretization/subcontrolvolumebase.hh
diff --git a/dumux/implicit/subcontrolvolumeface.hh b/dumux/discretization/subcontrolvolumefacebase.hh
similarity index 87%
rename from dumux/implicit/subcontrolvolumeface.hh
rename to dumux/discretization/subcontrolvolumefacebase.hh
index f0f964e81d54c8a73c1d9cf418c953e1ff959c88..ab732275d7e7c33d073146635603ee9b5d627879 100644
--- a/dumux/implicit/subcontrolvolumeface.hh
+++ b/dumux/discretization/subcontrolvolumefacebase.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief Base class for a sub control volume face
  */
-#ifndef DUMUX_SUBCONTROLVOLUMEFACE_HH
-#define DUMUX_SUBCONTROLVOLUMEFACE_HH
+#ifndef DUMUX_DISCRETIZATION_SUBCONTROLVOLUMEFACEBASE_HH
+#define DUMUX_DISCRETIZATION_SUBCONTROLVOLUMEFACEBASE_HH
 
 #include <utility>
 #include <dune/common/fvector.hh>
@@ -30,12 +30,12 @@ namespace Dumux
 {
 
 /*!
- * \ingroup ImplicitModel
+ * \ingroup Discretization
  * \brief Base class for a sub control volume face, i.e a part of the boundary
  *        of a sub control volume we computing a flux on.
  */
 template<class Geometry, typename IndexType>
-class SubControlVolumeFace
+class SubControlVolumeFaceBase
 {
     using Scalar = typename Geometry::ctype;
     static const int dim = Geometry::mydimension;
@@ -45,13 +45,13 @@ class SubControlVolumeFace
     using LocalPosition = Dune::FieldVector<Scalar, dim>;
 
 public:
-    SubControlVolumeFace(const Geometry& geometry,
-                         const GlobalPosition& ipGlobal,
-                         const GlobalPosition& unitOuterNormal,
-                         IndexType scvfIndex,
-                         IndexType indexInElement,
-                         const std::vector<IndexType>& scvIndices,
-                         bool boundary = false)
+    SubControlVolumeFaceBase(const Geometry& geometry,
+                             const GlobalPosition& ipGlobal,
+                             const GlobalPosition& unitOuterNormal,
+                             IndexType scvfIndex,
+                             IndexType indexInElement,
+                             const std::vector<IndexType>& scvIndices,
+                             bool boundary = false)
     : geometry_(geometry),
       ipGlobal_(ipGlobal),
       ipLocal_(geometry.local(ipGlobal)),
@@ -142,4 +142,5 @@ private:
 
 } // end namespace
 
+
 #endif
diff --git a/dumux/implicit/volumevariables.hh b/dumux/discretization/volumevariables.hh
similarity index 97%
rename from dumux/implicit/volumevariables.hh
rename to dumux/discretization/volumevariables.hh
index fe3b89b594ae12dc0395565bfceba8907482dc08..7acd2069a9a46354041ce230b7e3cdeb6b6cf5ff 100644
--- a/dumux/implicit/volumevariables.hh
+++ b/dumux/discretization/volumevariables.hh
@@ -21,10 +21,10 @@
  * \brief Base class for the model specific class which provides
  *        access to all volume averaged quantities.
  */
-#ifndef DUMUX_IMPLICIT_VOLUME_VARIABLES_HH
-#define DUMUX_IMPLICIT_VOLUME_VARIABLES_HH
+#ifndef DUMUX_DISCRETIZATION_VOLUME_VARIABLES_HH
+#define DUMUX_DISCRETIZATION_VOLUME_VARIABLES_HH
 
-#include "properties.hh"
+#include <dumux/implicit/properties.hh>
 
 #include <dumux/common/valgrind.hh>
 
diff --git a/dumux/implicit/box/propertydefaults.hh b/dumux/implicit/box/propertydefaults.hh
index 1eb743f5e83a6260b53966ce098d3e6776aeab14..3681ef036c5559759336bcbc84057b2c62819ed6 100644
--- a/dumux/implicit/box/propertydefaults.hh
+++ b/dumux/implicit/box/propertydefaults.hh
@@ -27,18 +27,20 @@
 #define DUMUX_BOX_PROPERTY_DEFAULTS_HH
 
 #include <dumux/implicit/propertydefaults.hh>
-#include <dumux/implicit/fvelementgeometry.hh>
+#include <dumux/discretization/fvelementgeometry.hh>
+#include <dumux/discretization/box/subcontrolvolume.hh>
+#include <dumux/discretization/box/subcontrolvolumeface.hh>
+#include <dumux/discretization/box/fluxvariablescachevector.hh>
+#include <dumux/discretization/box/volumevariablesvector.hh>
+#include <dumux/discretization/box/fvelementgeometryvector.hh>
+#include <dumux/discretization/box/stencils.hh>
 #include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
 
-#include "fluxvariablescachevector.hh"
-#include "volumevariablesvector.hh"
 #include "elementboundarytypes.hh"
-#include "fvelementgeometryvector.hh"
 #include "localresidual.hh"
 #include "localjacobian.hh"
 #include "assembler.hh"
 #include "properties.hh"
-#include "stencils.hh"
 
 namespace Dumux {
 
@@ -69,7 +71,7 @@ private:
     using ScvGeometry = Dune::MultiLinearGeometry<Scalar, dim, dimWorld>;
     using IndexType = typename GridView::IndexSet::IndexType;
 public:
-    using type = SubControlVolume<ScvGeometry, IndexType, /*isBox=*/true>;
+    using type = BoxSubControlVolume<ScvGeometry, IndexType, /*isBox=*/true>;
 };
 
 SET_PROP(BoxModel, SubControlVolumeFace)
@@ -82,7 +84,7 @@ private:
     using ScvfGeometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld>;
     using IndexType = typename GridView::IndexSet::IndexType;
 public:
-    using type = SubControlVolumeFace<ScvfGeometry, IndexType>;
+    using type = BoxSubControlVolumeFace<ScvfGeometry, IndexType>;
 };
 
 //! Set the default for the ElementBoundaryTypes
@@ -95,22 +97,22 @@ SET_TYPE_PROP(BoxModel, DofMapper, typename GET_PROP_TYPE(TypeTag, VertexMapper)
 SET_TYPE_PROP(BoxModel, StencilsVector, BoxStencilsVector<TypeTag>);
 
 //! The global current volume variables vector class
-SET_TYPE_PROP(BoxModel, CurrentVolumeVariablesVector, Dumux::BoxVolumeVariablesVector<TypeTag, false, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
+SET_TYPE_PROP(BoxModel, CurrentVolumeVariablesVector, BoxVolumeVariablesVector<TypeTag, false, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
 
 //! The global previous volume variables vector class
-SET_TYPE_PROP(BoxModel, PreviousVolumeVariablesVector, Dumux::BoxVolumeVariablesVector<TypeTag, true, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
+SET_TYPE_PROP(BoxModel, PreviousVolumeVariablesVector, BoxVolumeVariablesVector<TypeTag, true, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
 
 //! The global flux variables cache vector class
-SET_TYPE_PROP(BoxModel, FluxVariablesCacheVector, Dumux::BoxFluxVariablesCacheVector<TypeTag>);
+SET_TYPE_PROP(BoxModel, FluxVariablesCacheVector, BoxFluxVariablesCacheVector<TypeTag>);
 
 //! Set the BaseLocalResidual to BoxLocalResidual
 SET_TYPE_PROP(BoxModel, BaseLocalResidual, BoxLocalResidual<TypeTag>);
 
 //! Assembler for the global jacobian matrix
-SET_TYPE_PROP(BoxModel, JacobianAssembler, Dumux::BoxAssembler<TypeTag>);
+SET_TYPE_PROP(BoxModel, JacobianAssembler, BoxAssembler<TypeTag>);
 
 //! The local jacobian operator
-SET_TYPE_PROP(BoxModel, LocalJacobian, Dumux::BoxLocalJacobian<TypeTag>);
+SET_TYPE_PROP(BoxModel, LocalJacobian, BoxLocalJacobian<TypeTag>);
 
 //! indicate that this is a box discretization
 SET_BOOL_PROP(BoxModel, ImplicitIsBox, true);
diff --git a/dumux/implicit/cellcentered/propertydefaults.hh b/dumux/implicit/cellcentered/propertydefaults.hh
index 72cfb72f4587b484adb272ab05f3deff8274a39b..ef9432377183dafd17437e249897e20369114569 100644
--- a/dumux/implicit/cellcentered/propertydefaults.hh
+++ b/dumux/implicit/cellcentered/propertydefaults.hh
@@ -28,13 +28,13 @@
 #define DUMUX_CC_PROPERTY_DEFAULTS_HH
 
 #include <dumux/implicit/propertydefaults.hh>
+#include <dumux/discretization/cellcentered/fluxvariablescachevector.hh>
+#include <dumux/discretization/cellcentered/volumevariablesvector.hh>
+#include <dumux/discretization/cellcentered/stencils.hh>
 
-#include "fluxvariablescachevector.hh"
-#include "volumevariablesvector.hh"
 #include "elementboundarytypes.hh"
 #include "localresidual.hh"
 #include "properties.hh"
-#include "stencils.hh"
 #include "localjacobian.hh"
 #include "assembler.hh"
 
diff --git a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh
index 48af42db4dd5b4220a6669bccd7acee75cce5cf4..cc904e2b56783dc6a1d67589904e457e8811a284 100644
--- a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh
+++ b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh
@@ -28,9 +28,10 @@
 #define DUMUX_CCTPFA_PROPERTY_DEFAULTS_HH
 
 #include <dumux/implicit/propertydefaults.hh>
-#include <dumux/implicit/fvelementgeometry.hh>
 #include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
-#include <dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh>
+#include <dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh>
+#include <dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh>
+#include <dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh>
 #include <dumux/implicit/cellcentered/properties.hh>
 
 namespace Dumux {
@@ -53,7 +54,7 @@ private:
     using ScvGeometry = typename Grid::template Codim<0>::Geometry;
     using IndexType = typename Grid::LeafGridView::IndexSet::IndexType;
 public:
-    typedef Dumux::SubControlVolume<ScvGeometry, IndexType, /*isBox=*/false> type;
+    typedef Dumux::CCTpfaSubControlVolume<ScvGeometry, IndexType> type;
 };
 
 SET_PROP(CCTpfaModel, SubControlVolumeFace)
@@ -63,7 +64,7 @@ private:
     using ScvfGeometry = typename Grid::template Codim<1>::Geometry;
     using IndexType = typename Grid::LeafGridView::IndexSet::IndexType;
 public:
-    typedef Dumux::SubControlVolumeFace<ScvfGeometry, IndexType> type;
+    typedef Dumux::CCTpfaSubControlVolumeFace<ScvfGeometry, IndexType> type;
 };
 
 } // namespace Properties
diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh
index 6fa48bf48053184452c7c15d76b4e46f63e1c1ce..2d4219e401879e44f97545bcfb3f40e5e81f5644 100644
--- a/dumux/implicit/propertydefaults.hh
+++ b/dumux/implicit/propertydefaults.hh
@@ -38,15 +38,15 @@
 #include <dumux/porousmediumflow/implicit/fluxvariables.hh>
 #include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
 
-#include <dumux/porousmediumflow/constitutivelaws/darcyslaw.hh>
-#include <dumux/porousmediumflow/constitutivelaws/fickslaw.hh>
+#include <dumux/discretization/volumevariables.hh>
+#include <dumux/discretization/fvelementgeometry.hh>
+#include <dumux/discretization/darcyslaw.hh>
+#include <dumux/discretization/fickslaw.hh>
 
 #include "properties.hh"
 #include "model.hh"
 #include "assembler.hh"
 #include "localjacobian.hh"
-#include "volumevariables.hh"
-#include "fvelementgeometry.hh"
 
 namespace Dumux {
 
diff --git a/dumux/porousmediumflow/1p/implicit/volumevariables.hh b/dumux/porousmediumflow/1p/implicit/volumevariables.hh
index f83411592953f7779e8feb4281647004dbad4274..d888b631cdc5b59eb26ee2cc5624bc545f6840af 100644
--- a/dumux/porousmediumflow/1p/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/1p/implicit/volumevariables.hh
@@ -25,7 +25,7 @@
 #define DUMUX_1P_VOLUME_VARIABLES_HH
 
 #include "properties.hh"
-#include <dumux/implicit/volumevariables.hh>
+#include <dumux/discretization/volumevariables.hh>
 
 #include <dumux/material/fluidstates/immiscible.hh>
 
diff --git a/dumux/porousmediumflow/2p/implicit/volumevariables.hh b/dumux/porousmediumflow/2p/implicit/volumevariables.hh
index 1f698eb780b0ee308a427adf0eba403e1e6febf7..5e001fc651a23bd9101f8c98eda7099ca7fe5999 100644
--- a/dumux/porousmediumflow/2p/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2p/implicit/volumevariables.hh
@@ -27,7 +27,7 @@
 
 #include "properties.hh"
 
-#include <dumux/implicit/volumevariables.hh>
+#include <dumux/discretization/volumevariables.hh>
 
 #include <dune/common/fvector.hh>
 
diff --git a/dumux/porousmediumflow/CMakeLists.txt b/dumux/porousmediumflow/CMakeLists.txt
index e48629a5b22547c44f1c6e3b59cba7dce77561e7..de81dec4ffd6e92f1655c5556218d8115d01c1ba 100644
--- a/dumux/porousmediumflow/CMakeLists.txt
+++ b/dumux/porousmediumflow/CMakeLists.txt
@@ -11,10 +11,9 @@ add_subdirectory("3p")
 add_subdirectory("3p3c")
 add_subdirectory("co2")
 add_subdirectory("compositional")
-add_subdirectory("constitutivelaws")
 add_subdirectory("immiscible")
 add_subdirectory("implicit")
 add_subdirectory("mpnc")
 add_subdirectory("nonisothermal")
 add_subdirectory("richards")
-add_subdirectory("sequential")
\ No newline at end of file
+add_subdirectory("sequential")
diff --git a/dumux/porousmediumflow/constitutivelaws/CMakeLists.txt b/dumux/porousmediumflow/constitutivelaws/CMakeLists.txt
deleted file mode 100644
index e5e8f670beb046f6464be5725217a3c84cf147cd..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/constitutivelaws/CMakeLists.txt
+++ /dev/null
@@ -1,5 +0,0 @@
-#install headers
-install(FILES
-darcyslaw.hh
-fickslaw.hh
-DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/constitutivelaws)
\ No newline at end of file
diff --git a/dumux/porousmediumflow/constitutivelaws/fickslaw.hh b/dumux/porousmediumflow/constitutivelaws/fickslaw.hh
deleted file mode 100644
index a1829d1898416bee0bd9a48ac5c304bb14d78fbf..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/constitutivelaws/fickslaw.hh
+++ /dev/null
@@ -1,379 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program is free software: you can redistribute it and/or modify    *
- *   it under the terms of the GNU General Public License as published by    *
- *   the Free Software Foundation, either version 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program is distributed in the hope that it will be useful,         *
- *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
- *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
- *   GNU General Public License for more details.                            *
- *                                                                           *
- *   You should have received a copy of the GNU General Public License       *
- *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- * \brief This file contains the data which is required to calculate
- *        diffusive mass fluxes due to molecular diffusion with Fick's law.
- */
-#ifndef DUMUX_POROUSMEDIUMFLOW_FICKS_LAW_HH
-#define DUMUX_POROUSMEDIUMFLOW_FICKS_LAW_HH
-
-#include <dune/common/float_cmp.hh>
-
-#include <dumux/common/math.hh>
-#include <dumux/common/parameters.hh>
-
-#include <dumux/implicit/properties.hh>
-
-
-namespace Dumux
-{
-
-namespace Properties
-{
-// forward declaration of properties
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(FluidSystem);
-NEW_PROP_TAG(EffectiveDiffusivityModel);
-}
-
-/*!
- * \ingroup CCTpfaFicksLaw
- * \brief Evaluates the diffusive mass flux according to Fick's law
- */
-template <class TypeTag, typename DiscretizationMethod = void>
-class FicksLaw
-{};
-
-// Specialization for the CC-Tpfa Method
-template <class TypeTag>
-class FicksLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::CCTpfa>::type >
-{
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-    typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffDiffModel;
-    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
-    typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GridView::IndexSet::IndexType IndexType;
-    typedef typename std::vector<IndexType> Stencil;
-
-    using Element = typename GridView::template Codim<0>::Entity;
-
-    enum { dim = GridView::dimension} ;
-    enum { dimWorld = GridView::dimensionworld} ;
-    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
-
-    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
-    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
-public:
-
-    static Scalar flux(const Problem& problem,
-                       const Element& element,
-                       const SubControlVolumeFace& scvFace,
-                       const int phaseIdx,
-                       const int compIdx)
-    {
-        // diffusion tensors are always solution dependent
-        Scalar tij = calculateTransmissibility_(problem, scvFace, phaseIdx, compIdx);
-
-        // Get the inside volume variables
-        const auto insideScvIdx = scvFace.insideScvIdx();
-        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
-        const auto& insideVolVars = problem.model().curVolVars(insideScv);
-
-        // and the outside volume variables
-        const auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx());
-
-        // compute the diffusive flux
-        const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
-        const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
-        const auto rho = 0.5*(insideVolVars.molarDensity(phaseIdx) + outsideVolVars.molarDensity(phaseIdx));
-
-        return rho*tij*(xInside - xOutside);
-    }
-
-    static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace)
-    {
-        std::vector<IndexType> stencil;
-        stencil.clear();
-        if (!scvFace.boundary())
-        {
-            stencil.push_back(scvFace.insideScvIdx());
-            stencil.push_back(scvFace.outsideScvIdx());
-        }
-        else
-            stencil.push_back(scvFace.insideScvIdx());
-
-        return stencil;
-    }
-
-private:
-
-
-    static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace, const int phaseIdx, const int compIdx)
-    {
-        Scalar tij;
-
-        const auto insideScvIdx = scvFace.insideScvIdx();
-        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
-        const auto& insideVolVars = problem.model().curVolVars(insideScvIdx);
-
-        auto insideD = insideVolVars.diffusionCoefficient(phaseIdx, compIdx);
-        insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), insideD);
-        Scalar ti = calculateOmega_(problem, scvFace, insideD, insideScv);
-
-        if (!scvFace.boundary())
-        {
-            const auto outsideScvIdx = scvFace.outsideScvIdx();
-            const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
-            const auto& outsideVolVars = problem.model().curVolVars(outsideScvIdx);
-
-            auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx, compIdx);
-            outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), outsideD);
-            Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideD, outsideScv);
-
-            tij = scvFace.area()*(ti * tj)/(ti + tj);
-        }
-        else
-        {
-            tij = scvFace.area()*ti;
-        }
-
-        return tij;
-    }
-
-    static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &D, const SubControlVolume &scv)
-    {
-        GlobalPosition Dnormal;
-        D.mv(scvFace.unitOuterNormal(), Dnormal);
-
-        auto distanceVector = scvFace.center();
-        distanceVector -= scv.center();
-        distanceVector /= distanceVector.two_norm2();
-
-        Scalar omega = Dnormal * distanceVector;
-        omega *= problem.model().curVolVars(scv).extrusionFactor();
-
-        return omega;
-    }
-
-    static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, Scalar D, const SubControlVolume &scv)
-    {
-        auto distanceVector = scvFace.center();
-        distanceVector -= scv.center();
-        distanceVector /= distanceVector.two_norm2();
-
-        Scalar omega = D * (distanceVector * scvFace.unitOuterNormal());
-        omega *= problem.model().curVolVars(scv).extrusionFactor();
-
-        return omega;
-    }
-};
-
-// Specialization for the Box Method
-template <class TypeTag>
-class FicksLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::Box>::type >
-{
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
-    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
-    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
-    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-    using IndexType = typename GridView::IndexSet::IndexType;
-    using Element = typename GridView::template Codim<0>::Entity;
-
-    enum { dim = GridView::dimension} ;
-    enum { dimWorld = GridView::dimensionworld} ;
-    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
-
-    using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
-    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
-
-public:
-
-    void update(const Problem &problem,
-                const Element& element,
-                const SubControlVolumeFace &scvFace,
-                int phaseIdx, int compIdx)
-    {
-        DUNE_THROW(Dune::NotImplemented, "Fick's law for the Box method is not yet implemented!");
-
-        problemPtr_ = &problem;
-        scvFacePtr_ = &scvFace;
-
-        phaseIdx_ = phaseIdx;
-        compIdx_ = compIdx;
-
-        updateStencil_();
-
-        // TODO for non solution dependent diffusion tensors...
-    }
-
-    void update(const Problem &problem, const SubControlVolumeFace &scvFace,
-                int phaseIdx, int compIdx,
-                VolumeVariables* boundaryVolVars)
-    {
-        boundaryVolVars_ = boundaryVolVars;
-        update(problem, scvFace, phaseIdx, compIdx);
-    }
-
-    void beginFluxComputation(bool boundaryVolVarsUpdated = false)
-    {
-        // diffusion tensors are always solution dependent
-        updateTransmissibilities_();
-
-        // Get the inside volume variables
-        const auto insideScvIdx = scvFace_().insideScvIdx();
-        const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx);
-        const auto* insideVolVars = &problem_().model().curVolVars(insideScv);
-
-        // and the outside volume variables
-        const VolumeVariables* outsideVolVars;
-        if (!scvFace_().boundary())
-            outsideVolVars = &problem_().model().curVolVars(scvFace_().outsideScvIdx());
-        else
-        {
-            outsideVolVars = boundaryVolVars_;
-            if (!boundaryVolVarsUpdated)
-            {
-                // update the boudary volvars for Dirichlet boundaries
-                const auto element = problem_().model().fvGeometries().element(insideScv);
-                const auto dirichletPriVars = problem_().dirichlet(element, scvFace_());
-                boundaryVolVars_->update(dirichletPriVars, problem_(), element, insideScv);
-            }
-        }
-
-        // compute the diffusive flux
-        const auto xInside = insideVolVars->moleFraction(phaseIdx_, compIdx_);
-        const auto xOutside = outsideVolVars->moleFraction(phaseIdx_, compIdx_);
-        const auto rho = 0.5*(insideVolVars->molarDensity(phaseIdx_) + outsideVolVars->molarDensity(phaseIdx_));
-
-        rhoDGradXNormal_ = rho*tij_*(xInside - xOutside);
-    }
-
-
-
-    /*!
-     * \brief A function to calculate the mass flux over a sub control volume face
-     *
-     * \param phaseIdx The index of the phase of which the flux is to be calculated
-     * \param compIdx The index of the transported component
-     */
-    Scalar flux() const
-    {
-        return rhoDGradXNormal_;
-    }
-
-    std::set<IndexType> stencil() const
-    {
-        return stencil_;
-    }
-
-protected:
-
-
-    void updateTransmissibilities_()
-    {
-        const auto insideScvIdx = scvFace_().insideScvIdx();
-        const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx);
-        const auto& insideVolVars = problem_().model().curVolVars(insideScvIdx);
-
-        auto insideD = insideVolVars.diffusionCoefficient(phaseIdx_, compIdx_);
-        insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx_), insideD);
-        Scalar ti = calculateOmega_(insideD, insideScv);
-
-        if (!scvFace_().boundary())
-        {
-            const auto outsideScvIdx = scvFace_().outsideScvIdx();
-            const auto& outsideScv = problem_().model().fvGeometries().subControlVolume(outsideScvIdx);
-            const auto& outsideVolVars = problem_().model().curVolVars(outsideScvIdx);
-
-            auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx_, compIdx_);
-            outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx_), outsideD);
-            Scalar tj = -1.0*calculateOmega_(outsideD, outsideScv);
-
-            tij_ = scvFace_().area()*(ti * tj)/(ti + tj);
-        }
-        else
-        {
-            tij_ = scvFace_().area()*ti;
-        }
-    }
-
-    void updateStencil_()
-    {
-        // fill the stencil
-        if (!scvFace_().boundary())
-            stencil_= {scvFace_().insideScvIdx(), scvFace_().outsideScvIdx()};
-        else
-            // fill the stencil
-            stencil_ = {scvFace_().insideScvIdx()};
-    }
-
-    Scalar calculateOmega_(const DimWorldMatrix &D, const SubControlVolume &scv) const
-    {
-        GlobalPosition Dnormal;
-        D.mv(scvFace_().unitOuterNormal(), Dnormal);
-
-        auto distanceVector = scvFace_().center();
-        distanceVector -= scv.center();
-        distanceVector /= distanceVector.two_norm2();
-
-        Scalar omega = Dnormal * distanceVector;
-        omega *= problem_().model().curVolVars(scv).extrusionFactor();
-
-        return omega;
-    }
-
-    Scalar calculateOmega_(Scalar D, const SubControlVolume &scv) const
-    {
-        auto distanceVector = scvFace_().center();
-        distanceVector -= scv.center();
-        distanceVector /= distanceVector.two_norm2();
-
-        Scalar omega = D * (distanceVector * scvFace_().unitOuterNormal());
-        omega *= problem_().model().curVolVars(scv).extrusionFactor();
-
-        return omega;
-    }
-
-    const Problem &problem_() const
-    {
-        return *problemPtr_;
-    }
-
-    const SubControlVolumeFace& scvFace_() const
-    {
-        return *scvFacePtr_;
-    }
-
-    const Problem *problemPtr_;
-    const SubControlVolumeFace *scvFacePtr_; //!< Pointer to the sub control volume face for which the flux variables are created
-    std::set<IndexType> stencil_;         //!< Indices of the cells of which the pressure is needed for the flux calculation
-
-    //! Boundary volume variables (they only get updated on Dirichlet boundaries)
-    VolumeVariables* boundaryVolVars_;
-
-    IndexType phaseIdx_;
-    IndexType compIdx_;
-    //! Precomputed values
-    Scalar tij_; //!< transmissibility for the flux calculation
-    Scalar rhoDGradXNormal_; //! rho*D(grad(x))*n
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/porousmediumflow/implicit/fluxvariables.hh b/dumux/porousmediumflow/implicit/fluxvariables.hh
index 3038322ea0d75f3b6faec5cbf48c55d8a4687eb7..a528cd9872120c7b479e97c192de8befafdb2fbd 100644
--- a/dumux/porousmediumflow/implicit/fluxvariables.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariables.hh
@@ -24,7 +24,7 @@
 #define DUMUX_POROUSMEDIUMFLOW_IMPLICIT_FLUXVARIABLES_HH
 
 #include <dumux/implicit/properties.hh>
-#include <dumux/implicit/fluxvariablesbase.hh>
+#include <dumux/discretization/fluxvariablesbase.hh>
 
 namespace Dumux
 {