diff --git a/.gitignore b/.gitignore
index 52f825dc43e7d51c909701876fb721c6e71044f6..7a0886a7116cdc461fafb2654c5b9fe9b9d2d26f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,7 @@
 # build system clutter
 build-*
 Testing
+!bin/testing
 
 # auto-saved files
 *~
diff --git a/.gitlab-ci/default.yml b/.gitlab-ci/default.yml
index 012847136f8409bed2566eb1e896a7a229fa814c..83b9bbbe6c90d67ba61f38463579b7bbb12ee172 100644
--- a/.gitlab-ci/default.yml
+++ b/.gitlab-ci/default.yml
@@ -4,6 +4,7 @@ default:
 
 stages:
   - configure
+  - linting
   - select
   - build
   - test
@@ -31,6 +32,29 @@ configure:
     expire_in: 3 hours
 
 
+black (python):
+  stage: linting
+  image: registry.gitlab.com/pipeline-components/black:latest
+  script:
+  # only check the python folder for now (Python code related to bindings)
+  # TODO: maybe extend this to the utility scripts?
+    - black --check --verbose -- python
+    - black --check --verbose -- test/python
+
+
+pylint (python):
+  stage: linting
+  script:
+    - source bin/testing/ci-setup-python-env.sh
+    - |
+      if [ -d build-cmake/python/dumux ] ; then
+        pylint --rcfile=.pylintrc build-cmake/python/dumux
+      fi
+  needs:
+    - job: configure
+      artifacts: true
+
+
 select tests:
   stage: select
   script:
diff --git a/.pylintrc b/.pylintrc
new file mode 100644
index 0000000000000000000000000000000000000000..0706a7b78be43b211cca8ba591f8b6604b88fd3f
--- /dev/null
+++ b/.pylintrc
@@ -0,0 +1,626 @@
+[MASTER]
+
+# A comma-separated list of package or module names from where C extensions may
+# be loaded. Extensions are loading into the active Python interpreter and may
+# run arbitrary code.
+extension-pkg-allow-list=
+
+# A comma-separated list of package or module names from where C extensions may
+# be loaded. Extensions are loading into the active Python interpreter and may
+# run arbitrary code. (This is an alternative name to extension-pkg-allow-list
+# for backward compatibility.)
+extension-pkg-whitelist=
+
+# Return non-zero exit code if any of these messages/categories are detected,
+# even if score is above --fail-under value. Syntax same as enable. Messages
+# specified are enabled, while categories only check already-enabled messages.
+fail-on=
+
+# Specify a score threshold to be exceeded before program exits with error.
+fail-under=10.0
+
+# Files or directories to be skipped. They should be base names, not paths.
+ignore=CVS
+
+# Add files or directories matching the regex patterns to the ignore-list. The
+# regex matches against paths.
+ignore-paths=
+
+# Files or directories matching the regex patterns are skipped. The regex
+# matches against base names, not paths.
+ignore-patterns=
+
+# Python code to execute, usually for sys.path manipulation such as
+# pygtk.require().
+#init-hook=
+
+# Use multiple processes to speed up Pylint. Specifying 0 will auto-detect the
+# number of processors available to use.
+jobs=1
+
+# Control the amount of potential inferred values when inferring a single
+# object. This can help the performance when dealing with large functions or
+# complex, nested conditions.
+limit-inference-results=100
+
+# List of plugins (as comma separated values of python module names) to load,
+# usually to register additional checkers.
+load-plugins=
+
+# Pickle collected data for later comparisons.
+persistent=yes
+
+# When enabled, pylint would attempt to guess common misconfiguration and emit
+# user-friendly hints instead of false-positive error messages.
+suggestion-mode=yes
+
+# Allow loading of arbitrary C extensions. Extensions are imported into the
+# active Python interpreter and may run arbitrary code.
+unsafe-load-any-extension=no
+
+
+[MESSAGES CONTROL]
+
+# Only show warnings with the listed confidence levels. Leave empty to show
+# all. Valid levels: HIGH, INFERENCE, INFERENCE_FAILURE, UNDEFINED.
+confidence=
+
+# Disable the message, report, category or checker with the given id(s). You
+# can either give multiple identifiers separated by comma (,) or put this
+# option multiple times (only on the command line, not in the configuration
+# file where it should appear only once). You can also use "--disable=all" to
+# disable everything first and then reenable specific checks. For example, if
+# you want to run only the similarities checker, you can use "--disable=all
+# --enable=similarities". If you want to run only the classes checker, but have
+# no Warning level messages displayed, use "--disable=all --enable=classes
+# --disable=W".
+disable=print-statement,
+        parameter-unpacking,
+        unpacking-in-except,
+        old-raise-syntax,
+        backtick,
+        long-suffix,
+        old-ne-operator,
+        old-octal-literal,
+        import-star-module-level,
+        non-ascii-bytes-literal,
+        raw-checker-failed,
+        bad-inline-option,
+        locally-disabled,
+        file-ignored,
+        suppressed-message,
+        useless-suppression,
+        deprecated-pragma,
+        use-symbolic-message-instead,
+        apply-builtin,
+        basestring-builtin,
+        buffer-builtin,
+        cmp-builtin,
+        coerce-builtin,
+        execfile-builtin,
+        file-builtin,
+        long-builtin,
+        raw_input-builtin,
+        reduce-builtin,
+        standarderror-builtin,
+        unicode-builtin,
+        xrange-builtin,
+        coerce-method,
+        delslice-method,
+        getslice-method,
+        setslice-method,
+        no-absolute-import,
+        old-division,
+        dict-iter-method,
+        dict-view-method,
+        next-method-called,
+        metaclass-assignment,
+        indexing-exception,
+        raising-string,
+        reload-builtin,
+        oct-method,
+        hex-method,
+        nonzero-method,
+        cmp-method,
+        input-builtin,
+        round-builtin,
+        intern-builtin,
+        unichr-builtin,
+        map-builtin-not-iterating,
+        zip-builtin-not-iterating,
+        range-builtin-not-iterating,
+        filter-builtin-not-iterating,
+        using-cmp-argument,
+        eq-without-hash,
+        div-method,
+        idiv-method,
+        rdiv-method,
+        exception-message-attribute,
+        invalid-str-codec,
+        sys-max-int,
+        bad-python3-import,
+        deprecated-string-function,
+        deprecated-str-translate-call,
+        deprecated-itertools-function,
+        deprecated-types-field,
+        next-method-defined,
+        dict-items-not-iterating,
+        dict-keys-not-iterating,
+        dict-values-not-iterating,
+        deprecated-operator-function,
+        deprecated-urllib-function,
+        xreadlines-attribute,
+        deprecated-sys-function,
+        exception-escape,
+        comprehension-escape,
+        too-few-public-methods
+
+# Enable the message, report, category or checker with the given id(s). You can
+# either give multiple identifier separated by comma (,) or put this option
+# multiple time (only on the command line, not in the configuration file where
+# it should appear only once). See also the "--disable" option for examples.
+enable=c-extension-no-member
+
+
+[REPORTS]
+
+# Python expression which should return a score less than or equal to 10. You
+# have access to the variables 'error', 'warning', 'refactor', and 'convention'
+# which contain the number of messages in each category, as well as 'statement'
+# which is the total number of statements analyzed. This score is used by the
+# global evaluation report (RP0004).
+evaluation=10.0 - ((float(5 * error + warning + refactor + convention) / statement) * 10)
+
+# Template used to display messages. This is a python new-style format string
+# used to format the message information. See doc for all details.
+#msg-template=
+
+# Set the output format. Available formats are text, parseable, colorized, json
+# and msvs (visual studio). You can also give a reporter class, e.g.
+# mypackage.mymodule.MyReporterClass.
+output-format=text
+
+# Tells whether to display a full report or only the messages.
+reports=no
+
+# Activate the evaluation score.
+score=yes
+
+
+[REFACTORING]
+
+# Maximum number of nested blocks for function / method body
+max-nested-blocks=5
+
+# Complete name of functions that never returns. When checking for
+# inconsistent-return-statements if a never returning function is called then
+# it will be considered as an explicit return statement and no message will be
+# printed.
+never-returning-functions=sys.exit,argparse.parse_error
+
+
+[LOGGING]
+
+# The type of string formatting that logging methods do. `old` means using %
+# formatting, `new` is for `{}` formatting.
+logging-format-style=old
+
+# Logging modules to check that the string format arguments are in logging
+# function parameter format.
+logging-modules=logging
+
+
+[SPELLING]
+
+# Limits count of emitted suggestions for spelling mistakes.
+max-spelling-suggestions=4
+
+# Spelling dictionary name. Available dictionaries: none. To make it work,
+# install the 'python-enchant' package.
+spelling-dict=
+
+# List of comma separated words that should be considered directives if they
+# appear and the beginning of a comment and should not be checked.
+spelling-ignore-comment-directives=fmt: on,fmt: off,noqa:,noqa,nosec,isort:skip,mypy:
+
+# List of comma separated words that should not be checked.
+spelling-ignore-words=
+
+# A path to a file that contains the private dictionary; one word per line.
+spelling-private-dict-file=
+
+# Tells whether to store unknown words to the private dictionary (see the
+# --spelling-private-dict-file option) instead of raising a message.
+spelling-store-unknown-words=no
+
+
+[MISCELLANEOUS]
+
+# List of note tags to take in consideration, separated by a comma.
+notes=FIXME,
+      XXX,
+      TODO
+
+# Regular expression of note tags to take in consideration.
+#notes-rgx=
+
+
+[TYPECHECK]
+
+# List of decorators that produce context managers, such as
+# contextlib.contextmanager. Add to this list to register other decorators that
+# produce valid context managers.
+contextmanager-decorators=contextlib.contextmanager
+
+# List of members which are set dynamically and missed by pylint inference
+# system, and so shouldn't trigger E1101 when accessed. Python regular
+# expressions are accepted.
+generated-members=
+
+# Tells whether missing members accessed in mixin class should be ignored. A
+# mixin class is detected if its name ends with "mixin" (case insensitive).
+ignore-mixin-members=yes
+
+# Tells whether to warn about missing members when the owner of the attribute
+# is inferred to be None.
+ignore-none=yes
+
+# This flag controls whether pylint should warn about no-member and similar
+# checks whenever an opaque object is returned when inferring. The inference
+# can return multiple potential results while evaluating a Python object, but
+# some branches might not be evaluated, which results in partial inference. In
+# that case, it might be useful to still emit no-member and other checks for
+# the rest of the inferred objects.
+ignore-on-opaque-inference=yes
+
+# List of class names for which member attributes should not be checked (useful
+# for classes with dynamically set attributes). This supports the use of
+# qualified names.
+ignored-classes=optparse.Values,thread._local,_thread._local
+
+# List of module names for which member attributes should not be checked
+# (useful for modules/projects where namespaces are manipulated during runtime
+# and thus existing member attributes cannot be deduced by static analysis). It
+# supports qualified module names, as well as Unix pattern matching.
+ignored-modules=
+
+# Show a hint with possible names when a member name was not found. The aspect
+# of finding the hint is based on edit distance.
+missing-member-hint=yes
+
+# The minimum edit distance a name should have in order to be considered a
+# similar match for a missing member name.
+missing-member-hint-distance=1
+
+# The total number of similar names that should be taken in consideration when
+# showing a hint for a missing member.
+missing-member-max-choices=1
+
+# List of decorators that change the signature of a decorated function.
+signature-mutators=
+
+
+[VARIABLES]
+
+# List of additional names supposed to be defined in builtins. Remember that
+# you should avoid defining new builtins when possible.
+additional-builtins=
+
+# Tells whether unused global variables should be treated as a violation.
+allow-global-unused-variables=yes
+
+# List of names allowed to shadow builtins
+allowed-redefined-builtins=
+
+# List of strings which can identify a callback function by name. A callback
+# name must start or end with one of those strings.
+callbacks=cb_,
+          _cb
+
+# A regular expression matching the name of dummy variables (i.e. expected to
+# not be used).
+dummy-variables-rgx=_+$|(_[a-zA-Z0-9_]*[a-zA-Z0-9]+?$)|dummy|^ignored_|^unused_
+
+# Argument names that match this expression will be ignored. Default to name
+# with leading underscore.
+ignored-argument-names=_.*|^ignored_|^unused_
+
+# Tells whether we should check for unused import in __init__ files.
+init-import=no
+
+# List of qualified module names which can have objects that can redefine
+# builtins.
+redefining-builtins-modules=six.moves,past.builtins,future.builtins,builtins,io
+
+
+[FORMAT]
+
+# Expected format of line ending, e.g. empty (any line ending), LF or CRLF.
+expected-line-ending-format=
+
+# Regexp for a line that is allowed to be longer than the limit.
+ignore-long-lines=^\s*(# )?<?https?://\S+>?$
+
+# Number of spaces of indent required inside a hanging or continued line.
+indent-after-paren=4
+
+# String used as indentation unit. This is usually "    " (4 spaces) or "\t" (1
+# tab).
+indent-string='    '
+
+# Maximum number of characters on a single line.
+max-line-length=100
+
+# Maximum number of lines in a module.
+max-module-lines=1000
+
+# Allow the body of a class to be on the same line as the declaration if body
+# contains single statement.
+single-line-class-stmt=no
+
+# Allow the body of an if to be on the same line as the test if there is no
+# else.
+single-line-if-stmt=no
+
+
+[SIMILARITIES]
+
+# Ignore comments when computing similarities.
+ignore-comments=yes
+
+# Ignore docstrings when computing similarities.
+ignore-docstrings=yes
+
+# Ignore imports when computing similarities.
+ignore-imports=no
+
+# Ignore function signatures when computing similarities.
+ignore-signatures=no
+
+# Minimum lines number of a similarity.
+min-similarity-lines=10
+
+
+[BASIC]
+
+# Naming style matching correct argument names.
+argument-naming-style=camelCase
+
+# Regular expression matching correct argument names. Overrides argument-
+# naming-style.
+#argument-rgx=
+
+# Naming style matching correct attribute names.
+attr-naming-style=camelCase
+
+# Regular expression matching correct attribute names. Overrides attr-naming-
+# style.
+#attr-rgx=
+
+# Bad variable names which should always be refused, separated by a comma.
+bad-names=foo,
+          bar,
+          baz,
+          toto,
+          tutu,
+          tata
+
+# Bad variable names regexes, separated by a comma. If names match any regex,
+# they will always be refused
+bad-names-rgxs=
+
+# Naming style matching correct class attribute names.
+class-attribute-naming-style=any
+
+# Regular expression matching correct class attribute names. Overrides class-
+# attribute-naming-style.
+#class-attribute-rgx=
+
+# Naming style matching correct class constant names.
+class-const-naming-style=UPPER_CASE
+
+# Regular expression matching correct class constant names. Overrides class-
+# const-naming-style.
+#class-const-rgx=
+
+# Naming style matching correct class names.
+class-naming-style=PascalCase
+
+# Regular expression matching correct class names. Overrides class-naming-
+# style.
+#class-rgx=
+
+# Naming style matching correct constant names.
+const-naming-style=UPPER_CASE
+
+# Regular expression matching correct constant names. Overrides const-naming-
+# style.
+#const-rgx=
+
+# Minimum line length for functions/classes that require docstrings, shorter
+# ones are exempt.
+docstring-min-length=-1
+
+# Naming style matching correct function names.
+function-naming-style=camelCase
+
+# Regular expression matching correct function names. Overrides function-
+# naming-style.
+#function-rgx=
+
+# Good variable names which should always be accepted, separated by a comma.
+good-names=i,
+           j,
+           k,
+           ex,
+           Run,
+           _
+
+# Good variable names regexes, separated by a comma. If names match any regex,
+# they will always be accepted
+good-names-rgxs=
+
+# Include a hint for the correct naming format with invalid-name.
+include-naming-hint=no
+
+# Naming style matching correct inline iteration names.
+inlinevar-naming-style=any
+
+# Regular expression matching correct inline iteration names. Overrides
+# inlinevar-naming-style.
+#inlinevar-rgx=
+
+# Naming style matching correct method names.
+method-naming-style=camelCase
+
+# Regular expression matching correct method names. Overrides method-naming-
+# style.
+#method-rgx=
+
+# Naming style matching correct module names.
+module-naming-style=snake_case
+
+# Regular expression matching correct module names. Overrides module-naming-
+# style.
+#module-rgx=
+
+# Colon-delimited sets of names that determine each other's naming style when
+# the name regexes allow several styles.
+name-group=
+
+# Regular expression which should only match function or class names that do
+# not require a docstring.
+no-docstring-rgx=^_
+
+# List of decorators that produce properties, such as abc.abstractproperty. Add
+# to this list to register other decorators that produce valid properties.
+# These decorators are taken in consideration only for invalid-name.
+property-classes=abc.abstractproperty
+
+# Naming style matching correct variable names.
+variable-naming-style=camelCase
+
+# Regular expression matching correct variable names. Overrides variable-
+# naming-style.
+#variable-rgx=
+
+
+[STRING]
+
+# This flag controls whether inconsistent-quotes generates a warning when the
+# character used as a quote delimiter is used inconsistently within a module.
+check-quote-consistency=no
+
+# This flag controls whether the implicit-str-concat should generate a warning
+# on implicit string concatenation in sequences defined over several lines.
+check-str-concat-over-line-jumps=no
+
+
+[IMPORTS]
+
+# List of modules that can be imported at any level, not just the top level
+# one.
+allow-any-import-level=
+
+# Allow wildcard imports from modules that define __all__.
+allow-wildcard-with-all=no
+
+# Analyse import fallback blocks. This can be used to support both Python 2 and
+# 3 compatible code, which means that the block might have code that exists
+# only in one or another interpreter, leading to false positives when analysed.
+analyse-fallback-blocks=no
+
+# Deprecated modules which should not be used, separated by a comma.
+deprecated-modules=
+
+# Output a graph (.gv or any supported image format) of external dependencies
+# to the given file (report RP0402 must not be disabled).
+ext-import-graph=
+
+# Output a graph (.gv or any supported image format) of all (i.e. internal and
+# external) dependencies to the given file (report RP0402 must not be
+# disabled).
+import-graph=
+
+# Output a graph (.gv or any supported image format) of internal dependencies
+# to the given file (report RP0402 must not be disabled).
+int-import-graph=
+
+# Force import order to recognize a module as part of the standard
+# compatibility libraries.
+known-standard-library=
+
+# Force import order to recognize a module as part of a third party library.
+known-third-party=enchant
+
+# Couples of modules and preferred modules, separated by a comma.
+preferred-modules=
+
+
+[CLASSES]
+
+# Warn about protected attribute access inside special methods
+check-protected-access-in-special-methods=no
+
+# List of method names used to declare (i.e. assign) instance attributes.
+defining-attr-methods=__init__,
+                      __new__,
+                      setUp,
+                      __post_init__
+
+# List of member names, which should be excluded from the protected access
+# warning.
+exclude-protected=_asdict,
+                  _fields,
+                  _replace,
+                  _source,
+                  _make,
+                  _typeName, # Dune-specific
+                  _includes # Dune-specific
+
+# List of valid names for the first argument in a class method.
+valid-classmethod-first-arg=cls
+
+# List of valid names for the first argument in a metaclass class method.
+valid-metaclass-classmethod-first-arg=cls
+
+
+[DESIGN]
+
+# Maximum number of arguments for function / method.
+max-args=5
+
+# Maximum number of attributes for a class (see R0902).
+max-attributes=7
+
+# Maximum number of boolean expressions in an if statement (see R0916).
+max-bool-expr=5
+
+# Maximum number of branch for function / method body.
+max-branches=12
+
+# Maximum number of locals for function / method body.
+max-locals=15
+
+# Maximum number of parents for a class (see R0901).
+max-parents=7
+
+# Maximum number of public methods for a class (see R0904).
+max-public-methods=20
+
+# Maximum number of return / yield for function / method body.
+max-returns=6
+
+# Maximum number of statements in function / method body.
+max-statements=50
+
+# Minimum number of public methods for a class (see R0903).
+min-public-methods=2
+
+
+[EXCEPTIONS]
+
+# Exceptions that will emit a warning when being caught. Defaults to
+# "BaseException, Exception".
+overgeneral-exceptions=BaseException,
+                       Exception
diff --git a/cmake.opts b/cmake.opts
index 540cfb464f9cc8cd6f4d21cae36c0765599d4bc8..73220b6ec05dc062697a12be7dc0a5ec44d6aff9 100644
--- a/cmake.opts
+++ b/cmake.opts
@@ -47,7 +47,7 @@ DUMUX_DISABLE_PROP_MACROS=""
 
 # to enable Python bindings
 DUMUX_ENABLE_PYTHON_BINDINGS=""
-#DUMUX_ENABLE_PYTHON_BINDINGS="-DDUNE_ENABLE_PYTHONBINDINGS=1 -DCMAKE_POSITION_INDEPENDENT_CODE=1"
+#DUMUX_ENABLE_PYTHON_BINDINGS="-DDUNE_ENABLE_PYTHONBINDINGS=1 -DBUILD_SHARED_LIBS=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
diff --git a/dumux/python/CMakeLists.txt b/dumux/python/CMakeLists.txt
index 4f95771a6614eff8b2224375adc22822a3c0c1ec..4065a887e9b43003d786799885741e4598d7e22a 100644
--- a/dumux/python/CMakeLists.txt
+++ b/dumux/python/CMakeLists.txt
@@ -1,2 +1,6 @@
+add_subdirectory(assembly)
 add_subdirectory(common)
 add_subdirectory(discretization)
+add_subdirectory(io)
+add_subdirectory(material)
+add_subdirectory(porousmediumflow)
diff --git a/dumux/python/assembly/CMakeLists.txt b/dumux/python/assembly/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9e57d351c17f6c341f52378466108d62419fb2ad
--- /dev/null
+++ b/dumux/python/assembly/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_PYTHON_ASSEMBLY_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_PYTHON_ASSEMBLY_HEADERS}
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/assembly)
diff --git a/dumux/python/assembly/fvassembler.hh b/dumux/python/assembly/fvassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1ede28dc7574a7a7f20bfa754c90a5489fa86afe
--- /dev/null
+++ b/dumux/python/assembly/fvassembler.hh
@@ -0,0 +1,77 @@
+// -*- 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_FVASSEMBLER_HH
+#define DUMUX_PYTHON_COMMON_FVASSEMBLER_HH
+
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/pybind11/stl.h>
+
+namespace Dumux::Python {
+
+// Python wrapper for the FVAssembler C++ class
+template<class FVAssembler, class... options>
+void registerFVAssembler(pybind11::handle scope, pybind11::class_<FVAssembler, options...> cls)
+{
+    using pybind11::operator""_a;
+
+    using Problem = typename FVAssembler::Problem;
+    using GridGeometry = typename FVAssembler::GridGeometry;
+    using GridVariables = typename FVAssembler::GridVariables;
+    using SolutionVector = typename FVAssembler::ResidualType;
+
+    static_assert(std::is_same_v<GridGeometry, typename Problem::GridGeometry>);
+    cls.def(pybind11::init([](std::shared_ptr<const Problem> problem,
+                              std::shared_ptr<const GridGeometry> gridGeometry,
+                              std::shared_ptr<GridVariables> gridVariables){
+        return std::make_shared<FVAssembler>(problem, gridGeometry, gridVariables);
+    }));
+
+    // TODO assembler with time loop
+
+    cls.def_property_readonly("numDofs", &FVAssembler::numDofs);
+    cls.def_property_readonly("problem", &FVAssembler::problem);
+    cls.def_property_readonly("gridGeometry", &FVAssembler::gridGeometry);
+    cls.def_property_readonly("gridView", &FVAssembler::gridView);
+    cls.def_property_readonly("jacobian", &FVAssembler::jacobian);
+    cls.def_property_readonly("residual", &FVAssembler::residual);
+    cls.def_property_readonly("prevSol", &FVAssembler::prevSol);
+    cls.def_property_readonly("isStationaryProblem", &FVAssembler::isStationaryProblem);
+    cls.def_property_readonly("gridVariables", [](FVAssembler& self) { return self.gridVariables(); });
+
+    cls.def("updateGridVariables", [](FVAssembler& self, const SolutionVector& curSol){
+        self.updateGridVariables(curSol);
+    });
+
+    cls.def("assembleResidual", [](FVAssembler& self, const SolutionVector& curSol){
+        self.assembleResidual(curSol);
+    });
+
+    cls.def("assembleJacobianAndResidual", [](FVAssembler& self, const SolutionVector& curSol){
+        self.assembleJacobianAndResidual(curSol);
+    });
+}
+
+} // end namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/common/boundarytypes.hh b/dumux/python/common/boundarytypes.hh
index 460ca8a155352fb7e678bcef254b92fc7f91a1b2..6a423610a3f5e44a2ffe5bee1f304751d28ad5dc 100644
--- a/dumux/python/common/boundarytypes.hh
+++ b/dumux/python/common/boundarytypes.hh
@@ -41,8 +41,8 @@ void registerBoundaryTypes(pybind11::handle scope, pybind11::class_<BoundaryType
     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);
