diff --git a/partitioned-heat-conduction-direct/README.md b/partitioned-heat-conduction-direct/README.md
new file mode 100644
index 000000000..4a66e09e1
--- /dev/null
+++ b/partitioned-heat-conduction-direct/README.md
@@ -0,0 +1,40 @@
+---
+title: Partitioned heat conduction (direct access setup)
+permalink: tutorials-partitioned-heat-conduction-direct.html
+keywords: Nutils, Heat conduction, Direct mesh access
+summary: This tutorial is a modified version of the "partitioned heat conduction" tutorial showcasing direct mesh access.
+---
+
+{% note %}
+Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/partitioned-heat-conduction-direct). Read how in the [tutorials introduction](https://www.precice.org/tutorials.html).
+{% endnote %}
+
+## Setup
+
+This case is a modified version of the [partitioned heat conduction tutorial](tutorials-partitioned-heat-conduction.html). Main modification is that we here use the [direct mesh access feature](couple-your-code-direct-access.html) to let the solvers compute the data mapping and not preCICE.
+
+Further minor modifications:
+
+- We use a parallel coupling scheme instead of a serial one to prevent running into the problem where we are trying to add a zero column to the quasi-Newton matrix. For serial coupling, this happens here because one data field converges much faster than the other.
+
+## Available solvers
+
+Currently only `nutils` is provided as a solver. The data mapping is computed by directly sampling the FEM function representation at the inquired locations.
+
+## Running the simulation
+
+Open two terminals and run:
+
+```bash
+cd nutils
+./run.sh -d
+```
+
+and
+
+```bash
+cd nutils
+./run.sh -n
+```
+
+See the [partitioned heat conduction tutorial](tutorials-partitioned-heat-conduction.html).
diff --git a/partitioned-heat-conduction-direct/clean-tutorial.sh b/partitioned-heat-conduction-direct/clean-tutorial.sh
new file mode 100755
index 000000000..471f8f887
--- /dev/null
+++ b/partitioned-heat-conduction-direct/clean-tutorial.sh
@@ -0,0 +1,8 @@
+#!/bin/sh
+set -e -u
+
+# shellcheck disable=SC1091
+. ../tools/cleaning-tools.sh
+
+clean_tutorial .
+
diff --git a/partitioned-heat-conduction-direct/nutils/clean.sh b/partitioned-heat-conduction-direct/nutils/clean.sh
new file mode 100755
index 000000000..6893c1ea5
--- /dev/null
+++ b/partitioned-heat-conduction-direct/nutils/clean.sh
@@ -0,0 +1,6 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+clean_nutils .
diff --git a/partitioned-heat-conduction-direct/nutils/heat.py b/partitioned-heat-conduction-direct/nutils/heat.py
new file mode 100644
index 000000000..e3ab91d98
--- /dev/null
+++ b/partitioned-heat-conduction-direct/nutils/heat.py
@@ -0,0 +1,156 @@
+#! /usr/bin/env python3
+
+from nutils import cli, mesh, function, solver, export
+import functools
+import treelog
+import numpy as np
+import precice
+
+
+def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.3):
+
+ if side == 'Dirichlet':
+ x_grid = np.linspace(0, 1, n)
+ elif side == 'Neumann':
+ x_grid = np.linspace(1, 2, n)
+ else:
+ raise Exception('invalid side {!r}'.format(side))
+ y_grid = np.linspace(0, 1, n)
+
+ # define the Nutils mesh
+ domain, geom = mesh.rectilinear([x_grid, y_grid])
+ coupling_boundary = domain.boundary['right' if side == 'Dirichlet' else 'left']
+ read_sample = coupling_boundary.sample('gauss', degree=degree * 2)
+
+ # Nutils namespace
+ ns = function.Namespace()
+ ns.x = geom
+ ns.basis = domain.basis('std', degree=degree)
+ ns.alpha = alpha # parameter of problem
+ ns.beta = beta # parameter of problem
+ ns.u = 'basis_n ?lhs_n' # solution
+ ns.dudt = 'basis_n (?lhs_n - ?lhs0_n) / ?dt' # time derivative
+ ns.flux = 'basis_n ?fluxdofs_n' # heat flux
+ ns.f = 'beta - 2 - 2 alpha' # rhs
+ ns.uexact = '1 + x_0 x_0 + alpha x_1 x_1 + beta ?t' # analytical solution
+ ns.readbasis = read_sample.basis()
+ ns.readfunc = 'readbasis_n ?readdata_n'
+
+ # define the weak form
+ res = domain.integral('(basis_n dudt - basis_n f + basis_n,i u_,i) d:x' @ ns, degree=degree * 2)
+
+ # set boundary conditions at non-coupling boundaries
+ # top and bottom boundary are non-coupling for both sides
+ sqr = domain.boundary['top,bottom,left' if side == 'Dirichlet'
+ else 'top,bottom,right'].integral('(u - uexact)^2 d:x' @ ns, degree=degree * 2)
+
+ if side == 'Dirichlet':
+ sqr += read_sample.integral('(u - readfunc)^2 d:x' @ ns)
+ else:
+ res += read_sample.integral('basis_n readfunc d:x' @ ns)
+
+ # preCICE setup
+ interface = precice.Interface(side, "../precice-config.xml", 0, 1)
+
+ mesh_id_read = interface.get_mesh_id("Dirichlet-Mesh" if side == "Dirichlet" else "Neumann-Mesh")
+ mesh_id_write = interface.get_mesh_id("Neumann-Mesh" if side == "Dirichlet" else "Dirichlet-Mesh")
+
+ vertex_ids_read = interface.set_mesh_vertices(mesh_id_read, read_sample.eval(ns.x))
+ interface.set_mesh_access_region(mesh_id_write, [.9, 1.1, -.1, 1.1])
+
+ precice_dt = interface.initialize()
+
+ vertex_ids_write, coords = interface.get_mesh_vertices_and_ids(mesh_id_write)
+ write_sample = domain.locate(ns.x, coords, eps=1e-10, tol=1e-10)
+ precice_write = functools.partial(interface.write_block_scalar_data, interface.get_data_id(
+ "Heat-Flux" if side == "Dirichlet" else "Temperature", mesh_id_write), vertex_ids_write)
+ precice_read = functools.partial(interface.read_block_scalar_data, interface.get_data_id(
+ "Temperature" if side == "Dirichlet" else "Heat-Flux", mesh_id_read), vertex_ids_read)
+
+ # helper functions to project heat flux to coupling boundary
+ if side == 'Dirichlet':
+ # To communicate the flux to the Neumann side we should not simply
+ # evaluate u_,i n_i as this is an unbounded term leading to suboptimal
+ # convergence. Instead we project ∀ v: ∫_Γ v flux = ∫_Γ v u_,i n_i and
+ # evaluate flux. While the right-hand-side contains the same unbounded
+ # term, we can use the strong identity du/dt - u_,ii = f to rewrite it
+ # to ∫_Ω [v (du/dt - f) + v_,i u_,i] - ∫_∂Ω\Γ v u_,k n_k, in which we
+ # recognize the residual and an integral over the exterior boundary.
+ # While the latter still contains the problematic unbounded term, we
+ # can use the fact that the flux is a known value at the top and bottom
+ # via the Dirichlet boundary condition, and impose it as constraints.
+ right_sqr = domain.boundary['right'].integral('flux^2 d:x' @ ns, degree=degree * 2)
+ right_cons = solver.optimize('fluxdofs', right_sqr, droptol=1e-10)
+ # right_cons is NaN in dofs that are NOT supported on the right boundary
+ flux_sqr = domain.boundary['right'].boundary['top,bottom'].integral(
+ '(flux - uexact_,0)^2 d:x' @ ns, degree=degree * 2)
+ flux_cons = solver.optimize('fluxdofs', flux_sqr, droptol=1e-10,
+ constrain=np.choose(np.isnan(right_cons), [np.nan, 0.]))
+ # flux_cons is NaN in dofs that are supported on ONLY the right boundary
+ flux_res = read_sample.integral('basis_n flux d:x' @ ns) - res
+
+ # write initial data
+ if interface.is_action_required(precice.action_write_initial_data()):
+ precice_write(write_sample.eval(0.))
+ interface.mark_action_fulfilled(precice.action_write_initial_data())
+
+ interface.initialize_data()
+
+ t = 0.
+ istep = 0
+
+ # initial condition
+ sqr0 = domain.integral('(u - uexact)^2' @ ns, degree=degree * 2)
+ lhs = solver.optimize('lhs', sqr0, arguments=dict(t=t))
+ bezier = domain.sample('bezier', degree * 2)
+
+ while interface.is_coupling_ongoing():
+
+ # save checkpoint
+ if interface.is_action_required(precice.action_write_iteration_checkpoint()):
+ checkpoint = lhs, t, istep
+ interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
+
+ # read data from interface
+ read_data = precice_read()
+
+ # prepare next timestep
+ lhs0 = lhs
+ istep += 1
+ dt = min(timestep, precice_dt)
+ t += dt
+
+ # update (time-dependent) boundary condition
+ cons = solver.optimize('lhs', sqr, droptol=1e-15, arguments=dict(t=t, readdata=read_data))
+
+ # solve nutils timestep
+ lhs = solver.solve_linear('lhs', res, constrain=cons, arguments=dict(lhs0=lhs0, dt=dt, t=t, readdata=read_data))
+
+ # write data to interface
+ if side == 'Dirichlet':
+ fluxdofs = solver.solve_linear(
+ 'fluxdofs', flux_res, arguments=dict(
+ lhs0=lhs0, lhs=lhs, dt=dt, t=t), constrain=flux_cons)
+ write_data = write_sample.eval('flux' @ ns, fluxdofs=fluxdofs)
+ else:
+ write_data = write_sample.eval('u' @ ns, lhs=lhs)
+ precice_write(write_data)
+
+ # do the coupling
+ precice_dt = interface.advance(dt)
+
+ # read checkpoint if required
+ if interface.is_action_required(precice.action_read_iteration_checkpoint()):
+ lhs, t, istep = checkpoint
+ interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
+ else:
+ # generate output
+ x, u, uexact = bezier.eval(['x_i', 'u', 'uexact'] @ ns, lhs=lhs, t=t)
+ with treelog.add(treelog.DataLog()):
+ export.vtk(side + "-" + str(istep), bezier.tri, x, Temperature=u, reference=uexact)
+
+ interface.finalize()
+
+
+if __name__ == '__main__':
+ cli.run(main)
diff --git a/partitioned-heat-conduction-direct/nutils/run.sh b/partitioned-heat-conduction-direct/nutils/run.sh
new file mode 100755
index 000000000..4b531fe9c
--- /dev/null
+++ b/partitioned-heat-conduction-direct/nutils/run.sh
@@ -0,0 +1,25 @@
+#!/bin/sh
+set -e -u
+
+usage() { echo "Usage: cmd [-d] [-n]" 1>&2; exit 1; }
+
+# Check if no input argument was provided
+if [ -z "$*" ] ; then
+ usage
+fi
+
+while getopts ":dn" opt; do
+ case ${opt} in
+ d)
+ rm -rf Dirichlet-*.vtk
+ NUTILS_RICHOUTPUT=no python3 heat.py --side=Dirichlet
+ ;;
+ n)
+ rm -rf Neumann-*.vtk
+ NUTILS_RICHOUTPUT=no python3 heat.py --side=Neumann
+ ;;
+ *)
+ usage
+ ;;
+ esac
+done
diff --git a/partitioned-heat-conduction-direct/precice-config.xml b/partitioned-heat-conduction-direct/precice-config.xml
new file mode 100644
index 000000000..b3287af65
--- /dev/null
+++ b/partitioned-heat-conduction-direct/precice-config.xml
@@ -0,0 +1,60 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<precice-configuration>
+ <log>
+ <sink
+ filter="%Severity% > debug and %Rank% = 0"
+ format="---[precice] %ColorizedSeverity% %Message%"
+ enabled="true" />
+ </log>
+
+ <solver-interface dimensions="2" experimental="yes">
+ <data:scalar name="Temperature" />
+ <data:scalar name="Heat-Flux" />
+
+ <mesh name="Dirichlet-Mesh">
+ <use-data name="Temperature" />
+ <use-data name="Heat-Flux" />
+ </mesh>
+
+ <mesh name="Neumann-Mesh">
+ <use-data name="Temperature" />
+ <use-data name="Heat-Flux" />
+ </mesh>
+
+ <participant name="Dirichlet">
+ <use-mesh name="Dirichlet-Mesh" provide="yes" />
+ <use-mesh name="Neumann-Mesh" from="Neumann" direct-access="true" />
+ <write-data name="Heat-Flux" mesh="Neumann-Mesh" />
+ <read-data name="Temperature" mesh="Dirichlet-Mesh" />
+ </participant>
+
+ <participant name="Neumann">
+ <use-mesh name="Neumann-Mesh" provide="yes" />
+ <use-mesh name="Dirichlet-Mesh" from="Dirichlet" direct-access="true" />
+ <write-data name="Temperature" mesh="Dirichlet-Mesh" />
+ <read-data name="Heat-Flux" mesh="Neumann-Mesh" />
+ </participant>
+
+ <m2n:sockets from="Dirichlet" to="Neumann" exchange-directory=".." />
+
+ <coupling-scheme:parallel-implicit>
+ <participants first="Dirichlet" second="Neumann" />
+ <max-time value="1.0" />
+ <time-window-size value="0.1" />
+ <max-iterations value="100" />
+ <exchange data="Heat-Flux" mesh="Neumann-Mesh" from="Dirichlet" to="Neumann" />
+ <exchange data="Temperature" mesh="Dirichlet-Mesh" from="Neumann" to="Dirichlet" />
+ <relative-convergence-measure data="Heat-Flux" mesh="Neumann-Mesh" limit="1e-5" />
+ <relative-convergence-measure data="Temperature" mesh="Dirichlet-Mesh" limit="1e-5" />
+ <acceleration:IQN-ILS>
+ <data name="Temperature" mesh="Dirichlet-Mesh" />
+ <data name="Heat-Flux" mesh="Neumann-Mesh" />
+ <initial-relaxation value="0.1" />
+ <max-used-iterations value="10" />
+ <time-windows-reused value="5" />
+ <preconditioner type="residual-sum" />
+ <filter type="QR2" limit="1e-3" />
+ </acceleration:IQN-ILS>
+ </coupling-scheme:parallel-implicit>
+ </solver-interface>
+</precice-configuration>
diff --git a/partitioned-heat-conduction/nutils/run.sh b/partitioned-heat-conduction/nutils/run.sh
old mode 100644
new mode 100755
index 60e3f3911..4b531fe9c
--- a/partitioned-heat-conduction/nutils/run.sh
+++ b/partitioned-heat-conduction/nutils/run.sh
@@ -12,11 +12,11 @@ while getopts ":dn" opt; do
case ${opt} in
d)
rm -rf Dirichlet-*.vtk
- NUTILS_RICHOUTPUT=no python3 heat.py side=Dirichlet
+ NUTILS_RICHOUTPUT=no python3 heat.py --side=Dirichlet
;;
n)
rm -rf Neumann-*.vtk
- NUTILS_RICHOUTPUT=no python3 heat.py side=Neumann
+ NUTILS_RICHOUTPUT=no python3 heat.py --side=Neumann
;;
*)
usage