diff --git a/CHANGELOG.md b/CHANGELOG.md
index b59f48585619847c4f9f0362845d166288cc2cc0..dc40f538db5e89dc51815e2d2466c1a615666750 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,10 @@ Differences Between DuMu<sup>x</sup> 3.3 and DuMu<sup>x</sup> 3.2
 - The DuMu<sup>x</sup> install script has been translated to Python to improve portability. The old shell script will be removed after release 3.3.
 - The velocity reconstruction for immiscible porous-media models has been improved, leading to slightly
   different velocity fields in the vicinity of Neumann boundaries.
+- Basic support for Python bindings has been added. Python bindings are an experimental feature
+  and might undergo unannounced API changes until further notice. This concerns the files in the folders `python` and `dumux/python`. To activate
+    - add `-DDUNE_ENABLE_PYTHONBINDINGS=TRUE` and `-DCMAKE_POSITION_INDEPENDENT_CODE=TRUE` to your CMAKE_FLAGS and run dunecontrol
+    - run `python3 dune-common/bin/setup_dunepy.py`
 
 ### Immediate interface changes not allowing/requiring a deprecation period:
 
@@ -20,8 +24,8 @@ Differences Between DuMu<sup>x</sup> 3.3 and DuMu<sup>x</sup> 3.2
 ### Deleted classes/files, property names, constants/enums:
 
 - Everything that has been deprecated before release 3.2 has been removed.
-- All of the geometry headers previously saved in `dumux/common/geometry` have been relocated to `dumux/geometry`. 
-  The headers in `dumux/common/geometry` are deprecated and will be removed in 3.3. The geometry tests have been moved from `test/common/geometry` 
+- All of the geometry headers previously saved in `dumux/common/geometry` have been relocated to `dumux/geometry`.
+  The headers in `dumux/common/geometry` are deprecated and will be removed in 3.3. The geometry tests have been moved from `test/common/geometry`
   and `test/common/boundingboxtree` to `test/geometry`.
 
 ### Other noteworthy changes:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8ad564c6feb74967dc13aabe4a82703392c234d7..a0e184b4066da101b01d324eb2ab4fba455d84f4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -33,5 +33,11 @@ add_subdirectory(dumux)
 add_subdirectory(test EXCLUDE_FROM_ALL)
 add_subdirectory(examples EXCLUDE_FROM_ALL)
 
+# if Python bindings are enabled, include necessary sub directories.
+if(DUNE_ENABLE_PYTHONBINDINGS)
+  add_subdirectory(python)
+  dune_python_install_package(PATH "python")
+endif()
+
 # finalize the dune project, e.g. generating config.h etc.
 finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
diff --git a/cmake.opts b/cmake.opts
index 8c5f184a1b0c7969bf624461a8e8247dd652d9ee..de36f45777d11d9586a27a754559cebd02bf77a5 100644
--- a/cmake.opts
+++ b/cmake.opts
@@ -46,6 +46,10 @@ DUMUX_ENABLE_HEADERCHECK=OFF
 DUMUX_DISABLE_PROP_MACROS=""
 #DUMUX_DISABLE_PROP_MACROS="-DDUMUX_ENABLE_OLD_PROPERTY_MACROS=0"
 
+# to enable Python bindings
+DUMUX_ENABLE_PYTHON_BINDINGS=""
+#DUMUX_ENABLE_PYTHON_BINDINGS="-DDUNE_ENABLE_PYTHONBINDINGS=1 -DCMAKE_POSITION_INDEPENDENT_CODE=1"
+
 # for debug opts you can set DCMAKE_BUILD_TYPE to "Debug" or "RelWithDebInfo"
 # you can also do this in any of the CMakeLists.txt in Dumux
 # just rerun cmake again afterwards (run cmake <path-to-build-dir>)
@@ -56,4 +60,5 @@ CMAKE_FLAGS="$SPECIFIC_COMPILER $SPECIFIC_GENERATOR $OPM_FLAGS
 -DCMAKE_CXX_FLAGS_RELWITHDEBINFO='$GXX_RELEASE_OPTS $GXX_RELEASE_WARNING_OPTS -g -ggdb -Wall $DUMUX_DISABLE_PROP_MACROS'
 -DCMAKE_BUILD_TYPE=Release
 -DENABLE_HEADERCHECK=$DUMUX_ENABLE_HEADERCHECK
+$DUMUX_ENABLE_PYTHON_BINDINGS
 "
diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt
index 819d019ad5f37d819bd128854c2892abba2684c9..8b27a751738d8a46aa6be7c793cfafad433a5d6d 100644
--- a/dumux/CMakeLists.txt
+++ b/dumux/CMakeLists.txt
@@ -13,3 +13,8 @@ add_subdirectory(multidomain)
 add_subdirectory(nonlinear)
 add_subdirectory(parallel)
 add_subdirectory(porousmediumflow)