+    cls.def_property_readonly("isDirichlet", &BoundaryTypes::hasDirichlet);
+    cls.def_property_readonly("isNeumann", &BoundaryTypes::hasNeumann);
 }
 
 template <class BoundaryTypes>
diff --git a/dumux/python/common/fvproblem.hh b/dumux/python/common/fvproblem.hh
index 2f9a1ec8bcdcc8b32292f9bdfe7f9c54672c26b6..21e6904314289d94844d44cfc21ef630f9b54422 100644
--- a/dumux/python/common/fvproblem.hh
+++ b/dumux/python/common/fvproblem.hh
@@ -42,7 +42,7 @@ namespace Dumux::Python {
  * \ingroup Common
  * \brief A C++ wrapper for a Python problem
  */
-template<class GridGeometry_, class PrimaryVariables>
+template<class GridGeometry_, class PrimaryVariables, bool enableInternalDirichletConstraints_>
 class FVProblem
 {
 public:
@@ -59,22 +59,38 @@ public:
     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
+    FVProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+              pybind11::object pyProblem)
+    : gridGeometry_(gridGeometry)
+    , pyProblem_(pyProblem)
+    , name_("python_problem")
+    , paramGroup_("")
     {
-        return pyProblem_.attr("name").template cast<std::string>();
+        if (pybind11::hasattr(pyProblem_, "name"))
+            name_ = pyProblem.attr("name")().template cast<std::string>();
+
+        if (pybind11::hasattr(pyProblem_, "paramGroup"))
+            paramGroup_ = pyProblem.attr("paramGroup")().template cast<std::string>();
     }
 
+    const std::string& name() const
+    { return name_; }
+
+    const std::string& paramGroup() const
+    { return paramGroup_; }
+
     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>();
+        {
+            if (pybind11::hasattr(pyProblem_, "boundaryTypes"))
+                return pyProblem_.attr("boundaryTypes")(element, scv).template cast<BoundaryTypes>();
+            else
+                return pyProblem_.attr("boundaryTypesAtPos")(scv.dofPosition()).template cast<BoundaryTypes>();
+        }
     }
 
     BoundaryTypes boundaryTypes(const Element &element,
@@ -83,7 +99,12 @@ public:
         if constexpr (isBox)
             DUNE_THROW(Dune::InvalidStateException, "boundaryTypes(..., scvf) called for box method.");
         else
-            return pyProblem_.attr("boundaryTypes")(element, scvf).template cast<BoundaryTypes>();
+        {
+            if (pybind11::hasattr(pyProblem_, "boundaryTypes"))
+                return pyProblem_.attr("boundaryTypes")(element, scvf).template cast<BoundaryTypes>();
+            else
+                return pyProblem_.attr("boundaryTypesAtPos")(scvf.ipGlobal()).template cast<BoundaryTypes>();
+        }
     }
 
     PrimaryVariables dirichlet(const Element &element,
@@ -92,7 +113,12 @@ public:
         if constexpr (!isBox)
             DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for cell-centered method.");
         else
-            return pyProblem_.attr("dirichlet")(element, scv).template cast<PrimaryVariables>();
+        {
+            if (pybind11::hasattr(pyProblem_, "dirichlet"))
+                return pyProblem_.attr("dirichlet")(element, scv).template cast<PrimaryVariables>();
+            else
+                return pyProblem_.attr("dirichletAtPos")(scv.dofPosition()).template cast<PrimaryVariables>();
+        }
     }
 
     PrimaryVariables dirichlet(const Element &element,
@@ -101,7 +127,12 @@ public:
         if constexpr (isBox)
             DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for box method.");
         else
-            return pyProblem_.attr("dirichlet")(element, scvf).template cast<PrimaryVariables>();
+        {
+            if (pybind11::hasattr(pyProblem_, "dirichlet"))
+                return pyProblem_.attr("dirichlet")(element, scvf).template cast<PrimaryVariables>();
+            else
+                return pyProblem_.attr("dirichletAtPos")(scvf.ipGlobal()).template cast<PrimaryVariables>();
+        }
     }
 
     template<class ElementVolumeVariables, class ElementFluxVariablesCache>
@@ -111,7 +142,10 @@ public:
                         const ElementFluxVariablesCache& elemFluxVarsCache,
                         const SubControlVolumeFace& scvf) const
     {
-        return pyProblem_.attr("neumann")(element, fvGeometry, scvf).template cast<NumEqVector>();
+        if (pybind11::hasattr(pyProblem_, "neumann"))
+            return pyProblem_.attr("neumann")(element, fvGeometry, scvf).template cast<NumEqVector>();
+        else
+            return pyProblem_.attr("neumannAtPos")(scvf.ipGlobal()).template cast<NumEqVector>();
     }
 
     template<class ElementVolumeVariables>
@@ -120,7 +154,10 @@ public:
                        const ElementVolumeVariables& elemVolVars,
                        const SubControlVolume &scv) const
     {
-        return pyProblem_.attr("source")(element, fvGeometry, scv).template cast<NumEqVector>();
+        if (pybind11::hasattr(pyProblem_, "source"))
+            return pyProblem_.attr("source")(element, fvGeometry, scv).template cast<NumEqVector>();
+        else
+            return pyProblem_.attr("sourceAtPos")(scv.dofPosition()).template cast<NumEqVector>();
     }
 
     NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
@@ -128,6 +165,18 @@ public:
         return pyProblem_.attr("sourceAtPos")(globalPos).template cast<NumEqVector>();
     }
 
+    template<class ElementVolumeVariables>
+    NumEqVector scvPointSources(const Element& element,
+                                const FVElementGeometry& fvGeometry,
+                                const ElementVolumeVariables& elemVolVars,
+                                const SubControlVolume& scv) const
+    {
+        if (pybind11::hasattr(pyProblem_, "scvPointSources"))
+            return pyProblem_.attr("scvPointSources")(element, fvGeometry, scv).template cast<NumEqVector>();
+        else
+            return NumEqVector(0.0);
+    }
+
     template<class Entity>
     PrimaryVariables initial(const Entity& entity) const
     {
@@ -142,12 +191,37 @@ public:
         return pyProblem_.attr("extrusionFactor")(element, scv).template cast<Scalar>();
     }
 
+    Scalar temperatureAtPos(const GlobalPosition& globalPos) const
+    {
+        return pyProblem_.attr("temperatureAtPos")(globalPos).template cast<Scalar>();
+    }
+
+    static constexpr bool enableInternalDirichletConstraints()
+    { return enableInternalDirichletConstraints_; }
+
+    /*!
+     * \brief Add source term derivative to the Jacobian
+     * \note Only needed in case of analytic differentiation and solution dependent sources
+     */
+    template<class MatrixBlock, class VolumeVariables>
+    void addSourceDerivatives(MatrixBlock& block,
+                              const Element& element,
+                              const FVElementGeometry& fvGeometry,
+                              const VolumeVariables& volVars,
+                              const SubControlVolume& scv) const
+    {
+        if (pybind11::hasattr(pyProblem_, "addSourceDerivatives"))
+            pyProblem_.attr("addSourceDerivatives")(block, element, fvGeometry, scv);
+    }
+
     const GridGeometry& gridGeometry() const
     { return *gridGeometry_; }
 
 private:
     std::shared_ptr<const GridGeometry> gridGeometry_;
     pybind11::object pyProblem_;
+    std::string name_;
+    std::string paramGroup_;
 };
 
 // Python wrapper for the above FVProblem C++ class
@@ -158,12 +232,14 @@ void registerFVProblem(pybind11::handle scope, pybind11::class_<Problem, options
     using namespace Dune::Python;
 
     using GridGeometry = typename Problem::GridGeometry;
-    cls.def(pybind11::init([](std::shared_ptr<const GridGeometry> gridGeometry, pybind11::object p){
+    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; });
+    cls.def_property_readonly("gridGeometry", &Problem::gridGeometry);
 
     using GridView = typename GridGeometry::GridView;
     using Element = typename GridView::template Codim<0>::Entity;
@@ -188,7 +264,6 @@ void registerFVProblem(pybind11::handle scope, pybind11::class_<Problem, options
     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
diff --git a/dumux/python/common/timeloop.hh b/dumux/python/common/timeloop.hh
index 9117c0ee48fd13557549fc2eb1385e3a407ebbb0..241150807d029e51f12a13f574d15a7d9d27b4c7 100644
--- a/dumux/python/common/timeloop.hh
+++ b/dumux/python/common/timeloop.hh
@@ -42,18 +42,18 @@ void registerTimeLoop(pybind11::handle scope,
         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_property_readonly("time", &TimeLoop::time);
+    cls.def_property_readonly("timeStepSize", &TimeLoop::timeStepSize);
+    cls.def_property_readonly("finished", &TimeLoop::finished);
+    cls.def_property_readonly("isCheckPoint", &TimeLoop::isCheckPoint);
     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);
     });
diff --git a/dumux/python/common/volumevariables.hh b/dumux/python/common/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cce94ab1cbe7cb55abdb0b34761472961d4d440b
--- /dev/null
+++ b/dumux/python/common/volumevariables.hh
@@ -0,0 +1,82 @@
+#ifndef DUMUX_PYTHON_COMMON_VOLUMEVARIABLES_HH
+#define DUMUX_PYTHON_COMMON_VOLUMEVARIABLES_HH
+
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/pybind11/stl.h>
+
+#include <dune/python/common/typeregistry.hh>
+#include <dune/common/classname.hh>
+
+namespace Dumux::Python::Impl {
+
+// helper structs and functions for detecting member functions in volVars
+template <class VolumeVariables>
+using PhaseTemperatureDetector = decltype(std::declval<VolumeVariables>().temperature(0));
+
+template<class VolumeVariables>
+static constexpr bool hasPhaseTemperature()
+{ return Dune::Std::is_detected<PhaseTemperatureDetector, VolumeVariables>::value; }
+
+template <class VolumeVariables>
+using MoleFractionDetector = decltype(std::declval<VolumeVariables>().moleFraction(0, 0));
+
+template<class VolumeVariables>
+static constexpr bool hasMoleFraction()
+{ return Dune::Std::is_detected<MoleFractionDetector, VolumeVariables>::value; }
+
+template <class VolumeVariables>
+using MassFractionDetector = decltype(std::declval<VolumeVariables>().massFraction(0, 0));
+
+template<class VolumeVariables>
+static constexpr bool hasMassFraction()
+{ return Dune::Std::is_detected<MassFractionDetector, VolumeVariables>::value; }
+
+template <class VolumeVariables>
+using SaturationDetector = decltype(std::declval<VolumeVariables>().saturation(0));
+
+template<class VolumeVariables>
+static constexpr bool hasSaturation()
+{ return Dune::Std::is_detected<SaturationDetector, VolumeVariables>::value; }
+
+template <class VolumeVariables>
+using PermeabilityDetector = decltype(std::declval<VolumeVariables>().permeability());
+
+template<class VolumeVariables>
+static constexpr bool hasPermeability()
+{ return Dune::Std::is_detected<PermeabilityDetector, VolumeVariables>::value; }
+
+template<class VolumeVariables>
+void registerVolumeVariables(pybind11::handle scope)
+{
+    using namespace Dune::Python;
+
+    auto [cls, addedToRegistry] = insertClass<VolumeVariables>(
+        scope, "VolumeVariables",
+        GenerateTypeName(Dune::className<VolumeVariables>()),
+        IncludeFiles{}
+    );
+
+    if (addedToRegistry)
+    {
+        using pybind11::operator""_a;
+
+        cls.def("pressure", &VolumeVariables::pressure, "phaseIdx"_a=0);
+        cls.def("density", &VolumeVariables::density, "phaseIdx"_a=0);
+        cls.def("temperature", &VolumeVariables::temperature);
+
+        if constexpr(hasSaturation<VolumeVariables>())
+            cls.def("saturation", &VolumeVariables::saturation, "saturation"_a=0);
+        if constexpr(hasPhaseTemperature<VolumeVariables>())
+            cls.def("temperature", &VolumeVariables::temperature, "phaseIdx"_a=0);
+        if constexpr(hasMoleFraction<VolumeVariables>())
+            cls.def("moleFraction", &VolumeVariables::moleFraction, "phaseIdx"_a, "compIdx"_a);
+        if constexpr(hasMassFraction<VolumeVariables>())
+            cls.def("massFraction", &VolumeVariables::massFraction, "phaseIdx"_a, "compIdx"_a);
+        if constexpr(hasPermeability<VolumeVariables>())
+            cls.def("permeability", &VolumeVariables::permeability);
+    }
+}
+
+} // namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/discretization/gridgeometry.hh b/dumux/python/discretization/gridgeometry.hh
index 69c6af1f09965ab30f07b0b893e938499bccb6f9..34a4c6a806197d20eac6ea97f60c51067e758e27 100644
--- a/dumux/python/discretization/gridgeometry.hh
+++ b/dumux/python/discretization/gridgeometry.hh
@@ -48,12 +48,12 @@ void registerSubControlVolume(pybind11::handle scope)
     {
         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);
+        cls.def_property_readonly("center", &SCV::center);
+        cls.def_property_readonly("volume", &SCV::volume);
+        cls.def_property_readonly("dofIndex", &SCV::dofIndex);
+        cls.def_property_readonly("localDofIndex", &SCV::localDofIndex);
+        cls.def_property_readonly("dofPosition", &SCV::dofPosition);
+        cls.def_property_readonly("elementIndex", &SCV::elementIndex);
 
     }
 }
@@ -73,14 +73,14 @@ void registerSubControlVolumeFace(pybind11::handle scope)
     {
         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);
+        cls.def_property_readonly("center", &SCVF::center);
+        cls.def_property_readonly("area", &SCVF::area);
+        cls.def_property_readonly("ipGlobal", &SCVF::ipGlobal);
+        cls.def_property_readonly("boundary", &SCVF::boundary);
+        cls.def_property_readonly("unitOuterNormal", &SCVF::unitOuterNormal);
+        cls.def_property_readonly("insideScvIdx", &SCVF::insideScvIdx);
+        cls.def_property_readonly("outsideScvIdx", [](SCVF& self){ return self.outsideScvIdx(); });
+        cls.def_property_readonly("index", &SCVF::index);
     }
 }
 