+
+# if Python bindings are enabled, include necessary sub directories.
+if(DUNE_ENABLE_PYTHONBINDINGS)
+  add_subdirectory(python)
+endif()
diff --git a/dumux/python/CMakeLists.txt b/dumux/python/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4f95771a6614eff8b2224375adc22822a3c0c1ec
--- /dev/null
+++ b/dumux/python/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_subdirectory(common)
+add_subdirectory(discretization)
diff --git a/dumux/python/common/CMakeLists.txt b/dumux/python/common/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..19d092e5cc5cdd6502a833d3f00f9654cf6c82f2
--- /dev/null
+++ b/dumux/python/common/CMakeLists.txt
@@ -0,0 +1,5 @@
+install(FILES
+boundarytypes.hh
+fvproblem.hh
+timeloop.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/common)
diff --git a/dumux/python/common/boundarytypes.hh b/dumux/python/common/boundarytypes.hh
new file mode 100644
index 0000000000000000000000000000000000000000..460ca8a155352fb7e678bcef254b92fc7f91a1b2
--- /dev/null
+++ b/dumux/python/common/boundarytypes.hh
@@ -0,0 +1,65 @@
+// -*- 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 3 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 TODO: docme!
+ */
+
+#ifndef DUMUX_PYTHON_COMMON_BOUNDARYTYPES_HH
+#define DUMUX_PYTHON_COMMON_BOUNDARYTYPES_HH
+
+#include <dune/common/classname.hh>
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/common/typeregistry.hh>
+
+#include <dumux/common/boundarytypes.hh>
+
+namespace Dumux::Python {
+
+template <class BoundaryTypes, class... Options>
+void registerBoundaryTypes(pybind11::handle scope, pybind11::class_<BoundaryTypes, Options...> cls)
+{
+    using pybind11::operator""_a;
+
+    cls.def(pybind11::init());
+    cls.def("reset", &BoundaryTypes::reset);
+    cls.def("setNeumann", &BoundaryTypes::setAllNeumann);
+    cls.def("setDirichlet", &BoundaryTypes::setAllDirichlet);
+    cls.def("isDirichlet", &BoundaryTypes::hasDirichlet);
+    cls.def("isNeumann", &BoundaryTypes::hasNeumann);
+}
+
+template <class BoundaryTypes>
+void registerBoundaryTypes(pybind11::handle scope)
+{
+    using namespace Dune::Python;
+
+    auto [cls, addedToRegistry] = insertClass<BoundaryTypes>(
+        scope, "BoundaryTypes",
+        GenerateTypeName(Dune::className<BoundaryTypes>()),
+        IncludeFiles{"dumux/python/common/boundarytypes.hh"}
+    );
+
+    if (addedToRegistry)
+        registerBoundaryTypes(scope, cls);
+}
+
+} // namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/common/fvproblem.hh b/dumux/python/common/fvproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bf8aaba7bf0d6af994d1c539594a730d3afed192
--- /dev/null
+++ b/dumux/python/common/fvproblem.hh
@@ -0,0 +1,195 @@
+// -*- 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 3 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 TODO: docme!
+ */
+
+#ifndef DUMUX_PYTHON_COMMON_FVPROBLEM_HH
+#define DUMUX_PYTHON_COMMON_FVPROBLEM_HH
+
+#include <string>
+#include <memory>
+#include <tuple>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/python/pybind11/pybind11.h>
+
+#include <dumux/common/boundarytypes.hh>
+#include <dumux/python/common/boundarytypes.hh>
+
+namespace Dumux::Python {
+
+/*!
+ * \ingroup Common
+ * \brief A C++ wrapper for a Python problem
+ */
+template<class GridGeometry_, class PrimaryVariables>
+class FVProblem
+{
+public:
+    using GridGeometry = GridGeometry_;
+    using Scalar = typename PrimaryVariables::value_type;
+    using NumEqVector = Dune::FieldVector<Scalar, PrimaryVariables::dimension>;
+    using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+    using FVElementGeometry = typename GridGeometry::LocalView;
+    using SubControlVolume = typename GridGeometry::SubControlVolume;
+    using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+    static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
+    static constexpr std::size_t numEq = static_cast<std::size_t>(PrimaryVariables::dimension);
+    using BoundaryTypes = Dumux::BoundaryTypes<PrimaryVariables::dimension>;
+
+    FVProblem(std::shared_ptr<const GridGeometry> gridGeometry, pybind11::object pyProblem)
+    : gridGeometry_(gridGeometry), pyProblem_(pyProblem)
+    {}
+
+    std::string name() const
+    {
+        return pyProblem_.attr("name").template cast<std::string>();
+    }
+
+    BoundaryTypes boundaryTypes(const Element &element,
+                                const SubControlVolume &scv) const
+    {
+        if constexpr (!isBox)
+            DUNE_THROW(Dune::InvalidStateException, "boundaryTypes(..., scv) called for cell-centered method.");
+        else
+            return pyProblem_.attr("boundaryTypes")(element, scv).template cast<BoundaryTypes>();
+    }
+
+    BoundaryTypes boundaryTypes(const Element &element,
+                                const SubControlVolumeFace &scvf) const
+    {
+        if constexpr (isBox)
+            DUNE_THROW(Dune::InvalidStateException, "boundaryTypes(..., scvf) called for box method.");
+        else
+            return pyProblem_.attr("boundaryTypes")(element, scvf).template cast<BoundaryTypes>();
+    }
+
+    PrimaryVariables dirichlet(const Element &element,
+                               const SubControlVolume &scv) const
+    {
+        if constexpr (!isBox)
+            DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for cell-centered method.");
+        else
+            return pyProblem_.attr("dirichlet")(element, scv).template cast<PrimaryVariables>();
+    }
+
+    PrimaryVariables dirichlet(const Element &element,
+                               const SubControlVolumeFace &scvf) const
+    {
+        if constexpr (isBox)
+            DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for box method.");
+        else
+            return pyProblem_.attr("dirichlet")(element, scvf).template cast<PrimaryVariables>();
+    }
+
+    template<class ElementVolumeVariables, class ElementFluxVariablesCache>
+    NumEqVector neumann(const Element& element,
+                        const FVElementGeometry& fvGeometry,
+                        const ElementVolumeVariables& elemVolVars,
+                        const ElementFluxVariablesCache& elemFluxVarsCache,
+                        const SubControlVolumeFace& scvf) const
+    {
+        return pyProblem_.attr("neumann")(element, fvGeometry, scvf).template cast<NumEqVector>();
+    }
+
+    template<class ElementVolumeVariables>
+    NumEqVector source(const Element &element,
+                       const FVElementGeometry& fvGeometry,
+                       const ElementVolumeVariables& elemVolVars,
+                       const SubControlVolume &scv) const
+    {
+        return pyProblem_.attr("source")(element, fvGeometry, scv).template cast<NumEqVector>();
+    }
+
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    {
+        return pyProblem_.attr("sourceAtPos")(globalPos).template cast<NumEqVector>();
+    }
+
+    template<class Entity>
+    PrimaryVariables initial(const Entity& entity) const
+    {
+        return pyProblem_.attr("intial")(entity).template cast<PrimaryVariables>();
+    }
+
+    template<class ElementSolution>
+    Scalar extrusionFactor(const Element& element,
+                           const SubControlVolume& scv,
+                           const ElementSolution& elemSol) const
+    {
+        return pyProblem_.attr("extrusionFactor")(element, scv).template cast<Scalar>();
+    }
+
+    const GridGeometry& gridGeometry() const
+    { return *gridGeometry_; }
+
+private:
+    std::shared_ptr<const GridGeometry> gridGeometry_;
+    pybind11::object pyProblem_;
+};
+
+// Python wrapper for the above FVProblem C++ class
+template<class Problem, class... options>
+void registerFVProblem(pybind11::handle scope, pybind11::class_<Problem, options...> cls)
+{
+    using pybind11::operator""_a;
+    using namespace Dune::Python;
+
+    using GridGeometry = typename Problem::GridGeometry;
+    cls.def(pybind11::init([](std::shared_ptr<const GridGeometry> gridGeometry, pybind11::object p){
+        return std::make_shared<Problem>(gridGeometry, p);
+    }));
+
+    cls.def_property_readonly("name", &Problem::name);
+    cls.def_property_readonly("numEq", [](Problem&){ return Problem::numEq; });
+
+    using GridView = typename GridGeometry::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using Vertex = typename GridView::template Codim<GridView::dimension>::Entity;
+
+    if constexpr (Problem::isBox)
+    {
+        using SCV = typename Problem::SubControlVolume;
+        cls.def("boundaryTypes", pybind11::overload_cast<const Element&, const SCV&>(&Problem::boundaryTypes, pybind11::const_), "element"_a, "scv"_a);
+        cls.def("dirichlet", pybind11::overload_cast<const Element&, const SCV&>(&Problem::dirichlet, pybind11::const_), "element"_a, "scv"_a);
+    }
+    else
+    {
+        using SCVF = typename Problem::SubControlVolumeFace;
+        cls.def("boundaryTypes", pybind11::overload_cast<const Element&, const SCVF&>(&Problem::boundaryTypes, pybind11::const_), "element"_a, "scvf"_a);
+        cls.def("dirichlet", pybind11::overload_cast<const Element&, const SCVF&>(&Problem::dirichlet, pybind11::const_), "element"_a, "scvf"_a);
+    }
+
+    cls.def("neumann", &Problem::template neumann<decltype(std::ignore), decltype(std::ignore)>);
+    cls.def("source", &Problem::template source<decltype(std::ignore)>);
+    cls.def("sourceAtPos", &Problem::sourceAtPos);
+    cls.def("initial", &Problem::template initial<Element>);
+    cls.def("initial", &Problem::template initial<Vertex>);
+    cls.def("extrusionFactor", &Problem::template extrusionFactor<decltype(std::ignore)>);
+    cls.def("gridGeometry", &Problem::gridGeometry);
+}
+
+} // end namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/common/timeloop.hh b/dumux/python/common/timeloop.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9117c0ee48fd13557549fc2eb1385e3a407ebbb0
--- /dev/null
+++ b/dumux/python/common/timeloop.hh
@@ -0,0 +1,71 @@
+// -*- 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 3 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 TODO: docme!
+ */
+
+#ifndef DUMUX_PYTHON_COMMON_TIMELOOP_HH
+#define DUMUX_PYTHON_COMMON_TIMELOOP_HH
+
+#include <dumux/common/timeloop.hh>
+
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/pybind11/stl.h>
+
+namespace Dumux::Python {
+
+template <class Scalar, class... options>
+void registerTimeLoop(pybind11::handle scope,
+                      pybind11::class_<CheckPointTimeLoop<Scalar>, options...> cls)
+{
+    using pybind11::operator""_a;
+
+    using TimeLoop = CheckPointTimeLoop<Scalar>;
+    cls.def(pybind11::init([](Scalar startTime, Scalar dt, Scalar endTime, bool verbose){
+        return new TimeLoop(startTime, dt, endTime, verbose);
+    }), "startTime"_a, "dt"_a, "endTime"_a, "verbose"_a=true);
+
+    cls.def("time", &TimeLoop::time);
+    cls.def("timeStepSize", &TimeLoop::timeStepSize);
+    cls.def("start", &TimeLoop::start);
+    cls.def("reset", &TimeLoop::reset, "startTime"_a, "dt"_a, "endTime"_a, "verbose"_a=true);
+    cls.def("advanceTimeStep", &TimeLoop::advanceTimeStep);
+    cls.def("setTimeStepSize", &TimeLoop::setTimeStepSize, "dt"_a);
+    cls.def("setMaxTimeStepSize", &TimeLoop::setMaxTimeStepSize, "dt"_a);
+    cls.def("finished", &TimeLoop::finished);
+    cls.def("reportTimeStep", &TimeLoop::reportTimeStep);
+    cls.def("finalize", [](TimeLoop& self){ self.finalize(); });
+    cls.def("setPeriodicCheckPoint", &TimeLoop::setPeriodicCheckPoint, "interval"_a, "offset"_a=0.0);
+    cls.def("isCheckPoint", &TimeLoop::isCheckPoint);
+    cls.def("setCheckPoints", [](TimeLoop& self, const std::vector<double>& checkPoints) {
+        self.setCheckPoint(checkPoints);
+    });
+}
+
+template<class Scalar>
+void registerTimeLoop(pybind11::handle scope, const char *clsName = "TimeLoop")
+{
+    pybind11::class_<CheckPointTimeLoop<Scalar>> cls(scope, clsName);
+    registerTimeLoop(scope, cls);
+}
+
+} // namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/discretization/CMakeLists.txt b/dumux/python/discretization/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..efe5f30fb1b4cb2ddd6ba7a4dd6e3775ea611bc3
--- /dev/null
+++ b/dumux/python/discretization/CMakeLists.txt
@@ -0,0 +1,3 @@
+install(FILES
+gridgeometry.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/discretization)
diff --git a/dumux/python/discretization/gridgeometry.hh b/dumux/python/discretization/gridgeometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..69c6af1f09965ab30f07b0b893e938499bccb6f9
--- /dev/null
+++ b/dumux/python/discretization/gridgeometry.hh
@@ -0,0 +1,151 @@
+// -*- 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 3 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 TODO: docme!
+ */
+
+#ifndef DUMUX_PYTHON_DISCRETIZATION_GRIDGEOMETRY_HH
+#define DUMUX_PYTHON_DISCRETIZATION_GRIDGEOMETRY_HH
+
+#include <memory>
+#include <dune/common/classname.hh>
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/common/typeregistry.hh>
+
+namespace Dumux::Python {
+
+namespace Impl {
+
+template<class SCV>
+void registerSubControlVolume(pybind11::handle scope)
+{
+    using namespace Dune::Python;
+
+    auto [cls, addedToRegistry] = insertClass<SCV>(
+        scope, "SubControlVolume",
+        GenerateTypeName(Dune::className<SCV>()),
+        IncludeFiles{"dumux/python/discretization/gridgeometry.hh"}
+    );
+
+    if (addedToRegistry)
+    {
+        using pybind11::operator""_a;
+
+        cls.def("center", &SCV::center);
+        cls.def("volume", &SCV::volume);
+        cls.def("dofIndex", &SCV::dofIndex);
+        cls.def("localDofIndex", &SCV::localDofIndex);
+        cls.def("dofPosition", &SCV::dofPosition);
+        cls.def("elementIndex", &SCV::elementIndex);
+
+    }
+}
+
+template<class SCVF>
+void registerSubControlVolumeFace(pybind11::handle scope)
+{
+    using namespace Dune::Python;
+
+    auto [cls, addedToRegistry] = insertClass<SCVF>(
+        scope, "SubControlVolumeFace",
+        GenerateTypeName(Dune::className<SCVF>()),
+        IncludeFiles{"dumux/python/discretization/gridgeometry.hh"}
+    );
+
+    if (addedToRegistry)
+    {
+        using pybind11::operator""_a;
+
+        cls.def("center", &SCVF::center);
+        cls.def("area", &SCVF::area);
+        cls.def("ipGlobal", &SCVF::ipGlobal);
+        cls.def("boundary", &SCVF::boundary);
+        cls.def("unitOuterNormal", &SCVF::unitOuterNormal);
+        cls.def("insideScvIdx", &SCVF::insideScvIdx);
+        cls.def("outsideScvIdx", [](SCVF& self){ return self.outsideScvIdx(); });
+        cls.def("index", &SCVF::index);
+    }
+}
+
+template<class FVEG>
+void registerFVElementGeometry(pybind11::handle scope)
+{
+    using namespace Dune::Python;
+
+    auto [cls, addedToRegistry] = insertClass<FVEG>(
+        scope, "FVElementGeometry",
+        GenerateTypeName(Dune::className<FVEG>()),
+        IncludeFiles{"dumux/python/discretization/gridgeometry.hh"}
+    );
+
+    if (addedToRegistry)
+    {
+        using pybind11::operator""_a;
+
+        cls.def("numScvf", &FVEG::numScv);
+        cls.def("numScv", &FVEG::numScv);
+        cls.def("bind", &FVEG::bind, "element"_a);
+        cls.def("hasBoundaryScvf", &FVEG::hasBoundaryScvf);
+        cls.def("scvs", [](FVEG& self){
+            const auto range = scvs(self);
+            return pybind11::make_iterator(range.begin(), range.end());
+        }, pybind11::keep_alive<0, 1>());
+        cls.def("scvfs", [](FVEG& self){
+            const auto range = scvfs(self);
+            return pybind11::make_iterator(range.begin(), range.end());
+        }, pybind11::keep_alive<0, 1>());
+    }
+}
+
+} // end namespace Impl
+
+// see python/dumux/discretization/__init__.py for how this is used for JIT compilation
+template <class GG, class... Options>
+void registerGridGeometry(pybind11::handle scope, pybind11::class_<GG, Options...> cls)
+{
+    using pybind11::operator""_a;
+
+    using GridView = typename GG::GridView;
+
+    cls.def(pybind11::init([](GridView gridView){
+        return std::make_shared<GG>(gridView);
+    }), "gridView"_a);
+
+    cls.def("update", &GG::update);
+    cls.def("numDofs", &GG::numDofs);
+    cls.def("numScv", &GG::numScv);
+    cls.def("numScvf", &GG::numScvf);
+
+    using SubControlVolume = typename GG::SubControlVolume;
+    Impl::registerSubControlVolume<SubControlVolume>(scope);
+
+    using SubControlVolumeFace = typename GG::SubControlVolumeFace;
+    Impl::registerSubControlVolumeFace<SubControlVolumeFace>(scope);
+
+    // also compile the corresponding local view
+    using LocalView = typename GG::LocalView;
+    Impl::registerFVElementGeometry<LocalView>(scope);
+    // and make it accessible
+    cls.def("localView", [](GG& self){ return localView(self); });
+}
+
+} // namespace Dumux::Python
+
+#endif
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0c4ddccf90e05122277901c3cfbf4c9b2957cb68
--- /dev/null
+++ b/python/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_subdirectory(dumux)
+configure_file(setup.py.in setup.py)
diff --git a/python/README.md b/python/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..851ef4e6fb803e4b9a424e1666b9f782935c8181
--- /dev/null
+++ b/python/README.md
@@ -0,0 +1,71 @@
+# Python bindings for DuMu<sup>x</sup>
+
+The Python bindings have two objectives
+* making it possible to use DuMu<sup>x</sup> from Python
+* use Python code in DuMu<sup>x</sup>
+
+The binding are experimental until further notice which means
+the API might undergo unannounced changes breaking backwards
+compatibility. Track changes by regularly checking for current merge requests
+and newly merged Dumux commits on GitLab. Python bindings require Python 3.
+
+##  Installation
+
+Python bindings are part of the dune core modules >= version 2.7.
+We recommend to use the Dune modules on the master branch to try out
+Python bindings.
+
+### Example development setup
+
+This section is rather new and experimental. Please help improving
+this documentation in the future.
+
+Checkout the `master` branch of the Dune core modules and DuMu<sup>x</sup>
+
+```
+git clone https://gitlab.dune-project.org/core/dune-common
+git clone https://gitlab.dune-project.org/core/dune-geometry
+git clone https://gitlab.dune-project.org/core/dune-grid
+git clone https://gitlab.dune-project.org/core/dune-localfunctions
+git clone https://gitlab.dune-project.org/core/dune-istl
+git clone https://git.iws.uni-stuttgart.de/dumux-repositories/dumux
+
+cp dumux/cmake.opts .
+```
+
+Enable the Python bindings by adding the flag in `cmake.opts` (see comments inside the `.opts` file).
+Then run dunecontrol which also builds the Dune Python bindings.
+
+```
+./dune-common/bin/dunecontrol --opts=cmake.opts all
+```
+
+Add the Python binding modules to your Python path like this and install them:
+
+```
+export PYTHONPATH=$(pwd)/dune-common/build-cmake/python:$(pwd)/dune-grid/build-cmake/python:$(pwd)/dune-geometry/build-cmake/python:$(pwd)/dune-istl/build-cmake/python:$(pwd)/dune-localfunctions/build-cmake/python:$(pwd)/dumux/build-cmake/python
+python3 dune-common/bin/setup-dunepy.py --opts=cmake.opts install
+```
+
+If you are getting error with loading MPI in Python you might need to preload the MPI library
+
+```
+export LD_PRELOAD=/usr/lib/x86_64-linux-gnu/openmpi/lib/libmpi.so
+```
+
+Replace the path with the path to the MPI library on your system.
+Run your first DuMu<sup>x</sup> Python test
+
+```
+python3 dumux/test/python/test_py_gridgeometry.py
+```
+
+The Python bindings are based on just-in-time compilation of C++ code,
+so the first execution might take considerably longer than subsequent executions.
+In particular compiling the Dune grid interface may a few minutes.
+
+You can run all currently existing DuMu<sup>x</sup> Python tests with
+```
+cd dumux/build-cmake
+ctest -L python
+```
diff --git a/python/dumux/CMakeLists.txt b/python/dumux/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4f95771a6614eff8b2224375adc22822a3c0c1ec
--- /dev/null
+++ b/python/dumux/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_subdirectory(common)
+add_subdirectory(discretization)
diff --git a/python/dumux/common/CMakeLists.txt b/python/dumux/common/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..38d491fa22b9d7a63e3345c5f40e09514407b43d
--- /dev/null
+++ b/python/dumux/common/CMakeLists.txt
@@ -0,0 +1,5 @@
+add_python_targets(common
+  __init__
+)
+dune_add_pybind11_module(NAME _common)
+set_property(TARGET _common PROPERTY LINK_LIBRARIES dunecommon dunegrid APPEND)
diff --git a/python/dumux/common/__init__.py b/python/dumux/common/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..23cdbfa5710456ff954c9281e27e7df725cb6ff4
--- /dev/null
+++ b/python/dumux/common/__init__.py
@@ -0,0 +1,48 @@
+from ._common import *
+
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+
+# A problem decorator generator for Python problems
+#
+# from dumux.common import FVProblem
+# @FVProblem(gridGeometry)
+# class MyProblem:
+#    ...
+#
+def FVProblem(gridGeometry):
+
+    def createModule(numEq):
+        priVarType = "Dune::FieldVector<double, {}>".format(numEq)
+        ggType = gridGeometry._typeName
+        problemType = "Dumux::Python::FVProblem<{}, {}>".format(ggType, priVarType)
+        includes = gridGeometry._includes + ["dumux/python/common/fvproblem.hh"]
+        moduleName = "fvproblem_" + hashIt(problemType)
+        holderType = "std::shared_ptr<{}>".format(problemType)
+        generator = SimpleGenerator("FVProblem", "Dumux::Python")
+        module = generator.load(includes, problemType, moduleName, options=[holderType])
+        return module
+
+    def FVProblemDecorator(Cls):
+        module = createModule(Cls.numEq)
+        def createFVProblem():
+            return module.FVProblem(gridGeometry, Cls())
+        return createFVProblem
+
+    return FVProblemDecorator
+
+
+# Function for JIT copmilation of Dumux::BoundaryTypes
+def BoundaryTypes(numEq=1):
+    # only copmile this once per numEq
+    cacheKey = "BoundaryTypes_{}".format(numEq)
+    try:
+        return globals()[cacheKey]()
+    except KeyError:
+        includes = ["dumux/python/common/boundarytypes.hh"]
+        typeName = "Dumux::BoundaryTypes<{}>".format(numEq)
+        moduleName = "boundarytypes_" + hashIt(typeName)
+        generator = SimpleGenerator("BoundaryTypes", "Dumux::Python")
+        module = generator.load(includes, typeName, moduleName)
+        globals().update({cacheKey : module.BoundaryTypes})
+    return globals()[cacheKey]()
diff --git a/python/dumux/common/_common.cc b/python/dumux/common/_common.cc
new file mode 100644
index 0000000000000000000000000000000000000000..49bca345331814d7741894631e309edbd3c166f2
--- /dev/null
+++ b/python/dumux/common/_common.cc
@@ -0,0 +1,27 @@
+// -*- 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 3 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/>.   *
+ *****************************************************************************/
+
+#include <dune/python/pybind11/pybind11.h>
+#include <dumux/python/common/timeloop.hh>
+
+PYBIND11_MODULE(_common, module)
+{
+    // export time loop
+    Dumux::Python::registerTimeLoop<double>(module);
+}
diff --git a/python/dumux/discretization/CMakeLists.txt b/python/dumux/discretization/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..22b582c564c01350d456389cbab7543022d79cb3
--- /dev/null
+++ b/python/dumux/discretization/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_python_targets(discretization
+  __init__
+)
diff --git a/python/dumux/discretization/__init__.py b/python/dumux/discretization/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..238240b2ac2d44089298b8e2d159c266f8ac406d
--- /dev/null
+++ b/python/dumux/discretization/__init__.py
@@ -0,0 +1,22 @@
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+
+# construct a GridGeometry from a gridView
+# the grid geometry is JIT compiled
+def GridGeometry(gridView, discMethod="cctpfa"):
+    includes = gridView._includes + ["dumux/python/discretization/gridgeometry.hh"]
+
+    if discMethod == "cctpfa":
+        includes += ["dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh"]
+        typeName = "Dumux::CCTpfaFVGridGeometry<" + gridView._typeName + ">"
+    elif discMethod == "box":
+        includes += ["dumux/discretization/box/fvgridgeometry.hh"]
+        typeName = "Dumux::BoxFVGridGeometry<double, " + gridView._typeName + ">"
+    else:
+        raise ValueError("Unknown discMethod {}".format(discMethod))
+
+    moduleName = "gridgeometry_" + hashIt(typeName)
+    holderType = "std::shared_ptr<{}>".format(typeName)
+    generator = SimpleGenerator("GridGeometry", "Dumux::Python")
+    module = generator.load(includes, typeName, moduleName, options=[holderType])
+    return module.GridGeometry(gridView)
diff --git a/python/setup.py.in b/python/setup.py.in
new file mode 100644
index 0000000000000000000000000000000000000000..8cb1413ea2a7a3be12ab510cdc446af4eb2ea4d7
--- /dev/null
+++ b/python/setup.py.in
@@ -0,0 +1,11 @@
+from setuptools import setup, find_namespace_packages
+
+setup(name="dumux",
+      description="Python lib for dumux",
+      version="${DUMUX_VERSION}",
+      author="Timo Koch",
+      packages = find_namespace_packages(include=['dumux.*']),
+      zip_safe = 0,
+      package_data = {'': ['*.so']},
+      install_requires = ['portalocker'],
+)
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index b6d855fcc0dcabd9bc4072ab6c4d1a76bba629f0..d11f39959f967d9f73e7f47ee592157cf3b30ae4 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -10,3 +10,8 @@ add_subdirectory(multidomain)
 add_subdirectory(nonlinear)
 add_subdirectory(porousmediumflow)
 add_subdirectory(discretization)
+
+# if Python bindings are enabled, include Python binding tests
+if(DUNE_ENABLE_PYTHONBINDINGS)
+  add_subdirectory(python)
+endif()
diff --git a/test/python/CMakeLists.txt b/test/python/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4f50f615d14a2a9c8f0c589cceb456fd5070786b
--- /dev/null
+++ b/test/python/CMakeLists.txt
@@ -0,0 +1,11 @@
+dune_symlink_to_source_files(test_py_gridgeometry.py)
+dune_python_add_test(NAME test_py_gridgeometry
+                     COMMAND ${PYTHON_EXECUTABLE} test_py_gridgeometry.py
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+                     LABELS python unit)
+
+dune_symlink_to_source_files(test_py_fvproblem.py)
+dune_python_add_test(NAME test_py_fvproblem
+                     COMMAND ${PYTHON_EXECUTABLE} test_py_fvproblem.py
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+                     LABELS python unit)
diff --git a/test/python/test_boundaryloop.hh b/test/python/test_boundaryloop.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9642f390b8843bf3859381009e8eb53aa130d40a
--- /dev/null
+++ b/test/python/test_boundaryloop.hh
@@ -0,0 +1,110 @@
+// -*- 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 3 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 TODO: docme!
+ */
+
+#ifndef DUMUX_TEST_PYTHON_FVPROBLEM_BOUNDARY_LOOP_HH
+#define DUMUX_TEST_PYTHON_FVPROBLEM_BOUNDARY_LOOP_HH
+
+#include <string>
+#include <memory>
+#include <tuple>
+
+#include <iostream>
+#include <iomanip>
+#include <dune/python/pybind11/iostream.h>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/python/pybind11/pybind11.h>
+
+#include <dumux/common/boundarytypes.hh>
+#include <dumux/python/common/boundarytypes.hh>
+
+namespace Dumux::Python {
+
+/////////////////////////////////////////////////////////
+/// Some debugging/testing  stuff
+///////////////////////////////////////////////////////////
+template<class Problem>
+void printProblemTest(const Problem& problem)
+{
+    std::size_t numNeumann = 0;
+    std::size_t numDirichlet = 0;
+    double totalSource = 0;
+    const auto& gg = problem.gridGeometry();
+    for (const auto& element : elements(gg.gridView()))
+    {
+        auto fvGeometry = localView(gg);
+        fvGeometry.bindElement(element);
+
+        for (const auto& scv : scvs(fvGeometry))
+        {
+            const auto boundaryTypes = problem.boundaryTypes(element, scv);
+
+            if (boundaryTypes.hasDirichlet())
+                numDirichlet += 1;
+            else if (boundaryTypes.hasNeumann())
+                numNeumann += 1;
+
+            totalSource += problem.sourceAtPos(scv.center())[0]*scv.volume();
+        }
+    }
+
+    std::cout << "[cpp] Found " << numNeumann << " Neumann faces and " << numDirichlet << " Dirichlet faces" << std::endl;
+    std::cout << "[cpp] Total source " << totalSource << std::fixed << std::setprecision(3) << " kg/s" << std::endl;
+}
+
+template<class P>
+class PrintProblemTest
+{
+public:
+    using Problem = P;
+
+    PrintProblemTest(std::shared_ptr<const Problem> problem)
+    : problem_(problem)
+    {}
+
+    void print()
+    {
+        printProblemTest(*problem_);
+    }
+
+private:
+    std::shared_ptr<const Problem> problem_;
+};
+
+template<class T, class... options>
+void registerPrintProblemTest(pybind11::handle scope, pybind11::class_<T, options...> cls)
+{
+    using Problem = typename T::Problem;
+    cls.def(pybind11::init([](std::shared_ptr<const Problem> problem){
+        return std::make_unique<T>(problem);
+    }));
+    cls.def("print", [](T& self){
+        pybind11::scoped_ostream_redirect stream(std::cout, pybind11::module::import("sys").attr("stdout"));
+        self.print();
+    });
+}
+
+} // end namespace Dumux::Python
+
+#endif
diff --git a/test/python/test_py_fvproblem.py b/test/python/test_py_fvproblem.py
new file mode 100644
index 0000000000000000000000000000000000000000..74ecbf10570cf3fd7fae73c9178b5abc47e5d448
--- /dev/null
+++ b/test/python/test_py_fvproblem.py
@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+
+#  Compile some C++ code for reference
+from dumux.common import *
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+
+# debugging/testing code
+def PrintProblemTest(problem):
+    includes = problem._includes + ["test/python/test_boundaryloop.hh"]
+    typeName = "Dumux::Python::PrintProblemTest<{}>".format(problem._typeName)
+    moduleName = moduleName = "printbs_" + hashIt(problem._typeName)
+    generator = SimpleGenerator("PrintProblemTest", "Dumux::Python")
+    module = generator.load(includes, typeName, moduleName)
+    return module.PrintProblemTest(problem)
+
+############################################################
+# The actual Python test
+############################################################
+from dune.grid import structuredGrid
+from dumux.discretization import GridGeometry
+from dumux.common import BoundaryTypes, FVProblem
+
+gridView = structuredGrid([0,0,0],[1,1,1],[3,3,3])
+
+gridGeometry = GridGeometry(gridView, discMethod="box")
+gridGeometry.update()
+
+@FVProblem(gridGeometry)
+class Problem:
+    numEq = 2
+    name = "python_problem"
+
+    def boundaryTypes(self, element, scv):
+        bTypes = BoundaryTypes(self.numEq)
+        bTypes.setDirichlet()
+        return bTypes
+
+    def dirichlet(self, element, scv):
+        if scv.center()[0] > 0.5:
+            return [0.5, 0.5]
+        else:
+            return [1.0, 0.0]
+
+    def sourceAtPos(self, globalPos):
+        return [globalPos[0]]
+
+problem = Problem()
+print("Name of the problem: {}".format(problem.name))
+print("-- Number of equations: {}".format(problem.numEq))
+print()
+
+# test C++/Python compatibility of the hybrid class
+
+# print some info from C++
+printer = PrintProblemTest(problem)
+printer.print()
+
+# print the same info from Python
+numNeumann = 0
+numDirichlet = 0
+totalSource = 0
+for e in gridView.elements:
+    fvGeometry = problem.gridGeometry().localView()
+    fvGeometry.bind(e)
+    for scv in fvGeometry.scvs():
+        bTypes = problem.boundaryTypes(element=e, scv=scv)
+        if bTypes.isDirichlet():
+            numDirichlet += 1
+        elif bTypes.isNeumann():
+            numNeumann += 1
+        totalSource += problem.sourceAtPos(scv.center())[0]*scv.volume()
+
+print("[python] Found {} Neumann faces and {} Dirichlet faces".format(numNeumann, numDirichlet))
+print("[python] Total source {:.2f} kg/s".format(totalSource))
diff --git a/test/python/test_py_gridgeometry.py b/test/python/test_py_gridgeometry.py
new file mode 100755
index 0000000000000000000000000000000000000000..ca25ba36a1b50bd4d463f7f44be0eb15cfaae392
--- /dev/null
+++ b/test/python/test_py_gridgeometry.py
@@ -0,0 +1,21 @@
+#!/usr/bin/env python3
+
+from dune.grid import structuredGrid
+from dumux.discretization import GridGeometry
+
+gridView = structuredGrid([0,0],[1,1],[5,5])
+
+gridGeometry = GridGeometry(gridView, discMethod="cctpfa")
+gridGeometry.update()
+
+print("The total number of scvs is {}".format(gridGeometry.numScv()))
+print("The total number of scvfs is {}".format(gridGeometry.numScvf()))
+
+for e in gridView.elements:
+    fvGeometry = gridGeometry.localView()
+    fvGeometry.bind(e)
+
+    for scv in fvGeometry.scvs():
+        print("scv dofIndex: {}".format(scv.dofIndex()))
+        print("scv center: {}".format(scv.center()))
+        print("scv volume: {}".format(scv.volume()))