@@ -99,15 +99,16 @@ void registerFVElementGeometry(pybind11::handle scope)
     {
         using pybind11::operator""_a;
 
-        cls.def("numScvf", &FVEG::numScv);
-        cls.def("numScv", &FVEG::numScv);
+        cls.def_property_readonly("numScvf", &FVEG::numScv);
+        cls.def_property_readonly("numScv", &FVEG::numScv);
+        cls.def_property_readonly("hasBoundaryScvf", &FVEG::hasBoundaryScvf);
         cls.def("bind", &FVEG::bind, "element"_a);
-        cls.def("hasBoundaryScvf", &FVEG::hasBoundaryScvf);
-        cls.def("scvs", [](FVEG& self){
+        cls.def("bindElement", &FVEG::bindElement, "element"_a);
+        cls.def_property_readonly("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){
+        cls.def_property_readonly("scvfs", [](FVEG& self){
             const auto range = scvfs(self);
             return pybind11::make_iterator(range.begin(), range.end());
         }, pybind11::keep_alive<0, 1>());
@@ -124,14 +125,40 @@ void registerGridGeometry(pybind11::handle scope, pybind11::class_<GG, Options..
 
     using GridView = typename GG::GridView;
 
-    cls.def(pybind11::init([](GridView gridView){
-        return std::make_shared<GG>(gridView);
+    cls.def(pybind11::init([](const GridView& gridView){
+        auto gg = std::make_shared<GG>(gridView);
+        // remove this once the interface is changed on the C++ side (see #1056)
+        gg->update();
+        return gg;
     }), "gridView"_a);
 
-    cls.def("update", &GG::update);
-    cls.def("numDofs", &GG::numDofs);
-    cls.def("numScv", &GG::numScv);
-    cls.def("numScvf", &GG::numScvf);
+    // update this once the interface is changed on the C++ side (see #1056)
+    cls.def("update", [](GG& self, const GridView& gridView){
+        return self.update();
+    }, "gridView"_a);
+
+    cls.def_property_readonly("numDofs", &GG::numDofs);
+    cls.def_property_readonly("numScv", &GG::numScv);
+    cls.def_property_readonly("numScvf", &GG::numScvf);
+    cls.def_property_readonly("bBoxMax", &GG::bBoxMax);
+    cls.def_property_readonly("bBoxMin", &GG::bBoxMin);
+    cls.def_property_readonly("gridView", &GG::gridView);
+
+    cls.def_property_readonly_static("discMethod", [](const pybind11::object&){
+        return toString(GG::discMethod);
+    });
+
+    cls.def_property_readonly("localView", [](GG& self){
+        return localView(self);
+    });
+
+    using LocalView = typename GG::LocalView;
+    using Element = typename LocalView::Element;
+    cls.def("boundLocalView", [](GG& self, const Element& element){
+        auto view = localView(self);
+        view.bind(element);
+        return view;
+    }, "element"_a);
 
     using SubControlVolume = typename GG::SubControlVolume;
     Impl::registerSubControlVolume<SubControlVolume>(scope);
@@ -140,10 +167,8 @@ void registerGridGeometry(pybind11::handle scope, pybind11::class_<GG, Options..
     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
diff --git a/dumux/python/discretization/gridvariables.hh b/dumux/python/discretization/gridvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f2edf6039a6f2406da4b8f2f9cc85f5a9a31cfa1
--- /dev/null
+++ b/dumux/python/discretization/gridvariables.hh
@@ -0,0 +1,109 @@
+// -*- 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_GRIDVARIABLES_HH
+#define DUMUX_PYTHON_DISCRETIZATION_GRIDVARIABLES_HH
+
+#include <memory>
+#include <dune/istl/bvector.hh>
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/common/typeregistry.hh>
+#include <dumux/discretization/elementsolution.hh>
+
+namespace Dumux::Python {
+
+namespace Impl {
+
+template<class ElementSolution>
+void registerElementSolution(pybind11::handle scope)
+{
+    using namespace Dune::Python;
+
+    auto [cls, addedToRegistry] = insertClass<ElementSolution>(
+        scope, "ElementSolution",
+        GenerateTypeName(Dune::className<ElementSolution>()),
+        IncludeFiles{"dumux/discretization/elementsolution.hh"}
+    );
+
+    if (addedToRegistry)
+    {
+        using pybind11::operator""_a;
+
+        cls.def("__getitem__", [](const ElementSolution& self, std::size_t i){
+            if (i >= self.size())
+                throw pybind11::index_error();
+            return self[i];
+        });
+
+        cls.def_property_readonly("size", &ElementSolution::size);
+    }
+}
+
+} // end namespace Impl
+
+
+// see python/dumux/discretization/__init__.py for how this is used for JIT compilation
+template <class GV, class... Options>
+void registerGridVariables(pybind11::handle scope, pybind11::class_<GV, Options...> cls)
+{
+    using pybind11::operator""_a;
+
+    using Problem = typename GV::GridVolumeVariables::Problem;
+    using PrimaryVariables = typename GV::GridVolumeVariables::VolumeVariables::PrimaryVariables;
+    using GridGeometry = typename GV::GridGeometry;
+    using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+    using SolutionVector = Dune::BlockVector<PrimaryVariables>;
+
+    cls.def(pybind11::init([](std::shared_ptr<const Problem> problem,
+                              std::shared_ptr<const GridGeometry> gridGeometry){
+        return std::make_shared<GV>(problem, gridGeometry);
+    }));
+
+    cls.def("init", [](GV& self, const SolutionVector& sol) { return self.init(sol); });
+    cls.def("advanceTimeStep", &GV::advanceTimeStep);
+    cls.def_property_readonly("curGridVolVars", [](GV& self) { return self.curGridVolVars(); });
+    cls.def_property_readonly("gridFluxVarsCache", [](GV& self) { return self.gridFluxVarsCache(); });
+    cls.def_property_readonly("prevGridVolVars", [](GV& self) { return self.prevGridVolVars(); });
+    cls.def_property_readonly("gridGeometry", &GV::gridGeometry);
+
+    cls.def("updateAfterGridAdaption", [](GV& self, const SolutionVector& sol){
+        return self.updateAfterGridAdaption(sol);
+    });
+
+    cls.def("resetTimeStep", [](GV& self, const SolutionVector& sol){
+        return self.resetTimeStep(sol);
+    });
+
+    cls.def("update", [](GV& self, const SolutionVector& sol, const bool forceFluxCacheUpdate = false){
+        return self.update(sol, forceFluxCacheUpdate);
+    });
+
+    using ElementSolution = std::decay_t<decltype(elementSolution(std::declval<Element>(),
+                                                                  std::declval<SolutionVector>(),
+                                                                  std::declval<GridGeometry>()))>;
+    Impl::registerElementSolution<ElementSolution>(scope);
+}
+
+} // namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/io/CMakeLists.txt b/dumux/python/io/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ba21bd4a24f5096697bc1295385dee239d5520fa
--- /dev/null
+++ b/dumux/python/io/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_PYTHON_IO_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_PYTHON_IO_HEADERS}
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/io)
diff --git a/dumux/python/io/vtkoutputmodule.hh b/dumux/python/io/vtkoutputmodule.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7272fc4e5ad88ce276190f83d4ab38ae1ab01f20
--- /dev/null
+++ b/dumux/python/io/vtkoutputmodule.hh
@@ -0,0 +1,87 @@
+// -*- 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_IO_VTK_OUTPUTMODULE_HH
+#define DUMUX_PYTHON_IO_VTK_OUTPUTMODULE_HH
+
+#include <dumux/io/vtkoutputmodule.hh>
+
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/pybind11/stl.h>
+
+#include <dumux/python/common/volumevariables.hh>
+#include <dumux/io/velocityoutput.hh>
+
+namespace Dumux::Python {
+
+template <class GridVariables, class SolutionVector, class... options>
+void registerVtkOutputModule(pybind11::handle scope,
+                             pybind11::class_<VtkOutputModule<GridVariables, SolutionVector>, options...> cls)
+{
+    using pybind11::operator""_a;
+
+    using VtkOutputModule = Dumux::VtkOutputModule<GridVariables, SolutionVector>;
+    using VolumeVariables = typename VtkOutputModule::VolumeVariables;
+    Dumux::Python::Impl::registerVolumeVariables<VolumeVariables>(scope);
+
+
+    cls.def(pybind11::init([](const GridVariables& gridVariables,
+                              const SolutionVector& sol,
+                              const std::string& name){
+        return new VtkOutputModule(gridVariables, sol, name);
+    }));
+
+    using Scalar = double;
+
+    cls.def("addField", [](VtkOutputModule& self, const SolutionVector& sol, const std::string& name){
+        self.addField(sol, name);
+    });
+
+    cls.def("write", [](VtkOutputModule& self, Scalar time){
+        self.write(time);
+    });
+
+    cls.def("addVolumeVariable", [](VtkOutputModule& self,
+                                    std::function<Scalar(const VolumeVariables&)>&& f,
+                                    const std::string& name){
+        self.addVolumeVariable(std::move(f), name);
+    });
+
+    using VelocityOutputType = Dumux::VelocityOutput<GridVariables>;
+    cls.def("addVelocityOutput", [](VtkOutputModule& self, std::shared_ptr<VelocityOutputType> velocityOutput){
+        self.addVelocityOutput(velocityOutput);
+    });
+};
+
+
+template<class GridVariables, class SolutionVector>
+void registerVtkOutputModule(pybind11::handle scope, const char *clsName = "VtkOutputModule")
+{
+    using VtkOutputModule = Dumux::VtkOutputModule<GridVariables, SolutionVector>;
+    pybind11::class_<VtkOutputModule> cls(scope, clsName);
+    registerVtkOutputModule(scope, cls);
+}
+
+} // namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/material/CMakeLists.txt b/dumux/python/material/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f1c168afc35b03db804df7829820e885bd607308
--- /dev/null
+++ b/dumux/python/material/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_subdirectory(components)
+add_subdirectory(fluidsystems)
+add_subdirectory(spatialparams)
diff --git a/dumux/python/material/components/CMakeLists.txt b/dumux/python/material/components/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7aaff78f8d4d2739a3c718c0abc092f5ff5daf51
--- /dev/null
+++ b/dumux/python/material/components/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_PYTHON_MATERIAL_COMPONENTS_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_PYTHON_MATERIAL_COMPONENTS_HEADERS}
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/material/components)
diff --git a/dumux/python/material/components/component.hh b/dumux/python/material/components/component.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9ad1b1ca87ce468bbf2cb32ed8b476cc4a7c11a1
--- /dev/null
+++ b/dumux/python/material/components/component.hh
@@ -0,0 +1,75 @@
+#ifndef DUMUX_PYTHON_MATERIAL_COMPONENT_HH
+#define DUMUX_PYTHON_MATERIAL_COMPONENT_HH
+
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/pybind11/stl.h>
+
+#include <dumux/material/components/componenttraits.hh>
+
+
+namespace Dumux::Python {
+
+
+template <class Comp, class... options>
+void registerComponent(pybind11::handle scope,
+                       pybind11::class_<Comp, options...> cls)
+{
+    using pybind11::operator""_a;
+
+    cls.def(pybind11::init());
+
+    cls.def_property_readonly_static("name", &Comp::name);
+    cls.def_property_readonly_static("molarMass", &Comp::molarMass);
+
+    if constexpr (ComponentTraits<Comp>::hasLiquidState)
+    {
+        cls.def_static("liquidDensity", &Comp::liquidDensity, "temperature"_a, "pressure"_a);
+        cls.def_static("liquidMolarDensity", &Comp::liquidMolarDensity, "temperature"_a, "pressure"_a);
+        cls.def_property_readonly_static("liquidIsCompressible", &Comp::liquidIsCompressible);
+        cls.def_static("liquidViscosity", &Comp::liquidViscosity, "temperature"_a, "pressure"_a);
+        cls.def_static("liquidEnthalpy", &Comp::liquidEnthalpy, "temperature"_a, "pressure"_a);
+        cls.def_static("liquidInternalEnergy", &Comp::liquidInternalEnergy, "temperature"_a, "pressure"_a);
+        cls.def_static("liquidHeatCapacity", &Comp::liquidHeatCapacity, "temperature"_a, "pressure"_a);
+        cls.def_static("liquidThermalConductivity", &Comp::liquidThermalConductivity, "temperature"_a, "pressure"_a);
+        cls.def_static("vaporPressure", &Comp::vaporPressure, "temperature"_a);
+    }
+
+    if constexpr (ComponentTraits<Comp>::hasGasState)
+    {
+        using Scalar = typename Comp::Scalar;
+        cls.def_static("gasDensity", &Comp::gasDensity, "temperature"_a, "pressure"_a);
+        cls.def_static("gasMolarDensity", &Comp::gasMolarDensity, "temperature"_a, "pressure"_a);
+        cls.def_property_readonly_static("gasIsCompressible", &Comp::gasIsCompressible);
+        cls.def_property_readonly_static("gasIsIdeal", &Comp::gasIsIdeal);
+        cls.def_static("gasPressure", &Comp::gasPressure, "temperature"_a, "pressure"_a);
+        cls.def_static("gasViscosity", &Comp::gasViscosity, "temperature"_a, "pressure"_a);
+        cls.def_static("gasEnthalpy", &Comp::gasEnthalpy, "temperature"_a, "pressure"_a);
+        cls.def_static("gasInternalEnergy", &Comp::gasInternalEnergy, "temperature"_a, "pressure"_a);
+        cls.def_static("gasHeatCapacity", &Comp::gasHeatCapacity, "temperature"_a, "pressure"_a);
+        cls.def_static("gasThermalConductivity", &Comp::gasThermalConductivity, "temperature"_a, "pressure"_a);
+    }
+
+    if constexpr (ComponentTraits<Comp>::hasSolidState)
+    {
+        cls.def_property_readonly_static("solidIsCompressible", &Comp::solidIsCompressible);
+        cls.def_static("solidDensity", &Comp::solidDensity, "temperature_a");
+        cls.def_static("solidThermalConductivity", &Comp::solidThermalConductivity, "temperature_a");
+        cls.def_static("solidHeatCapacity", &Comp::solidHeatCapacity, "temperature_a");
+    }
+
+    if constexpr (ComponentTraits<Comp>::isIon)
+    {
+        cls.def_property_readonly_static("charge", &Comp::charge);
+    }
+
+    if constexpr (ComponentTraits<Comp>::hasLiquidState || ComponentTraits<Comp>::hasGasState)
+    {
+        cls.def_property_readonly_static("criticalTemperature", &Comp::criticalTemperature);
+        cls.def_property_readonly_static("criticalPressure", &Comp::criticalPressure);
+    }
+}
+
+
+} // namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/material/fluidsystems/CMakeLists.txt b/dumux/python/material/fluidsystems/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..91ab2273ca4e04474252d75a848cf29211f1c414
--- /dev/null
+++ b/dumux/python/material/fluidsystems/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_PYTHON_MATERIAL_FLUIDSYSTEMS_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_PYTHON_MATERIAL_FLUIDSYSTEMS_HEADERS}
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/material/fluidsystems)
diff --git a/dumux/python/material/fluidsystems/fluidsystem.hh b/dumux/python/material/fluidsystems/fluidsystem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7cf246b9151470900b072b9562221f78cc6ffb80
--- /dev/null
+++ b/dumux/python/material/fluidsystems/fluidsystem.hh
@@ -0,0 +1,41 @@
+// -*- 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_MATERIAL_FLUIDSYSTEMS_FLUIDSYSTEM_HH
+#define DUMUX_PYTHON_MATERIAL_FLUIDSYSTEMS_FLUIDSYSTEM_HH
+
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/pybind11/stl.h>
+
+namespace Dumux::Python {
+
+template <class FS, class... options>
+void registerFluidSystem(pybind11::handle scope,
+                         pybind11::class_<FS, options...> cls)
+{
+    cls.def(pybind11::init());
+}
+
+} // namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/material/spatialparams/CMakeLists.txt b/dumux/python/material/spatialparams/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..08816c5b1a287597d371c037e0e4703c0496b7cf
--- /dev/null
+++ b/dumux/python/material/spatialparams/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_PYTHON_MATERIAL_SPATIALPARAMS_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_PYTHON_MATERIAL_SPATIALPARAMS_HEADERS}
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/material/spatialparams)
diff --git a/dumux/python/material/spatialparams/spatialparams.hh b/dumux/python/material/spatialparams/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..448a981491bcfca1747799614d79ee7be53d6b3b
--- /dev/null
+++ b/dumux/python/material/spatialparams/spatialparams.hh
@@ -0,0 +1,98 @@
+// -*- 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_MATERIAL_SPATIAL_PARAMS_HH
+#define DUMUX_PYTHON_MATERIAL_SPATIAL_PARAMS_HH
+
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/pybind11/stl.h>
+
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux::Python {
+
+template<class GridGeometry, class Scalar, class PT>
+class FVSpatialParamsOneP
+: public Dumux::FVSpatialParamsOneP<GridGeometry, Scalar, FVSpatialParamsOneP<GridGeometry, Scalar, PT>>
+{
+    using ThisType = FVSpatialParamsOneP<GridGeometry, Scalar, PT>;
+    using ParentType = Dumux::FVSpatialParamsOneP<GridGeometry, Scalar, ThisType>;
+    using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+    using SubControlVolume = typename GridGeometry::SubControlVolume;
+    using GlobalPosition = typename GridGeometry::GridView::template Codim<0>::Geometry::GlobalCoordinate;
+
+public:
+
+    using PermeabilityType = PT;
+
+    FVSpatialParamsOneP(std::shared_ptr<const GridGeometry> gridGeometry,
+                        pybind11::object pySpatialParams)
+    : ParentType(gridGeometry)
+    , pySpatialParams_(pySpatialParams)
+    {}
+
+    template<class ElementSolution>
+    PermeabilityType permeability(const Element& element,
+                                  const SubControlVolume& scv,
+                                  const ElementSolution& elemSol) const
+    {
+        if (pybind11::hasattr(pySpatialParams_, "permeability"))
+            return pySpatialParams_.attr("permeability")(element, scv, elemSol).template cast<PermeabilityType>();
+        else
+            return pySpatialParams_.attr("permeabilityAtPos")(scv.center()).template cast<PermeabilityType>();
+    }
+
+    template<class ElementSolution>
+    Scalar porosity(const Element& element,
+                              const SubControlVolume& scv,
+                              const ElementSolution& elemSol) const
+    {
+        if (pybind11::hasattr(pySpatialParams_, "porosity"))
+            return pySpatialParams_.attr("porosity")(element, scv, elemSol).template cast<Scalar>();
+        else
+            return pySpatialParams_.attr("porosityAtPos")(scv.center()).template cast<Scalar>();
+    }
+
+private:
+    pybind11::object pySpatialParams_;
+};
+
+
+
+template <class SP, class... options>
+void registerOnePSpatialParams(pybind11::handle scope,
+                               pybind11::class_<SP, options...> cls)
+{
+    using pybind11::operator""_a;
+    using GridGeometry = std::decay_t<decltype(std::declval<SP>().gridGeometry())>;
+
+    cls.def(pybind11::init([](std::shared_ptr<GridGeometry> gridGeometry, pybind11::object sp){
+        return std::make_shared<SP>(gridGeometry, sp);
+    }));
+
+    cls.def("gravity", &SP::gravity);
+}
+
+} // namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/porousmediumflow/CMakeLists.txt b/dumux/python/porousmediumflow/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0b1c7f97be7f856efd1886dcd6b1b499f250fb88
--- /dev/null
+++ b/dumux/python/porousmediumflow/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_PYTHON_POROUSMEDIUMFLOW_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_PYTHON_POROUSMEDIUMFLOW_HEADERS}
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/porousmediumflow)
diff --git a/dumux/python/porousmediumflow/problem.hh b/dumux/python/porousmediumflow/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..30c72f3b405c489559dc6d940f453f5c24cf2061
--- /dev/null
+++ b/dumux/python/porousmediumflow/problem.hh
@@ -0,0 +1,119 @@
+// -*- 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_POROUSMEDIUMFLOW_PROBLEM_HH
+#define DUMUX_PYTHON_POROUSMEDIUMFLOW_PROBLEM_HH
+
+
+#include <dune/python/pybind11/pybind11.h>
+
+#include <dumux/python/common/fvproblem.hh>
+
+namespace Dumux::Python {
+
+/*!
+ * \ingroup Common
+ * \brief A C++ wrapper for a Python PorousMediumFlow problem
+ */
+template<class GridGeometry_, class PrimaryVariables, class SpatialParams_,
+         bool enableInternalDirichletConstraints>
+class PorousMediumFlowProblem : public FVProblem<GridGeometry_, PrimaryVariables, enableInternalDirichletConstraints>
+{
+    using ParentType = FVProblem<GridGeometry_, PrimaryVariables, enableInternalDirichletConstraints>;
+public:
+    using GridGeometry = GridGeometry_;
+    using SpatialParams = SpatialParams_;
+    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>;
+
+    PorousMediumFlowProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+                            std::shared_ptr<const SpatialParams> spatialParams,
+                            pybind11::object pyProblem)
+    : ParentType(gridGeometry, pyProblem)
+    , spatialParams_(spatialParams)
+    {}
+
+    const SpatialParams& spatialParams() const
+    { return *spatialParams_; }
+
+private:
+    std::shared_ptr<const SpatialParams> spatialParams_;
+};
+
+// Python wrapper for the above FVProblem C++ class
+template<class Problem, class... options>
+void registerPorousMediumFlowProblem(pybind11::handle scope, pybind11::class_<Problem, options...> cls)
+{
+    using pybind11::operator""_a;
+    using namespace Dune::Python;
+
+    using GridGeometry = typename Problem::GridGeometry;
+    using SpatialParams = typename Problem::SpatialParams;
+    cls.def(pybind11::init([](std::shared_ptr<const GridGeometry> gridGeometry,
+                              std::shared_ptr<const SpatialParams> spatialParams,
+                              pybind11::object p){
+        return std::make_shared<Problem>(gridGeometry, spatialParams, 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);
+    cls.def("spatialParams", &Problem::spatialParams);
+}
+
+} // end namespace Dumux::Python
+
+#endif
diff --git a/dumux/python/porousmediumflow/velocityoutput.hh b/dumux/python/porousmediumflow/velocityoutput.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1b53c62d7c608bc7526b1561760ffb8b4285eca1
--- /dev/null
+++ b/dumux/python/porousmediumflow/velocityoutput.hh
@@ -0,0 +1,48 @@
+// -*- 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_POROUSMEDIUMFLOW_VELOCITYOUTPUT_HH
+#define DUMUX_PYTHON_POROUSMEDIUMFLOW_VELOCITYOUTPUT_HH
+
+#include <dune/python/pybind11/pybind11.h>
+#include <dumux/porousmediumflow/velocityoutput.hh>
+
+namespace Dumux::Python {
+
+// Python wrapper for the PorousMediumFlowVelocityOutput class
+template<class GridVariables, class FluxVariables, class... options>
+void registerPorousMediumFlowVelocityOutput(pybind11::handle scope,
+                                            pybind11::class_<Dumux::PorousMediumFlowVelocityOutput<GridVariables, FluxVariables>, options...> cls)
+{
+    using pybind11::operator""_a;
+    using namespace Dune::Python;
+    using VelocityOutput = Dumux::PorousMediumFlowVelocityOutput<GridVariables, FluxVariables>;
+
+    cls.def(pybind11::init([](const GridVariables& gridVariables){
+        return std::make_shared<VelocityOutput>(gridVariables);
+    }));
+}
+
+} // end namespace Dumux::Python
+
+#endif
diff --git a/pyproject.toml b/pyproject.toml
index 562fd5e7ec1eafe724c9483195224d8b50674ce7..ec40a5cb831ca393b986d3e4ac950d18ecaf9709 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,3 +1,6 @@
 [build-system]
 requires = ['setuptools', 'wheel', 'scikit-build', 'cmake', 'ninja', 'requests', 'dune-common<=2.8.0.dev20201218', 'dune-grid<=2.8.0.dev20201218', 'dune-localfunctions<=2.8.0.dev20201218', 'dune-istl<=2.8.0.dev20201218']
 build-backend = 'setuptools.build_meta'
+
+[tool.black]
+line-length = 100
diff --git a/python/README.md b/python/README.md
index 851ef4e6fb803e4b9a424e1666b9f782935c8181..5fb8c28584a770b70d46ab78ccff1182930e98ad 100644
--- a/python/README.md
+++ b/python/README.md
@@ -40,11 +40,11 @@ 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:
+Add the Python binding modules to your Python path like this and install them with
+the setup environment setup script
 
 ```
-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
+source dumux/python/setup-python-env.sh
 ```
 
 If you are getting error with loading MPI in Python you might need to preload the MPI library
@@ -69,3 +69,29 @@ You can run all currently existing DuMu<sup>x</sup> Python tests with
 cd dumux/build-cmake
 ctest -L python
 ```
+##  Development
+
+All Python files should be linted by the tool [`black`](https://pypi.org/project/black/).
+You can install `black` with `pip install black` and run it from the dumux top-directory
+
+```
+black ./python
+```
+
+You can also run it on a specific file (replace `./python` by file name)
+This will automatically format the Python files. Run black before every commit changing Python files.
+
+The `dumux` Python module should be get a score of `10` from
+the tool [`pylint`](https://pypi.org/project/pylint/).
+You can install `pylint` with `pip install pylint` and run it from the dumux top-directory
+
+```
+pylint build-cmake/python/dumux
+```
+
+Pylint needs to be able to check imports so the modules need to be properly set up
+with `setup-dunepy.py` (see above). The `pylint` configuration file `dumux/.pylintrc` can
+be used to configure `pylint`. Some exceptions or other parameters than the default
+might be sensible in the future but generally advice given by `pylint` leads to better code.
+Different from `black`, `pylint` does no itself fix the code, you need to do this yourself.
+Always run `black` before checking `pylint`.
diff --git a/python/dumux/CMakeLists.txt b/python/dumux/CMakeLists.txt
index e11f39eed6571d247f90fbad0bc03b5324db7e29..e45fd127636d0e2f73081f762b85657a63c44cb0 100644
--- a/python/dumux/CMakeLists.txt
+++ b/python/dumux/CMakeLists.txt
@@ -1,5 +1,10 @@
+add_subdirectory(assembly)
 add_subdirectory(common)
 add_subdirectory(discretization)
+add_subdirectory(material)
+add_subdirectory(io)
+add_subdirectory(porousmediumflow)
+add_subdirectory(wrapping)
 
 add_python_targets(dumux
   __init__
diff --git a/python/dumux/__init__.py b/python/dumux/__init__.py
index de40ea7ca058e07a399acf61529d418c07eeee61..bb4c6ae25606e22f49e4e73de4426d1721d0e5d3 100644
--- a/python/dumux/__init__.py
+++ b/python/dumux/__init__.py
@@ -1 +1,13 @@
-__import__('pkg_resources').declare_namespace(__name__)
+"""The DuMux Python-bindings for using DuMux entirely from Python
+
+DuMux is
+* short for Dune for Multi-{Phase, Component, Scale, Physics, …} flow and transport in porous media
+* a free and open-source simulator for flow and transport processes in porous media
+* a research code written in C++
+* based on Dune (Distributed and Unified Numerics Environment)
+* a Dune user module in the Dune environment
+
+https://dumux.org/
+"""
+
+__import__("pkg_resources").declare_namespace(__name__)
diff --git a/python/dumux/assembly/CMakeLists.txt b/python/dumux/assembly/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..aba7cb7775de7b49adce1f7d1f0a1db3c68300db
--- /dev/null
+++ b/python/dumux/assembly/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_python_targets(assembly
+  __init__
+)
diff --git a/python/dumux/assembly/__init__.py b/python/dumux/assembly/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..b9e9ddcd3a9c3fdedf863ddd3c9c2050974744cf
--- /dev/null
+++ b/python/dumux/assembly/__init__.py
@@ -0,0 +1,56 @@
+"""Classes and function related to the assembly of linear systems"""
+
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
+
+
+@cppWrapperCreator
+def _createFVAssembler(*, problem, gridVariables, model, diffMethod="numeric", isImplicit=True):
+    """
+    Create an FVAssembler object
+
+    Args:
+        problem: A problem instance (boundary & initial conditions)
+        gridVariables: A grid variable instance (primary & secondary variables defined on a grid)
+        model: A DuMux model configuration instance
+        diffMethod (str): The method to compute derivatives of the residual (numeric or analytic)
+        isImplicit (bool): If the time discretization method is implicit or explicit
+
+    Returns:
+        An assembler object (Python-bindings of the corresponding C++ type)
+
+    Usage:
+        assembler = FVAssembler(
+            problem=problem, gridVariables=gridVars, model=model, diffMethod=diffMethod
+        )
+    """
+
+    if diffMethod == "numeric":
+        cppDiffMethod = "Dumux::DiffMethod::numeric"
+    elif diffMethod == "analytic":
+        cppDiffMethod = "Dumux::DiffMethod::analytic"
+    else:
+        raise ValueError(f"Unknown diffMethod {diffMethod}")
+
+    assemblerType = f"Dumux::FVAssembler<{model.cppType}, {cppDiffMethod}, {int(isImplicit)}>"
+    includes = (
+        problem._includes + problem.gridGeometry()._includes + ["dumux/assembly/fvassembler.hh"]
+    )
+    includes += ["dumux/python/assembly/fvassembler.hh"]
+
+    moduleName = "fvassembler_" + hashIt(assemblerType)
+    generator = SimpleGenerator("FVAssembler", "Dumux::Python")
+    module = generator.load(
+        includes,
+        assemblerType,
+        moduleName,
+        holder="std::shared_ptr",
+        preamble=model.cppHeader,
+    )
+    return module.FVAssembler(problem, problem.gridGeometry(), gridVariables)
+
+
+@cppWrapperClassAlias(creator=_createFVAssembler)
+class FVAssembler:
+    """Class alias used to instantiate an FVAssembler"""
diff --git a/python/dumux/common/CMakeLists.txt b/python/dumux/common/CMakeLists.txt
index dc0edac4c667c3749dba21046d8675d5bead2b08..ffc35a5823ad73956efb36ba72ae5434047f4913 100644
--- a/python/dumux/common/CMakeLists.txt
+++ b/python/dumux/common/CMakeLists.txt
@@ -1,5 +1,6 @@
 add_python_targets(common
   __init__
+  properties
 )
 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
index 23cdbfa5710456ff954c9281e27e7df725cb6ff4..872be2ffe593253705e7ee24b1d33ae9b183044c 100644
--- a/python/dumux/common/__init__.py
+++ b/python/dumux/common/__init__.py
@@ -1,21 +1,34 @@
-from ._common import *
+"""
+The DuMux common module
+containing classes and functions needed for most simulations
+"""
 
 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):
+from dumux.common.properties import Model, Property
+from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
+
+from ._common import *
+
+
+@cppWrapperCreator
+def _createFVProblemDecorator(gridGeometry, enableInternalDirichletConstraints=False):
+    """A problem decorator generator for Python problems
+
+    from dumux.common import FVProblem
+    @FVProblem(gridGeometry)
+    class MyProblem:
+        ...
+    """
 
     def createModule(numEq):
         priVarType = "Dune::FieldVector<double, {}>".format(numEq)
         ggType = gridGeometry._typeName
-        problemType = "Dumux::Python::FVProblem<{}, {}>".format(ggType, priVarType)
+        enableIntDirConstraint = "true" if enableInternalDirichletConstraints else "false"
+        problemType = "Dumux::Python::FVProblem<{}, {}, {}>".format(
+            ggType, priVarType, enableIntDirConstraint
+        )
         includes = gridGeometry._includes + ["dumux/python/common/fvproblem.hh"]
         moduleName = "fvproblem_" + hashIt(problemType)
         holderType = "std::shared_ptr<{}>".format(problemType)
@@ -23,18 +36,27 @@ def FVProblem(gridGeometry):
         module = generator.load(includes, problemType, moduleName, options=[holderType])
         return module
 
-    def FVProblemDecorator(Cls):
-        module = createModule(Cls.numEq)
+    def decorateFVProblem(cls):
+        module = createModule(cls.numEq)
+
         def createFVProblem():
-            return module.FVProblem(gridGeometry, Cls())
+            return module.FVProblem(gridGeometry, cls())
+
         return createFVProblem
 
-    return FVProblemDecorator
+    return decorateFVProblem
+
+
+@cppWrapperClassAlias(creator=_createFVProblemDecorator)
+class FVProblem:
+    """Class alias used to decorate a Python finite volume problem"""
 
 
-# Function for JIT copmilation of Dumux::BoundaryTypes
-def BoundaryTypes(numEq=1):
-    # only copmile this once per numEq
+@cppWrapperCreator
+def _createBoundaryTypes(numEq=1):
+    """Create BoundaryTypes instances"""
+
+    # only compile this once per numEq
     cacheKey = "BoundaryTypes_{}".format(numEq)
     try:
         return globals()[cacheKey]()
@@ -44,5 +66,10 @@ def BoundaryTypes(numEq=1):
         moduleName = "boundarytypes_" + hashIt(typeName)
         generator = SimpleGenerator("BoundaryTypes", "Dumux::Python")
         module = generator.load(includes, typeName, moduleName)
-        globals().update({cacheKey : module.BoundaryTypes})
+        globals().update({cacheKey: module.BoundaryTypes})
     return globals()[cacheKey]()
+
+
+@cppWrapperClassAlias(creator=_createBoundaryTypes)
+class BoundaryTypes:
+    """Class alias used to create a BoundaryTypes instance"""
diff --git a/python/dumux/common/_common.cc b/python/dumux/common/_common.cc
index 49bca345331814d7741894631e309edbd3c166f2..b1a27adffabdbdc77f070a3fd367b565f7d32b0b 100644
--- a/python/dumux/common/_common.cc
+++ b/python/dumux/common/_common.cc
@@ -18,10 +18,57 @@
  *****************************************************************************/
 
 #include <dune/python/pybind11/pybind11.h>
+#include <dune/python/pybind11/stl.h>
+
+#include <dumux/common/parameters.hh>
 #include <dumux/python/common/timeloop.hh>
 
 PYBIND11_MODULE(_common, module)
 {
+    using namespace Dumux;
+    using pybind11::operator""_a;
+
     // export time loop
-    Dumux::Python::registerTimeLoop<double>(module);
+    Python::registerTimeLoop<double>(module);
+
+    // initialize the parameters
+    Parameters::init();
+
+    // register parameter initializer for command line arguments
+    // Usage in Python: initParameters(sys.argv, {"bla": "blubb",})
+    module.def("initParameters", [](
+        std::vector<std::string>& argv,
+        pybind11::dict params,
+        const std::string& inputFileName
+    ){
+        // convert command line args to compatible format
+        std::vector<char *> argv_;
+        for (auto& arg : argv)
+            argv_.push_back(arg.data());
+
+        // convert command line args to parameter tree
+        const auto cmdParams = Parameters::parseCommandLine(argv_.size(), argv_.data());
+
+        const auto setParams = [&](auto& paramTree){
+            paramTree = cmdParams;
+            for (auto [key, value] : params)
+                paramTree[pybind11::cast<std::string>(key)] = pybind11::cast<std::string>(pybind11::str(value));
+        };
+
+        if (inputFileName.empty())
+            Parameters::init(setParams);
+        else
+            Parameters::init(inputFileName, setParams, /*input overwrites params*/false);
+
+    }, "argv"_a, "params"_a = pybind11::dict{}, "inputFileName"_a = "");
+
+    module.def("getParam", [](const std::string& name){
+        return getParam<std::string>(name);
+    }, "name"_a);
+
+    module.def("getParam", [](const std::string& name, const std::string& defaultValue){
+        return getParam<std::string>(name, defaultValue);
+    }, "name"_a, "default"_a);
+
+    module.def("printParameters", &Parameters::print);
 }
diff --git a/python/dumux/common/properties.py b/python/dumux/common/properties.py
new file mode 100644
index 0000000000000000000000000000000000000000..55d5d0a5ac62f731db36644aba8f5fb262f60469
--- /dev/null
+++ b/python/dumux/common/properties.py
@@ -0,0 +1,313 @@
+"""
+The DuMux property system in Python consisting of
+Property, TypeTag and Model
+"""
+
+import os
+from dataclasses import dataclass
+from typing import List, Union
+from dune.common.hashit import hashIt
+
+
+@dataclass
+class Property:
+    """
+    Properties are used to construct a model
+    Instances of Property are created with
+    Property.fromInstance(...), Property.fromCppType(...),
+    or  Property.fromValue(...).
+    """
+
+    cppType: str = None
+    cppIncludes: List[str] = None
+    requiredPropertyTypes: List[str] = None
+    requiredPropertyValues: List[str] = None
+    value: Union[bool, int, float] = None
+
+    @classmethod
+    def fromInstance(cls, inst):
+        """Create a Property from an instance of a wrapper"""
+        if not hasattr(inst, "_typeName"):
+            raise TypeError(
+                "The given instance {inst} does not have an attribute _typeName. "
+                "Only generated C++ objects (with Python bindings) are accepted."
+            )
+
+        return cls(cppType=inst._typeName, cppIncludes=inst._includes)
+
+    @classmethod
+    def fromCppType(
+        cls,
+        cppType: str,
+        *,
+        cppIncludes: List[str] = None,
+        requiredPropertyTypes: List[str] = None,
+        requiredPropertyValues: List[str] = None,
+    ):
+        """Create a Property from a given C++ type and includes"""
+        return cls(
+            cppType=cppType,
+            cppIncludes=cppIncludes,
+            requiredPropertyTypes=requiredPropertyTypes,
+            requiredPropertyValues=requiredPropertyValues,
+        )
+
+    @classmethod
+    def fromValue(cls, value):
+        """Create a Property from a given value"""
+        return cls(value=value)
+
+
+def typePropertyToString(propertyName, typeTagName, typeArg):
+    """Converts a Type Property to a string"""
+
+    propertyString = "template<class TypeTag>\n"
+    propertyString += "struct {}<TypeTag, TTag::{}>\n{{\n".format(propertyName, typeTagName)
+
+    if isinstance(typeArg, (float)):
+        propertyString += "    using type = {};\n".format("double")
+    elif isinstance(typeArg, (int)):
+        propertyString += "    using type = {};\n".format("int")
+    elif isinstance(typeArg, Property):
+        if typeArg.requiredPropertyTypes or typeArg.requiredPropertyValues:
+            propertyString += "private:\n"
+        if typeArg.requiredPropertyTypes is not None:
+            for reqProp in typeArg.requiredPropertyTypes:
+                propertyString += "    using {} = {};\n".format(
+                    reqProp, "GetPropType<TypeTag, Properties::{}>".format(reqProp)
+                )
+        if typeArg.requiredPropertyValues is not None:
+            for reqProp in typeArg.requiredPropertyValues:
+                reqPropLowerCase = reqProp[0].lower() + reqProp[1:]
+                propertyString += "    static constexpr auto {} = {};\n".format(
+                    reqPropLowerCase, "getPropValue<TypeTag, Properties::{}>()".format(reqProp)
+                )
+
+        propertyString += "public:\n"
+        propertyString += "    using type = {};\n".format(typeArg.cppType)
+
+    propertyString += "};"
+
+    return propertyString
+
+
+def valuePropertyToString(propertyName, typeTagName, value):
+    """Converts a Value Property to a string"""
+
+    propertyString = "template<class TypeTag>\n"
+    propertyString += "struct {}<TypeTag, TTag::{}>\n{{".format(propertyName, typeTagName)
+
+    # make sure to get the correct C++ types and values
+    if isinstance(value, bool):
+        value = str(value).lower()
+        cppType = "bool"
+    elif isinstance(value, int):
+        cppType = "int"
+    elif isinstance(value, float):
+        cppType = "Scalar"
+        propertyString += "\nprivate:\n"
+        propertyString += "    using Scalar = GetPropType<TypeTag, Properties::Scalar>;\n"
+        propertyString += "public:"
+    else:
+        raise ValueError(f"Invalid argument {value}. Expects bool, int or float.")
+
+    propertyString += "\n"
+    propertyString += f"    static constexpr {cppType} value = {value};"
+    propertyString += "\n"
+    propertyString += "};"
+
+    return propertyString
+
+
+_typeTags = {
+    "CCTpfaModel": {
+        "include": "dumux/discretization/cctpfa.hh",
+        "description": "A cell-centered two-point flux finite volume discretization scheme.",
+    },
+    "BoxModel": {
+        "include": "dumux/discretization/box.hh",
+        "description": "A node-centered two-point flux finite volume discretization scheme",
+    },
+    "OneP": {
+        "include": "dumux/porousmediumflow/1p/model.hh",
+        "description": "A model for single-phase flow in porous media.",
+    },
+}
+
+
+def listTypeTags():
+    """List all available TypeTags/Models that can be inherited from"""
+
+    print("\n**********************************\n")
+    print("The following TypeTags are availabe:")
+    for key, value in _typeTags.items():
+        print(key, ":", value["description"])
+    print("\n**********************************")
+
+
+def predefinedProperties():
+    """Create a list of properties defined in properties.hh"""
+
+    propertiesHeader = os.path.abspath(
+        os.path.dirname(__file__) + "/../../../../dumux/common/properties.hh"
+    )
+    with open(propertiesHeader) as header:
+        properties = []
+        for line in header:
+            if line.startswith("struct"):
+                properties.append(line.split(" ")[1])
+        return properties
+
+
+class TypeTag:
+    """TypeTags are inheritable collections of properties"""
+
+    knownProperties = predefinedProperties()
+
+    def __init__(self, name, *, inheritsFrom=None, gridGeometry=None):
+        self.inheritsFrom = inheritsFrom
+        self.includes = []
+        self.properties = {}
+        self.newPropertyDefinitions = []
+        self.gridGeometry = gridGeometry
+        self.name = name
+
+        if self.name in _typeTags.keys():
+            if inheritsFrom is not None:
+                raise ValueError(
+                    f"Existing TypeTag {name} cannot inherit from other TypeTags."
+                    f" Use TypeTag({name}) only."
+                )
+            self.isExistingTypeTag = True
+            self.includes = [_typeTags[self.name]["include"]]
+        else:
+            self.isExistingTypeTag = False
+
+        if self.inheritsFrom is not None:
+            # treat existing TypeTags by converting the given string to a real TypeTag instance
+            for idx, parentTypeTag in enumerate(self.inheritsFrom):
+                if not isinstance(parentTypeTag, TypeTag):
+                    if not isinstance(parentTypeTag, str):
+                        raise ValueError(
+                            "Unknown parent TypeTag {}. Use either argument of type TypeTag "
+                            "or a string for an existing TypeTag. "
+                            "List of existing TypeTags: {}".format(parentTypeTag, _typeTags.keys())
+                        )
+                    if parentTypeTag not in _typeTags.keys():
+                        raise ValueError(
+                            "Unknown TypeTag {}. List of existing TypeTags: {}".format(
+                                parentTypeTag, _typeTags.keys()
+                            )
+                        )
+                    self.inheritsFrom[idx] = TypeTag(parentTypeTag)
+
+            # pick up the properties and includes of the parent TypeTag
+            for parentTypeTag in reversed(self.inheritsFrom):
+                self.newPropertyDefinitions += parentTypeTag.newPropertyDefinitions
+                for key in parentTypeTag.properties:
+                    self.properties[key] = parentTypeTag.properties[key]
+                if parentTypeTag.includes is not None:
+                    for include in parentTypeTag.includes:
+                        self.includes.append(include)
+
+    def __setitem__(self, key, value: Property):
+        """the [] operator for setting values"""
+        if not isinstance(value, Property):
+            raise ValueError("Only values of type Property can be assigned to a model")
+
+        if key not in self.knownProperties:
+            print(f"Adding {key} as new property")
+            self.newPropertyDefinitions += [key]
+
+        self.properties[key] = value
+        if value.cppIncludes is not None:
+            for include in value.cppIncludes:
+                self.includes.append(include)
+
+    def __getitem__(self, key):
+        """the [] operator for getting values"""
+        return self.properties[key]
+
+    @property
+    def cppType(self):
+        """Returns the TypeTag as a string"""
+        return "Dumux::Properties::TTag::" + self.name
+
+    @property
+    def cppHeader(self):
+        """creates a string resembling a properties.hh file"""
+
+        file = "#ifndef DUMUX_{}_PROPERTIES_HH\n".format(self.name.upper())
+        file += "#define DUMUX_{}_PROPERTIES_HH\n\n".format(self.name.upper())
+
+        file += "#include <dumux/common/properties.hh>\n"
+
+        for include in self.includes:
+            assert "<" not in include
+            file += "#include <{}>\n".format(include)
+
+        file += "\n"
+
+        file += "namespace Dumux::Properties {\n\n"
+        file += "namespace TTag {\n"
+
+        if self.inheritsFrom is not None:
+            for otherTypeTag in self.inheritsFrom:
+                if not otherTypeTag.isExistingTypeTag:
+                    file += "struct {}{{}}; \n".format(otherTypeTag.name)
+
+            args = ",".join(x.name for x in self.inheritsFrom)
+        else:
+            args = ""
+
+        file += f"struct {self.name}\n{{\n"
+        file += f"    using InheritsFrom = std::tuple<{args}>;\n"
+        if self.gridGeometry is not None:
+            file += f"    using GridGeometry = {self.gridGeometry._typeName};\n"
+        file += "};\n} // end namespace TTag\n\n"
+
+        for newDef in self.newPropertyDefinitions:
+            file += "template<class TypeTag, class MyTypeTag>\n"
+            file += "struct " + newDef + " { using type =  UndefinedProperty; };\n\n"
+
+        for prop in self.properties:
+            if self[prop].value is not None:
+                file += valuePropertyToString(prop, self.name, self[prop].value) + "\n\n"
+            else:
+                file += typePropertyToString(prop, self.name, self[prop]) + "\n\n"
+
+        if self.gridGeometry is not None:
+            file += (
+                typePropertyToString(
+                    "Grid", self.name, Property.fromCppType("typename TypeTag::GridGeometry::Grid")
+                )
+                + "\n\n"
+            )
+
+        file += "} // end namespace Dumux::Properties \n\n"
+        file += "#endif"
+
+        return file
+
+
+class Model(TypeTag):
+    """A DuMux model specifies all properties necessary to build a DuMux simulator"""
+
+    def __init__(self, *, inheritsFrom: List[str], gridGeometry, scalar: str = "double"):
+        # generate a generic name
+        genericName = "TypeTag" + hashIt("".join(inheritsFrom) + gridGeometry._typeName + scalar)
+
+        # deduce the discretization tag from the grid geometry
+        discretizationMethod = gridGeometry.discMethod
+        discretizationMap = {
+            "box": "BoxModel",
+            "cctpfa": "CCTpfaModel",
+        }
+        if discretizationMethod in discretizationMap:
+            inheritsFrom += [discretizationMap[discretizationMethod]]
+
+        # call parent constructor
+        super().__init__(name=genericName, inheritsFrom=inheritsFrom, gridGeometry=gridGeometry)
+
+        # set the scalar type
+        self.__setitem__("Scalar", Property.fromCppType(scalar))
diff --git a/python/dumux/discretization/__init__.py b/python/dumux/discretization/__init__.py
index 238240b2ac2d44089298b8e2d159c266f8ac406d..f9ab3d47ec0839603a8e5f41ed462a7e55f64aea 100644
--- a/python/dumux/discretization/__init__.py
+++ b/python/dumux/discretization/__init__.py
@@ -1,9 +1,14 @@
+"""Classes and functions related to discretization methods"""
+
 from dune.generator.generator import SimpleGenerator
 from dune.common.hashit import hashIt
+from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
+
+
+@cppWrapperCreator
+def _createGridGeometry(gridView, discMethod="cctpfa"):
+    """Construct a GridGeometry from a gridView"""
 
-# 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":
@@ -20,3 +25,36 @@ def GridGeometry(gridView, discMethod="cctpfa"):
     generator = SimpleGenerator("GridGeometry", "Dumux::Python")
     module = generator.load(includes, typeName, moduleName, options=[holderType])
     return module.GridGeometry(gridView)
+
+
+@cppWrapperClassAlias(creator=_createGridGeometry)
+class GridGeometry:
+    """Class alias used to instantiate a GridGeometry"""
+
+
+@cppWrapperCreator
+def _createGridVariables(*, problem, model):
+    """Construct a GridGeometry from problem and the model"""
+
+    includes = [
+        "dumux/discretization/fvgridvariables.hh",
+        "dumux/python/discretization/gridvariables.hh",
+    ]
+    ggeo = f"Dumux::GetPropType<{model.cppType}, Dumux::Properties::GridGeometry>"
+    gvv = f"Dumux::GetPropType<{model.cppType}, Dumux::Properties::GridVolumeVariables>"
+    gfc = f"Dumux::GetPropType<{model.cppType}, Dumux::Properties::GridFluxVariablesCache>"
+    typeName = "Dumux::FVGridVariables<{}, {}, {}>".format(ggeo, gvv, gfc)
+
+    moduleName = "gridvariables_" + hashIt(typeName)
+    holderType = "std::shared_ptr<{}>".format(typeName)
+    generator = SimpleGenerator("GridVariables", "Dumux::Python")
+    module = generator.load(
+        includes, typeName, moduleName, options=[holderType], preamble=model.cppHeader
+    )
+    module.GridVariables.model = property(lambda self: model)
+    return module.GridVariables(problem, problem.gridGeometry())
+
+
+@cppWrapperClassAlias(creator=_createGridVariables)
+class GridVariables:
+    """Class alias used to instantiate GridVariables"""
diff --git a/python/dumux/io/CMakeLists.txt b/python/dumux/io/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f4c5c0ef90f208f15f91b450b7c616be4b668233
--- /dev/null
+++ b/python/dumux/io/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_python_targets(io
+  __init__
+)
diff --git a/python/dumux/io/__init__.py b/python/dumux/io/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..c92e22f9e63bf43f96cead7c71cff2e1bda05a8d
--- /dev/null
+++ b/python/dumux/io/__init__.py
@@ -0,0 +1,25 @@
+"""DuMux input-output library"""
+
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
+
+
+@cppWrapperCreator
+def _createVtkOutputModule(*, gridVariables, solutionVector, name):
+    """Construct a VtkOutputModule"""
+
+    includes = gridVariables._includes + solutionVector._includes
+    includes += ["dumux/python/io/vtkoutputmodule.hh", "dumux/io/vtkoutputmodule.hh"]
+    typeName = "Dumux::VtkOutputModule<{}, {}>".format(
+        gridVariables._typeName, solutionVector._typeName
+    )
+    moduleName = "vtkoutputmodule_" + hashIt(typeName)
+    generator = SimpleGenerator("VtkOutputModule", "Dumux::Python")
+    module = generator.load(includes, typeName, moduleName, preamble=gridVariables.model.cppHeader)
+    return module.VtkOutputModule(gridVariables, solutionVector, name)
+
+
+@cppWrapperClassAlias(creator=_createVtkOutputModule)
+class VtkOutputModule:
+    """Class alias used to instantiate a VtkOutputModule"""
diff --git a/python/dumux/material/CMakeLists.txt b/python/dumux/material/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c6683fdb0b9fc87f479c79f405fa95c1e2703360
--- /dev/null
+++ b/python/dumux/material/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_subdirectory(components)
+add_subdirectory(fluidsystems)
+add_subdirectory(spatialparams)
+
+add_python_targets(material
+  __init__
+)
diff --git a/python/dumux/material/__init__.py b/python/dumux/material/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..01534f242281cf3b905f3a0a44349f4076050a0d
--- /dev/null
+++ b/python/dumux/material/__init__.py
@@ -0,0 +1,5 @@
+"""DuMux material framework providing constitutive laws"""
+
+from dumux.material.fluidsystems import FluidSystem
+from dumux.material.components import Component
+from dumux.material.spatialparams import OnePSpatialParams
diff --git a/python/dumux/material/components/CMakeLists.txt b/python/dumux/material/components/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e02d4be73abc3753a51c9aadb95d408b2578d270
--- /dev/null
+++ b/python/dumux/material/components/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_python_targets(components
+  __init__
+)
diff --git a/python/dumux/material/components/__init__.py b/python/dumux/material/components/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..12b1989a070419c3eace1e3baa12814f0b99870f
--- /dev/null
+++ b/python/dumux/material/components/__init__.py
@@ -0,0 +1,38 @@
+"""Components are building blocks of fluid and solid systems"""
+
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
+
+
+_components = {
+    "Air": "dumux/material/components/air.hh",
+    "H2O": "dumux/material/components/h2o.hh",
+    "SimpleH2O": "dumux/material/components/simpleh2o.hh",
+    "N2": "dumux/material/components/n2.hh",
+    "Calcite": "dumux/material/components/calcite.hh",
+}
+
+
+def listComponents():
+    """List all available component names"""
+    print("The following components are availabe:")
+    print(_components.keys())
+
+
+@cppWrapperCreator
+def _createComponent(name, *, scalar="double"):
+    """Create a new component of the given name"""
+
+    typeName = f"Dumux::Components::{name} <{scalar}>"
+    moduleName = f"{name.lower()}_{hashIt(typeName)}"
+    includes = ["dumux/python/material/components/component.hh"]
+    includes += [_components[name]]
+    generator = SimpleGenerator("Component", "Dumux::Python")
+    module = generator.load(includes, typeName, moduleName)
+    return module.Component()
+
+
+@cppWrapperClassAlias(creator=_createComponent)
+class Component:
+    """Class alias used to instantiate a Component"""
diff --git a/python/dumux/material/fluidsystems/CMakeLists.txt b/python/dumux/material/fluidsystems/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..172368eadc24dffd751aff799c0b44341e5cc590
--- /dev/null
+++ b/python/dumux/material/fluidsystems/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_python_targets(fluidsystems
+  __init__
+)
diff --git a/python/dumux/material/fluidsystems/__init__.py b/python/dumux/material/fluidsystems/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..29c73efd8dcb6426cb0da6ef5f251175106b22d3
--- /dev/null
+++ b/python/dumux/material/fluidsystems/__init__.py
@@ -0,0 +1,46 @@
+"""Fluid systems"""
+
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+from dumux.common import Property
+from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
+
+
+@cppWrapperCreator
+def _createFluidSystem(name, *, component, scalar: Property = None):
+    """Construct a FluidSystem"""
+
+    includes = component._includes + ["dumux/python/material/fluidsystems/fluidsystem.hh"]
+
+    if scalar is None:
+        scalar = Property.fromCppType("double")
+    scalarType = scalar.cppType
+
+    availableFluidSystems = {
+        "OnePLiquid": {
+            "includes": ["dumux/material/fluidsystems/1pliquid.hh"],
+            "type": f"Dumux::FluidSystems::OnePLiquid<{scalarType}, {component._typeName}>",
+        },
+        "OnePGas": {
+            "includes": ["dumux/material/fluidsystems/1pgas.hh"],
+            "type": f"Dumux::FluidSystems::OnePGas<{scalarType}, {component._typeName}>",
+        },
+    }
+    if name not in availableFluidSystems:
+        raise NotImplementedError(
+            "FluidSystem of type " + name + " not implemented.\n"
+            "Available types are " + ", ".join(availableFluidSystems.keys())
+        )
+
+    includes += availableFluidSystems[name]["includes"]
+    typeName = availableFluidSystems[name]["type"]
+
+    moduleName = "fluidsystem_" + hashIt(typeName)
+    generator = SimpleGenerator("FluidSystem", "Dumux::Python")
+    module = generator.load(includes, typeName, moduleName)
+    return module.FluidSystem()
+
+
+@cppWrapperClassAlias(creator=_createFluidSystem)
+class FluidSystem:
+    """Class alias used to instantiate a FluidSystem"""
diff --git a/python/dumux/material/spatialparams/CMakeLists.txt b/python/dumux/material/spatialparams/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..db81731a161e443a1c2c0f686cbf84d16af1d2cb
--- /dev/null
+++ b/python/dumux/material/spatialparams/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_python_targets(spatialparams
+  __init__
+)
diff --git a/python/dumux/material/spatialparams/__init__.py b/python/dumux/material/spatialparams/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..3440a500027b40e53b87c2b2a77d58c802a83e9f
--- /dev/null
+++ b/python/dumux/material/spatialparams/__init__.py
@@ -0,0 +1,77 @@
+"""Spatial parameter classes"""
+
+import numpy as np
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+from dumux.common import Property
+from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
+
+
+@cppWrapperCreator
+def _createOnePSpatialParamsDecorator(*, gridGeometry, scalar: Property = None):
+    """Turn a Python spatial parameter class into an C++/Python hybrid class"""
+
+    dim = gridGeometry.gridView.dimension
+    includes = gridGeometry._includes + [
+        "dumux/python/material/spatialparams/spatialparams.hh",
+        "dune/common/fmatrix.hh",
+    ]
+
+    if scalar is None:
+        scalar = Property.fromCppType("double")
+    scalarType = scalar.cppType
+
+    permType = f"Dune::FieldMatrix<{scalarType}, {dim}, {dim}>"
+    typeName = (
+        f"Dumux::Python::FVSpatialParamsOneP<{gridGeometry._typeName}, {scalarType}, {permType}>"
+    )
+    moduleName = f"spatialparams_{hashIt(typeName)}"
+
+    def decorateOnePSpatialParams(cls):
+        generator = SimpleGenerator("OnePSpatialParams", "Dumux::Python")
+        module = generator.load(includes, typeName, moduleName, holder="std::shared_ptr")
+
+        def maybeConvertScalarToMatrix(permeabilityValue):
+            if isinstance(permeabilityValue, float):
+                matrix = np.zeros(shape=(dim, dim))
+                np.fill_diagonal(matrix, permeabilityValue)
+                return matrix.tolist()
+
+            return permeabilityValue
+
+        class Permeability:
+            """Permeability decorator to make sure permeability has correct type"""
+
+            def __init__(self, permeabilityFunction):
+                self.permeabilityFunction = permeabilityFunction
+
+            def __call__(self, element, scv, elemSol):
+                result = self.permeabilityFunction(element, scv, elemSol)
+                return maybeConvertScalarToMatrix(result)
+
+        class PermeabilityAtPos:
+            """PermeabilityAtPos decorator to make sure permeability has correct type"""
+
+            def __init__(self, permeabilityFunction):
+                self.permeabilityFunction = permeabilityFunction
+
+            def __call__(self, globalPos):
+                result = self.permeabilityFunction(globalPos)
+                return maybeConvertScalarToMatrix(result)
+
+        def createSpatialParams():
+            spatialParams = cls()
+            if hasattr(cls, "permeability"):
+                cls.permeability = Permeability(spatialParams.permeability)
+            if hasattr(cls, "permeabilityAtPos"):
+                cls.permeabilityAtPos = PermeabilityAtPos(spatialParams.permeabilityAtPos)
+            return module.OnePSpatialParams(gridGeometry, spatialParams)
+
+        return createSpatialParams
+
+    return decorateOnePSpatialParams
+
+
+@cppWrapperClassAlias(creator=_createOnePSpatialParamsDecorator)
+class OnePSpatialParams:
+    """Class alias used to decorate Python spatial parameter implementations"""
diff --git a/python/dumux/porousmediumflow/CMakeLists.txt b/python/dumux/porousmediumflow/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9b6c72a3884109241688de5f2cd672285b153f04
--- /dev/null
+++ b/python/dumux/porousmediumflow/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_python_targets(porousmediumflow
+  __init__
+)
diff --git a/python/dumux/porousmediumflow/__init__.py b/python/dumux/porousmediumflow/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..15e8eed9ab3ee0f68422cc5024757f5e8be37bd1
--- /dev/null
+++ b/python/dumux/porousmediumflow/__init__.py
@@ -0,0 +1,81 @@
+"""Classes and functions related to the porousmedium flow models"""
+
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
+
+
+@cppWrapperCreator
+def _createPorousMediumFlowProblemDecorator(
+    gridGeometry, spatialParams, enableInternalDirichletConstraints=False
+):
+    """A problem decorator generator for Python problems
+
+    Usage:
+        from dumux.common import PorousMediumFlowProblem
+        @PorousMediumFlowProblem(gridGeometry)
+        class MyProblem:
+            ...
+    """
+
+    def createModule(numEq):
+        priVars = f"Dune::FieldVector<double, {numEq}>"
+        ggType = gridGeometry._typeName
+        spType = spatialParams._typeName
+        enableIDC = "true" if enableInternalDirichletConstraints else "false"
+        problemType = (
+            "Dumux::Python::PorousMediumFlowProblem" f"<{ggType}, {priVars}, {spType}, {enableIDC}>"
+        )
+        includes = (
+            gridGeometry._includes
+            + spatialParams._includes
+            + ["dumux/python/porousmediumflow/problem.hh"]
+        )
+        moduleName = "fvproblem_" + hashIt(problemType)
+        generator = SimpleGenerator("PorousMediumFlowProblem", "Dumux::Python")
+        module = generator.load(includes, problemType, moduleName, holder="std::shared_ptr")
+        return module
+
+    def decoratePorousMediumFlowProblem(cls):
+        module = createModule(cls.numEq)
+
+        def createPorousMediumFlowProblem():
+            return module.PorousMediumFlowProblem(gridGeometry, spatialParams, cls())
+
+        return createPorousMediumFlowProblem
+
+    return decoratePorousMediumFlowProblem
+
+
+@cppWrapperClassAlias(creator=_createPorousMediumFlowProblemDecorator)
+class PorousMediumFlowProblem:
+    """A class alias used to create a problem decorator Python problems"""
+
+
+@cppWrapperCreator
+def _createPorousMediumFlowVelocityOutput(*, gridVariables):
+    """Create a PorousMediumFlowVelocityOutput"""
+
+    includes = gridVariables._includes
+    includes += ["dumux/python/porousmediumflow/velocityoutput.hh", "dumux/io/velocityoutput.hh"]
+    fluxVarsType = (
+        f"Dumux::GetPropType<{gridVariables.model.cppType}, Dumux::Properties::FluxVariables>"
+    )
+    typeName = f"Dumux::PorousMediumFlowVelocityOutput<{gridVariables._typeName}, {fluxVarsType}>"
+    moduleName = "porousmediumflowvelocityoutput_" + hashIt(typeName)
+    baseClass = [f"Dumux::VelocityOutput<{gridVariables._typeName}>"]
+    generator = SimpleGenerator("PorousMediumFlowVelocityOutput", "Dumux::Python")
+    module = generator.load(
+        includes,
+        typeName,
+        moduleName,
+        holder="std::shared_ptr",
+        preamble=gridVariables.model.cppHeader,
+        baseClasses=baseClass,
+    )
+    return module.PorousMediumFlowVelocityOutput(gridVariables)
+
+
+@cppWrapperClassAlias(creator=_createPorousMediumFlowVelocityOutput)
+class PorousMediumFlowVelocityOutput:
+    """A class alias used to create PorousMediumFlowVelocityOutput instances"""
diff --git a/python/dumux/wrapping/CMakeLists.txt b/python/dumux/wrapping/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9b2e59b4e0914fdfc0524f75fa839238e4dcb3e5
--- /dev/null
+++ b/python/dumux/wrapping/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_python_targets(wrapping
+  __init__
+)
diff --git a/python/dumux/wrapping/__init__.py b/python/dumux/wrapping/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..c32b764f0582ca29864f5ef0ebd241fe7396f48e
--- /dev/null
+++ b/python/dumux/wrapping/__init__.py
@@ -0,0 +1,40 @@
+"""Helper classes and function related Python <-> C++ interaction"""
+
+import functools
+from typing import Callable, Type
+
+
+def cppWrapperCreator(creator: Callable) -> Callable:
+    """
+    Decorator for creator functions that return a C++ type with Python bindings
+    resulting from C++ code generation and just-in-time compilation
+    """
+
+    def makeCreator(aliasClass: Type) -> Callable:
+        @functools.wraps(creator)
+        def _makeCreator(*args, **kwargs):
+            return creator(*args, **kwargs)
+
+        # make the creator assume the name of the alias class
+        _makeCreator.__name__ = aliasClass.__name__
+        return _makeCreator
+
+    return makeCreator
+
+
+def cppWrapperClassAlias(creator: Callable) -> Callable:
+    """
+    Decorator for a class alias corresponding to a creator function
+    This makes the creator function appear like a class constructor
+
+    Args:
+        creator (Callable): The corresponding creator function the
+                            decorated class is an alias for
+    Returns:
+        Callable: The creator function with the alias name
+    """
+
+    def makeCreator(aliasClass: Type) -> Callable:
+        return creator(aliasClass)
+
+    return makeCreator
diff --git a/python/setup-python-env.sh b/python/setup-python-env.sh
new file mode 100755
index 0000000000000000000000000000000000000000..f99c549c7abd1844ab413e45a81e67a5048d2f96
--- /dev/null
+++ b/python/setup-python-env.sh
@@ -0,0 +1,7 @@
+#!/bin/bash
+# Adds all Python modules found in other Dune modules to the PYTHONPATH
+./dune-common/bin/dunecontrol bexec "echo -n :\$(pwd)/python >> $(pwd)/pythonpath.txt"
+export PYTHONPATH=$PYTHONPATH$(cat pythonpath.txt)
+rm pythonpath.txt
+# Sets up the dune-py module for JIT compilation of Python binding code
+./dune-common/bin/setup-dunepy.py --opts=cmake.opts install
diff --git a/python/setup.py.in b/python/setup.py.in
index b59c513c242be131aa68204d3bd85c45016279cd..970d78de29b51be6fba1af93e6b12acec5f4f350 100644
--- a/python/setup.py.in
+++ b/python/setup.py.in
@@ -1,12 +1,13 @@
 from setuptools import setup, find_namespace_packages
 
-setup(name="${ProjectName}",
-      description="${ProjectDescription}",
-      version="${ProjectVersionString}",
-      author="${ProjectAuthor}",
-      author_email="${ProjectMaintainerEmail}",
-      packages = find_namespace_packages(include=['dumux.*']),
-      zip_safe = 0,
-      package_data = {'': ['*.so']},
-      install_requires = "${ProjectPythonRequires}".split(' ')
+setup(
+    name="${ProjectName}",
+    description="${ProjectDescription}",
+    version="${ProjectVersionString}",
+    author="${ProjectAuthor}",
+    author_email="${ProjectMaintainerEmail}",
+    packages=find_namespace_packages(include=["dumux.*"]),
+    zip_safe=0,
+    package_data={"": ["*.so"]},
+    install_requires="${ProjectPythonRequires}".split(" "),
 )
diff --git a/test/python/CMakeLists.txt b/test/python/CMakeLists.txt
index 4526eb14bf40c6220ae91d16a6c1042517ca0d15..43d416d6f44d465fe65f9cfbf250bc43953b9d90 100644
--- a/test/python/CMakeLists.txt
+++ b/test/python/CMakeLists.txt
@@ -1,43 +1,64 @@
 # We execute Python tests in the binary directory so possible output files
 # do not get dumped into the source tree
 
-if(${dune-common_VERSION} VERSION_LESS 2.8)
-  dune_python_add_test(NAME test_python_gridgeometry
-                       COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_gridgeometry.py
-                       WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-                       LABELS python unit)
+dune_python_add_test(NAME test_python_gridgeometry
+                     SCRIPT ${CMAKE_CURRENT_SOURCE_DIR}/test_gridgeometry.py
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+                     LABELS python unit)
+set_tests_properties(test_python_gridgeometry PROPERTIES TIMEOUT 1200)
 
-  dune_python_add_test(NAME test_python_fvproblem
-                       COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_fvproblem.py
-                       WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-                       LABELS python unit)
+dune_python_add_test(NAME test_python_fvproblem
+                     SCRIPT ${CMAKE_CURRENT_SOURCE_DIR}/test_fvproblem.py
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+                     LABELS python unit)
+set_tests_properties(test_python_fvproblem PROPERTIES TIMEOUT 1200)
 
-  dune_python_add_test(NAME test_python_explicit_transport_cctpfa
-                       LABELS python
-                       WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-                       COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
-                               --script fuzzy
-                               --files ${CMAKE_SOURCE_DIR}/test/references/test_python_explicit-transport-cctpfa-100.vtu
-                                       ${CMAKE_CURRENT_BINARY_DIR}/finitevolume-solution-100.vtu
-                               --command "${CMAKE_CURRENT_SOURCE_DIR}/test_explicit_transport_cctpfa.py")
-else()
-  dune_python_add_test(NAME test_python_gridgeometry
-                       SCRIPT ${CMAKE_CURRENT_SOURCE_DIR}/test_gridgeometry.py
-                       WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-                       LABELS python unit)
-
-  dune_python_add_test(NAME test_python_fvproblem
-                       SCRIPT ${CMAKE_CURRENT_SOURCE_DIR}/test_fvproblem.py
-                       WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-                       LABELS python unit)
-
-  dune_python_add_test(NAME test_python_explicit_transport_cctpfa
-                       LABELS python
-                       WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-                       SCRIPT ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+dune_python_add_test(NAME test_python_explicit_transport_cctpfa
+                     LABELS python
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+                     SCRIPT ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
                               --script fuzzy
                               --files ${CMAKE_SOURCE_DIR}/test/references/test_python_explicit-transport-cctpfa-100.vtu
                                       ${CMAKE_CURRENT_BINARY_DIR}/finitevolume-solution-100.vtu
                               --command "${CMAKE_CURRENT_SOURCE_DIR}/test_explicit_transport_cctpfa.py")
+set_tests_properties(test_python_explicit_transport_cctpfa PROPERTIES TIMEOUT 1200)
+
+dune_python_add_test(NAME test_python_1p_incompressible_box_numdiff
+                     LABELS python
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+                     SCRIPT ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+                                --script fuzzy
+                                --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_box-reference.vtu
+                                        ${CMAKE_CURRENT_BINARY_DIR}/test_1p_box_numeric-00000.vtu
+                                --command "${CMAKE_CURRENT_SOURCE_DIR}/test_1p.py -DiscMethod box -DiffMethod numeric")
+set_tests_properties(test_python_1p_incompressible_box_numdiff PROPERTIES TIMEOUT 1200)
+
+dune_python_add_test(NAME test_python_1p_incompressible_box_anadiff
+                     LABELS python
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+                     SCRIPT ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+                                --script fuzzy
+                                --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_box-reference.vtu
+                                        ${CMAKE_CURRENT_BINARY_DIR}/test_1p_box_analytic-00000.vtu
+                                --command "${CMAKE_CURRENT_SOURCE_DIR}/test_1p.py -DiscMethod box -DiffMethod analytic")
+set_tests_properties(test_python_1p_incompressible_box_anadiff PROPERTIES TIMEOUT 1200)
+
+dune_python_add_test(NAME test_python_1p_incompressible_cctpfa_numdiff
+                     LABELS python
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+                     SCRIPT ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+                                --script fuzzy
+                                --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_cc-reference.vtu
+                                        ${CMAKE_CURRENT_BINARY_DIR}/test_1p_cctpfa_numeric-00000.vtu
+                                --command "${CMAKE_CURRENT_SOURCE_DIR}/test_1p.py -DiscMethod cctpfa -DiffMethod numeric")
+set_tests_properties(test_python_1p_incompressible_cctpfa_numdiff PROPERTIES TIMEOUT 1200)
 
-endif()
+dune_python_add_test(NAME test_python_1p_incompressible_cctpfa_anadiff
+                     LABELS python
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+                     SCRIPT ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+                                --script fuzzy
+                                --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_cc-reference.vtu
+                                        ${CMAKE_CURRENT_BINARY_DIR}/test_1p_cctpfa_analytic-00000.vtu
+                                --command "${CMAKE_CURRENT_SOURCE_DIR}/test_1p.py -DiscMethod cctpfa -DiffMethod analytic")
+set_tests_properties(test_python_1p_incompressible_cctpfa_anadiff PROPERTIES TIMEOUT 1200)
diff --git a/test/python/test_1p.py b/test/python/test_1p.py
new file mode 100755
index 0000000000000000000000000000000000000000..80ff3088c16cd59e2ac085371dcc5716d01dd257
--- /dev/null
+++ b/test/python/test_1p.py
@@ -0,0 +1,162 @@
+#!/usr/bin/env python3
+
+import sys
+
+from dune.grid import structuredGrid
+from dune.istl import blockVector, CGSolver, SeqSSOR
+
+from dumux.common import initParameters, printParameters, getParam
+from dumux.common import BoundaryTypes, Model, Property
+from dumux.discretization import GridGeometry, GridVariables
+from dumux.assembly import FVAssembler
+from dumux.porousmediumflow import PorousMediumFlowProblem, PorousMediumFlowVelocityOutput
+from dumux.material import FluidSystem, Component, OnePSpatialParams
+from dumux.io import VtkOutputModule
+
+# Initialize parameters
+initParameters(
+    argv=sys.argv,
+    params={
+        "Problem.EnableGravity": True,
+        "Vtk.AddVelocity": False,
+        "Assembly.NumericDifference.PriVarMagnitude": 1e5,
+    },
+)
+
+discMethod = getParam("DiscMethod", default="box")
+if discMethod not in ["box", "cctpfa"]:
+    raise NotImplementedError(discMethod + " not supported yet. Use cctpfa or box")
+
+diffMethod = getParam("DiffMethod", default="numeric")
+if diffMethod not in ["analytic", "numeric"]:
+    raise NotImplementedError(diffMethod + " must be analytic or numeric")
+
+
+# Set up the grid and the grid geometry
+gridView = structuredGrid([0, 0], [1, 1], [10, 10])
+gridGeometry = GridGeometry(gridView=gridView, discMethod=discMethod)
+
+# Set up the model
+model = Model(inheritsFrom=["OneP"], gridGeometry=gridGeometry)
+
+# Tell model to use a particular local residual type
+model["LocalResidual"] = Property.fromCppType(
+    "OnePIncompressibleLocalResidual<TypeTag>",
+    cppIncludes=["dumux/porousmediumflow/1p/incompressiblelocalresidual.hh"],
+)
+
+# Setup the fluid system
+h20 = Component("SimpleH2O")
+onePLiquid = FluidSystem("OnePLiquid", component=h20, scalar=model["Scalar"])
+model["FluidSystem"] = Property.fromInstance(onePLiquid)
+
+
+# Define the spatial parameters
+@OnePSpatialParams(gridGeometry=gridGeometry)
+class SpatialParams:
+    dimWorld = gridGeometry.gridView.dimWorld
+    lensLowerLeft = [0.2, 0.2]
+    lensUpperRight = [0.8, 0.8]
+
+    def isLens(self, globalPos):
+        eps = 1.5e-7
+        for i in range(self.dimWorld):
+            if (globalPos[i] < self.lensLowerLeft[i] + eps) or (
+                globalPos[i] > self.lensUpperRight[i] - eps
+            ):
+                return False
+        return True
+
+    def permeability(self, element, scv, elemSol):
+        globalPos = scv.dofPosition
+        # permeability can be either given
+        # as scalar or tensor
+        if self.isLens(globalPos):
+            return [[1e-12, 0], [0, 1e-12]]
+
+        return 1e-10
+
+    def porosityAtPos(self, globalPos):
+        return 0.4
+
+
+# and set them as a model property
+spatialParams = SpatialParams()
+model["SpatialParams"] = Property.fromInstance(spatialParams)
+
+
+# Define the problem
+@PorousMediumFlowProblem(gridGeometry, spatialParams)
+class Problem:
+    numEq = 1
+
+    def boundaryTypesAtPos(self, globalPos):
+        bTypes = BoundaryTypes(self.numEq)
+        eps = 1e-6
+
+        if globalPos[1] < eps or globalPos[1] > gridGeometry.bBoxMax[1] - eps:
+            bTypes.setDirichlet()
+        else:
+            bTypes.setNeumann()
+
+        return bTypes
+
+    def dirichletAtPos(self, globalPos):
+        dp_dy_ = -1.0e5
+        return 1.0e5 + dp_dy_ * (globalPos[1] - gridGeometry.bBoxMax[1])
+
+    def sourceAtPos(self, globalPos):
+        return 0.0
+
+    def extrusionFactor(self, element, scv):
+        return 1.0
+
+    def temperatureAtPos(self, globalPos):
+        return 283.15
+
+    def name(self):
+        return "python_problem"
+
+    def paramGroup(self):
+        return ""
+
+    def neumann(self, element, fvGeometry, scvf):
+        return 0
+
+
+# and set it as a model property
+problem = Problem()
+model["Problem"] = Property.fromInstance(problem)
+
+# Initialize the GridVariables and the Assembler
+gridVars = GridVariables(problem=problem, model=model)
+assembler = FVAssembler(problem=problem, gridVariables=gridVars, model=model, diffMethod=diffMethod)
+sol = blockVector(assembler.numDofs)
+gridVars.init(sol)
+assembler.updateGridVariables(sol)
+print("numdofs", assembler.numDofs)
+
+# Assemble the Jacobian and the residual
+assembler.assembleJacobianAndResidual(sol)
+print("Assembly done")
+res = assembler.residual
+jac = assembler.jacobian
+
+# Solve the linear system
+S = CGSolver(jac.asLinearOperator(), SeqSSOR(jac), 1e-10)
+res *= -1
+_, _, converged, _, _ = S(sol, res)
+
+if not converged:
+    raise Exception("CGSolver has not converged")
+
+# Write to vtk
+testName = "test_1p_" + discMethod + "_" + diffMethod
+output = VtkOutputModule(gridVariables=gridVars, solutionVector=sol, name=testName)
+velocityoutput = PorousMediumFlowVelocityOutput(gridVariables=gridVars)
+output.addVelocityOutput(velocityoutput)
+output.addVolumeVariable(lambda vv: vv.pressure(), "p")
+output.write(1.0)
+
+# Print used and unused parameters
+printParameters()
diff --git a/test/python/test_explicit_transport_cctpfa.py b/test/python/test_explicit_transport_cctpfa.py
index 06b18e56055e80aff1db347106911369b4e956ff..4ad605c3a95cd79476cb861f7da1c3289be88567 100755
--- a/test/python/test_explicit_transport_cctpfa.py
+++ b/test/python/test_explicit_transport_cctpfa.py
@@ -9,6 +9,7 @@ import numpy as np
 plotting = True
 try:
     import matplotlib.pyplot as plt
+
     plt.ion()
 except ImportError:
     print("Warning: Plots are not generated as matplotlib could not be found.")
@@ -21,10 +22,9 @@ except ImportError:
 dimension = 2
 cells = 20
 
-gridView = structuredGrid([0]*dimension, [1]*dimension, [cells]*dimension)
+gridView = structuredGrid([0] * dimension, [1] * dimension, [cells] * dimension)
 
 gridGeometry = GridGeometry(gridView, discMethod="cctpfa")
-gridGeometry.update()
 
 elementMapper = gridView.indexSet
 
@@ -33,10 +33,13 @@ elementMapper = gridView.indexSet
 # Define problem (inital/boundary condtions) #
 ##############################################
 
+
 @FVProblem(gridGeometry)
 class Problem:
     numEq = 1
-    name = "finitevolume"
+
+    def name(self):
+        return "finitevolume"
 
     def boundaryTypes(self, element, scv):
         bTypes = BoundaryTypes(self.numEq)
@@ -44,7 +47,7 @@ class Problem:
         return bTypes
 
     def dirichlet(self, element, scvf):
-        if scvf.ipGlobal().two_norm < 0.25:
+        if scvf.ipGlobal.two_norm < 0.25:
             return 1.0
         else:
             return 0.0
@@ -59,7 +62,7 @@ problem = Problem()
 # Transport equation #
 ######################
 
-velocity = FieldVector([1]*dimension)
+velocity = FieldVector([1] * dimension)
 upwindWeight = 1.0
 
 
@@ -69,21 +72,23 @@ def advectiveFlux(insideConcentration, outsideConcentration, normal):
     downwindConcentration = outsideConcentration
     if normalVelocity < 0.0:
         upwindConcentration, downwindConcentration = downwindConcentration, upwindConcentration
-    return normalVelocity*(upwindWeight*upwindConcentration
-                           + (1.0-upwindWeight)*downwindConcentration)
+    return normalVelocity * (
+        upwindWeight * upwindConcentration + (1.0 - upwindWeight) * downwindConcentration
+    )
 
 
 ##########################
 # Define solution vector #
 ##########################
 
-solution = np.zeros(gridGeometry.numDofs())
+solution = np.zeros(gridGeometry.numDofs)
 
 
 ###################
 # Enable plotting #
 ###################
 
+
 @gridFunction(gridView)
 def solutionGridFunction(element, x):
     elementIdx = elementMapper.index(element)
@@ -94,13 +99,15 @@ def plot(time):
     if plotting and dimension == 2:
         fig = plt.figure()
         ax = fig.add_subplot(1, 1, 1)
-        solutionGridFunction.plot(figure=fig, clim=[0, 1+1e-6], gridLines=None)
-        ax.set_title("t = "+"{:0.2f}".format(time))
+        solutionGridFunction.plot(figure=fig, clim=[0, 1 + 1e-6], gridLines=None)
+        ax.set_title("t = " + "{:0.2f}".format(time))
         plt.show()
         plt.pause(1e-3)
-    gridView.writeVTK(problem.name + "-solution-{:0.2f}".format(time).replace(".", ""),
-                      celldata={"solution": solutionGridFunction},
-                      outputType=OutputType.ascii)
+    gridView.writeVTK(
+        problem.name + "-solution-{:0.2f}".format(time).replace(".", ""),
+        celldata={"solution": solutionGridFunction},
+        outputType=OutputType.ascii,
+    )
 
 
 #######################
@@ -117,36 +124,38 @@ plot(time=0)
 # Implement update #
 ####################
 
+
 def assembleUpdate():
-    update = np.zeros(gridGeometry.numDofs())
+    update = np.zeros(gridGeometry.numDofs)
 
     for element in gridView.elements:
-        fvGeometry = gridGeometry.localView()
+        fvGeometry = gridGeometry.localView
         fvGeometry.bind(element)
 
-        for scvf in fvGeometry.scvfs():
-            insideScvIdx = scvf.insideScvIdx()
+        for scvf in fvGeometry.scvfs:
+            insideScvIdx = scvf.insideScvIdx
             insideConcentration = solution[insideScvIdx]
 
-            if scvf.boundary():
+            if scvf.boundary:
                 bndType = problem.boundaryTypes(element, scvf)
-                if bndType.isDirichlet():
+                if bndType.isDirichlet:
                     outsideConcentration = problem.dirichlet(element, scvf)[0]
                 else:
-                    raise Exception('Only Dirichlet BCs are implemented!')
+                    raise Exception("Only Dirichlet BCs are implemented!")
 
             else:
-                outsideScvIdx = scvf.outsideScvIdx()
+                outsideScvIdx = scvf.outsideScvIdx
                 outsideConcentration = solution[outsideScvIdx]
 
-            flux = advectiveFlux(insideConcentration,
-                                 outsideConcentration,
-                                 scvf.unitOuterNormal())*scvf.area()
+            flux = (
+                advectiveFlux(insideConcentration, outsideConcentration, scvf.unitOuterNormal)
+                * scvf.area
+            )
 
             update[insideScvIdx] -= flux
 
-        for scv in fvGeometry.scvs():
-            update[scv.dofIndex()] /= scv.volume()
+        for scv in fvGeometry.scvs:
+            update[scv.dofIndex] /= scv.volume
 
     return update
 
@@ -160,13 +169,13 @@ timeLoop = TimeLoop(startTime=0.0, dt=dt, endTime=1.0, verbose=True)
 timeLoop.setPeriodicCheckPoint(0.25)
 
 timeLoop.start()
-while not timeLoop.finished():
+while not timeLoop.finished:
     update = assembleUpdate()
-    solution += timeLoop.timeStepSize() * update
+    solution += timeLoop.timeStepSize * update
     timeLoop.advanceTimeStep()
     timeLoop.reportTimeStep()
 
-    if timeLoop.isCheckPoint():
-        plot(time=timeLoop.time())
+    if timeLoop.isCheckPoint:
+        plot(time=timeLoop.time)
 
 timeLoop.finalize()
diff --git a/test/python/test_fvproblem.py b/test/python/test_fvproblem.py
old mode 100644
new mode 100755
index 74ecbf10570cf3fd7fae73c9178b5abc47e5d448..12adceaad7533c68e212c4251a6e669dc2cc569f
--- a/test/python/test_fvproblem.py
+++ b/test/python/test_fvproblem.py
@@ -14,6 +14,7 @@ def PrintProblemTest(problem):
     module = generator.load(includes, typeName, moduleName)
     return module.PrintProblemTest(problem)
 
+
 ############################################################
 # The actual Python test
 ############################################################
@@ -21,15 +22,17 @@ 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])
+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 name(self):
+        return "python_problem"
 
     def boundaryTypes(self, element, scv):
         bTypes = BoundaryTypes(self.numEq)
@@ -37,7 +40,7 @@ class Problem:
         return bTypes
 
     def dirichlet(self, element, scv):
-        if scv.center()[0] > 0.5:
+        if scv.center[0] > 0.5:
             return [0.5, 0.5]
         else:
             return [1.0, 0.0]
@@ -45,6 +48,7 @@ class Problem:
     def sourceAtPos(self, globalPos):
         return [globalPos[0]]
 
+
 problem = Problem()
 print("Name of the problem: {}".format(problem.name))
 print("-- Number of equations: {}".format(problem.numEq))
@@ -61,15 +65,14 @@ numNeumann = 0
 numDirichlet = 0
 totalSource = 0
 for e in gridView.elements:
-    fvGeometry = problem.gridGeometry().localView()
-    fvGeometry.bind(e)
-    for scv in fvGeometry.scvs():
+    fvGeometry = problem.gridGeometry.boundLocalView(e)  # test problem interface
+    for scv in fvGeometry.scvs:
         bTypes = problem.boundaryTypes(element=e, scv=scv)
-        if bTypes.isDirichlet():
+        if bTypes.isDirichlet:
             numDirichlet += 1
-        elif bTypes.isNeumann():
+        elif bTypes.isNeumann:
             numNeumann += 1
-        totalSource += problem.sourceAtPos(scv.center())[0]*scv.volume()
+        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))
+print(f"[python] Found {numNeumann} Neumann faces and {numDirichlet} Dirichlet faces")
+print(f"[python] Total source {totalSource:.2f} kg/s")
diff --git a/test/python/test_gridgeometry.py b/test/python/test_gridgeometry.py
index ca25ba36a1b50bd4d463f7f44be0eb15cfaae392..03e0a4249645411f619356865c2628e02dd6293a 100755
--- a/test/python/test_gridgeometry.py
+++ b/test/python/test_gridgeometry.py
@@ -3,19 +3,16 @@
 from dune.grid import structuredGrid
 from dumux.discretization import GridGeometry
 
-gridView = structuredGrid([0,0],[1,1],[5,5])
+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()))
+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()))
+    fvGeometry = gridGeometry.boundLocalView(e)
+    for scv in fvGeometry.scvs:
+        print(f"scv dofIndex: {scv.dofIndex}")
+        print(f"scv center: {scv.center}")
+        print(f"scv volume: {scv.volume}")