From 51023a23a91e24c231e007d66f0d4288dc60c389 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 17 Jul 2020 17:05:48 +0200
Subject: [PATCH] setup new TPR layered soil examples

---
 Plotskripte_Jim/links                         |    1 +
 Plotskripte_Jim/plot2d_vtu.py                 |  502 ++++++++
 Plotskripte_Jim/plot3d_vtu.py                 |  510 +++++++++
 ...patch-realistic-different-permeabilties.py |   50 +-
 ...soil_with_inner_patch-realistic-pure-dd.py | 1015 +++++++++++++++++
 ...layered_soil_with_inner_patch-realistic.py |    1 -
 6 files changed, 2053 insertions(+), 26 deletions(-)
 create mode 100644 Plotskripte_Jim/links
 create mode 100644 Plotskripte_Jim/plot2d_vtu.py
 create mode 100644 Plotskripte_Jim/plot3d_vtu.py
 create mode 100755 Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic-pure-dd.py

diff --git a/Plotskripte_Jim/links b/Plotskripte_Jim/links
new file mode 100644
index 0000000..bffbaf4
--- /dev/null
+++ b/Plotskripte_Jim/links
@@ -0,0 +1 @@
+https://lorensen.github.io/VTKExamples/site/Python/
diff --git a/Plotskripte_Jim/plot2d_vtu.py b/Plotskripte_Jim/plot2d_vtu.py
new file mode 100644
index 0000000..08d5b09
--- /dev/null
+++ b/Plotskripte_Jim/plot2d_vtu.py
@@ -0,0 +1,502 @@
+#!/usr/bin/python
+
+# Python script to plot .vtk and .vtu data
+# and adding streamlines.
+
+import os
+from os.path import basename
+import sys
+from optparse import OptionParser
+import subprocess
+
+import re
+
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import matplotlib.tri as mtri
+from matplotlib import cm
+import matplotlib.patches
+
+from matplotlib.backends.backend_pdf import PdfPages
+from matplotlib.colors import LinearSegmentedColormap
+
+import fileinput
+from glob import glob
+
+from scipy.interpolate import griddata
+import numpy as np
+from vtk import *
+from vtk.util.numpy_support import vtk_to_numpy
+from vtk import vtkUnstructuredGridReader
+from vtk import vtkXMLUnstructuredGridReader
+
+
+def cm2inch(*tupl):
+    # Convert centimeters to inch
+    inch = 2.54
+    if isinstance(tupl[0], tuple):
+        return tuple(i / inch for i in tupl[0])
+    else:
+        return tuple(i / inch for i in tupl)
+
+
+quiver_pos_x = None
+quiver_pos_y = None
+
+#========================================================================
+
+def plot_vtk_file(filename, interface_filename=None, scalar_variable=None, datarange = None, start_points = None, cmap_binedges = None, plot_velocity=False, output_folder='default', cmap=plt.cm.Spectral_r):
+
+    quivercolor = 'white'
+    interfacecolor = 'black'
+    # cmap = plt.cm.Spectral_r
+    # # cmap = plt.cm.cividis
+    # cmap = plt.cm.RdBu_r
+    # # cmap = plt.cm.Blues
+
+    # Get the filetype
+    extension = os.path.splitext(filename)[1]
+
+    if (extension == '.vtu'):
+        reader = vtkXMLUnstructuredGridReader()
+    elif extension == '.vtk':
+        reader = vtkUnstructuredGridReader()
+
+    reader.SetFileName(filename)
+    # reader.ReadAllVectorsOn()
+    # reader.ReadAllScalarsOn()
+    reader.Update()
+
+
+    if interface_filename is not None:
+        # Get the filetype
+        interface_extension = os.path.splitext(interface_filename)[1]
+        if (interface_extension == '.vtu'):
+            interface_reader = vtkXMLUnstructuredGridReader()
+        elif interface_extension == '.vtk':
+            interface_reader = vtkUnstructuredGridReader()
+
+        interface_reader.SetFileName(interface_filename)
+        interface_reader.Update()
+
+    # Get the coordinates of nodes in the mesh
+    nodes_vtk_array = reader.GetOutput().GetPoints().GetData()
+
+    print('--- Read cell values...')
+    n_cell_arrays = reader.GetOutput().GetCellData().GetNumberOfArrays()
+    n_point_arrays = reader.GetOutput().GetPointData().GetNumberOfArrays()
+
+    scalar_idx = -1
+    for i in range(n_cell_arrays):
+        var_name = reader.GetOutput().GetCellData().GetArrayName(i)
+        if (var_name == scalar_variable):
+            scalar_idx = i
+    # print("scalar_idx = {}".format(scalar_idx))
+    if scalar_idx >= 0:
+        scalar_vtk_array = reader.GetOutput().GetCellData().GetArray(scalar_idx)
+
+    # time_field = reader.GetOutput().GetFieldData().GetArray(0)
+    # simulation_time = vtk_to_numpy(time_field)[0]
+
+    #####################################
+    nodes_numpy_array = vtk_to_numpy(nodes_vtk_array)
+    x, y, z = nodes_numpy_array[:, 0], nodes_numpy_array[
+        :, 1], nodes_numpy_array[:, 2]
+
+    if scalar_idx >= 0:
+        scalar_numpy_array = vtk_to_numpy(scalar_vtk_array)
+        if (extension == '.vtk'):
+            scalar_numpy_array = scalar_numpy_array[:,scalar_idx]
+
+    try:
+        if (np.shape(scalar_numpy_array)[1] > 1):
+            scalar_numpy_array = np.linalg.norm(scalar_numpy_array, axis=1)
+    except Exception as e:
+        print(e)
+        # raise
+
+    tri_vtk_array = reader.GetOutput().GetCells().GetData()
+    tri = vtk_to_numpy(tri_vtk_array)
+    triangles = tri
+
+    ntri = np.shape(triangles)
+    ntri = (ntri[0]) / (3 + 1)
+    ntri = int(ntri)
+
+    ia = np.zeros(ntri, dtype=int)
+    ib = np.zeros(ntri, dtype=int)
+    ic = np.zeros(ntri, dtype=int)
+
+    cell_centers_x = np.zeros(ntri)
+    cell_centers_y = np.zeros(ntri)
+
+    # triangles are formatted like (nPoints, id1, id2, id3, ...., idnPoints)
+    for i in range(0, ntri):
+        ia[i] = int(triangles[4 * i + 1])
+        ib[i] = int(triangles[4 * i + 2])
+        ic[i] = int(triangles[4 * i + 3])
+        cell_centers_x[i] = (x[ia[i]] + x[ib[i]] + x[ic[i]]) / 3.0
+        cell_centers_y[i] = (y[ia[i]] + y[ib[i]] + y[ic[i]]) / 3.0
+
+    triangles = np.vstack((ia, ib, ic))
+    triangles = triangles.T
+    #triangles = np.reshape(triangles, (-1, M))
+    # print(triangles)
+    # print(np.shape(triangles))
+
+    ############################################
+
+    triangulation = mtri.Triangulation(x=x, y=y, triangles=triangles)
+
+    # ############################################
+    #
+    # print('--- Map to grid...')
+    #
+    # # make velocity mesh grid:
+    # if(plot_velocity):
+    #     resolution = 200
+    # else:
+    #     resolution = 100
+    # xi, yi = np.meshgrid(np.linspace(x.min(), x.max(), resolution),
+    #                      np.linspace(y.min(), y.max(), resolution))
+    #
+    # # 'cubic', 'nearest', 'linear'
+    # if (plot_velocity):
+    #     ui = griddata((cell_centers_x, cell_centers_y), u, (xi, yi), method='cubic')
+    #     vi = griddata((cell_centers_x, cell_centers_y), v, (xi, yi), method='cubic')
+    #
+    # print('--- Mask values...')
+    # print('--- (Note: Please adapt this manually to your application)')
+    # # Mask Values:
+    # mask = grid_mask(xi, yi)
+    #
+    # if (plot_velocity):
+    #     vi[mask] = np.nan
+    #     ui[mask] = np.nan
+    #
+    #     vel_magnitude = np.sqrt(ui * ui + vi * vi)
+    #
+    # ############################################
+
+    print('--- Plot...')
+    matplotlib.style.use('ggplot')
+    # matplotlib.style.use('default')
+    #matplotlib.style.use(['seaborn-deep', 'seaborn-talk'])
+    # matplotlib.style.use('seaborn-darkgrid')
+    # matplotlib.style.use('seaborn-whitegrid')
+
+    matplotlib.rc('font', family='sans-serif')
+    matplotlib.rc('font', serif='Biolinum')
+    matplotlib.rc('text', usetex='false')
+
+    # plt.figure(figsize=(12,12))
+    #plt.figure(figsize=cm2inch(30, 30))
+    # plt.figure(figsize=cm2inch(28, 28))
+    plt.figure(figsize=cm2inch(15, 15))
+
+    ax = plt.gca()
+    ax.set_xlim(x.min(), x.max())
+    ax.set_ylim(y.min(), y.max())
+    # plt.gca().set_aspect('equal')
+    ax.set_aspect('equal')
+
+    # cmap_lin = plt.cm.Spectral
+    # cmap = plt.cm.Spectral_r
+    # cmap_lin = plt.cm.RdBu
+    # cmap_lin = plt.get_cmap('plasma')
+    # cmap_lin = plt.get_cmap('viridis')
+    # cmap_lin = matplotlib.colors.ListedColormap(cmap_lin.colors[::-1])
+    if cmap_binedges is not None:
+        norm = matplotlib.colors.BoundaryNorm(boundaries=cmap_binedges, ncolors=256)
+    else:
+        norm = None
+
+
+    #plt.tripcolor(x, y, triangles, facecolors=scalar_numpy_array, edgecolors='k', cmap=plt.cm.Spectral)
+    # im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=plt.cm.Spectral_r)
+    # im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=plt.cm.Spectral)
+    if scalar_idx >= 0:
+        im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=cmap, norm=norm)
+    # im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=plt.get_cmap('viridis'))
+    #im = plt.tripcolor(triangulation, facecolors=scalar_numpy_array, cmap=plt.cm.viridis)
+    #plt.tripcolor(x, y, triangles, facecolors=scalar_numpy_array, cmap=plt.cm.Spectral, shading='gouraud')
+
+    #im = plt.triplot(x, y, triangles, color='grey', alpha=0.2, linewidth=0.25)
+    ax.triplot(triangulation, color='black', alpha=0.15, linewidth=0.2)
+    # plt.title('t = {:.2e}'.format(simulation_time))
+    plt.xlabel('x')
+    plt.ylabel('y')
+
+    if datarange is not None:
+        if (len(datarange) == 2):
+            plt.clim(datarange[0], datarange[1])
+
+
+    if interface_filename is not None:
+        nodes_vtk_array = interface_reader.GetOutput().GetPoints().GetData()
+        nodes_numpy_array = vtk_to_numpy(nodes_vtk_array)
+        xin, yin, zin = nodes_numpy_array[:, 0], nodes_numpy_array[:, 1], nodes_numpy_array[:, 2]
+
+        int_vtk_array = interface_reader.GetOutput().GetCells().GetData()
+        interface_array = vtk_to_numpy(int_vtk_array)
+
+        ninterface = np.shape(interface_array)
+        ninterface = (ninterface[0]) / (2 + 1)
+        ninterface = int(ninterface)
+
+        # cells are formatted like (nPoints, id1, id2, id3, ...., idnPoints)
+        for i in range(0, ninterface):
+            ia = int(interface_array[3 * i + 1])
+            ib = int(interface_array[3 * i + 2])
+            ax.plot([xin[ia], xin[ib]], [yin[ia], yin[ib]], alpha=1, lw=0.5, color=interfacecolor, solid_capstyle='round')
+            # interface_centers_x[i] = (xin[ia[i]] + xin[ib[i]]) / 2.0
+            # interface_centers_y[i] = (yin[ia[i]] + yin[ib[i]]) / 2.0
+
+    # create an axes on the right side of ax. The width of cax will be 5%
+    # of ax and the padding between cax and ax will be fixed at 0.1 inch.
+    if scalar_idx >= 0:
+        divider = make_axes_locatable(ax)
+        cax = divider.append_axes("right", size="3%", pad=0.3)
+
+        cbar = plt.colorbar(im, cax=cax, spacing='proportional')
+        # cbar = plt.colorbar(im, cax=cax)
+
+
+        # Name the colorbar automatically:
+        scalar_label = reader.GetOutput().GetCellData().GetArrayName(scalar_idx)
+        cbar.set_label(scalar_label, labelpad=5)
+        # cbar.set_label(scalar_variable, labelpad=5)
+
+    if (plot_velocity):
+        # plot_streams = True
+        plot_streams = False
+        # plot_quiver = False
+        plot_quiver = True
+        plot_contour = False
+
+        for i in range(n_cell_arrays):
+            var_name = reader.GetOutput().GetCellData().GetArrayName(i)
+            if (var_name == 'velocity'):
+                vector_idx = i
+
+        # print("vector_idx = {}".format(vector_idx))
+        vector_vtk_array = reader.GetOutput().GetCellData().GetArray(vector_idx)
+        v_numpy_array = vtk_to_numpy(vector_vtk_array)
+        u = v_numpy_array[:,0]
+        v = v_numpy_array[:,1]
+
+        if(plot_streams):
+            print('--- add streamlines...')
+
+            #ax.streamplot(xi, yi, ui, vi, color=vel_magnitude, linewidth=1, cmap=plt.cm.YlGnBl_r)
+            lw = vel_magnitude
+            lw[np.abs(lw) < 1e-9] = 0
+            lw = lw / np.nanmax(lw)
+            if(np.nanmax(lw) != np.nanmin(lw)):
+                lw = 2 * lw
+
+
+            # cont = ax.contour(xi, yi, vel_magnitude, 3, cmap=plt.get_cmap('plasma'), linewidths=0.5, alpha=0.33)
+            if (start_points is None):
+                stream_density = 0.66
+                # stream_density = 0.5
+                stream_objects = ax.streamplot(
+                xi, yi, ui, vi, color=vel_magnitude, linewidth=lw, arrowsize=3.0, cmap=plt.get_cmap('plasma'), density=stream_density)
+                    # xi, yi, ui, vi, color=vel_magnitude, linewidth=1.2, arrowsize=1.33, cmap=plt.get_cmap('plasma'), density=0.75)
+            else:
+                xxi = np.linspace(x.min(), x.max(), resolution)
+                yyi = np.linspace(y.min(), y.max(), resolution)
+                stream_objects = ax.streamplot(
+                xxi, yyi, ui, vi, color=vel_magnitude, linewidth=lw, arrowsize=3.0, cmap=plt.get_cmap('plasma'), start_points=start_points.T, density=200)
+                    # xi, yi, ui, vi, color=vel_magnitude, linewidth=1.2, arrowsize=1.33, cmap=plt.get_cmap('plasma'), density=0.75)
+
+            cax_streams = divider.append_axes("right", size="3%", pad=0.7)
+            cbar_streams = plt.colorbar(stream_objects.lines, cax=cax_streams)
+
+            #cbar_streams.ax.set_ylabel('Velocity', rotation=270)
+            cbar_streams.set_label('velocity', labelpad=5)
+
+            stream_alpha = 0.8
+            # stream_alpha = 0.5
+            stream_objects.lines.set_alpha(stream_alpha)
+            stream_objects.arrows.set_alpha(stream_alpha)
+
+            # for x in plt.gca().get_children():
+            # if type(x) == matplotlib.patches.FancyArrowPatch:
+            #     x.set_alpha(stream_alpha)  # or x.set_visible(False)
+            #     start_points = [start_points, 0.5 * (x.posA + x.posB)]
+            #     print(start_points)
+
+        if(plot_contour):
+            n_contourlines = 3
+            cont = ax.contour(xi, yi, vel_magnitude, n_contourlines, cmap=plt.get_cmap('plasma'), linewidths=0.5, alpha=0.33)
+            cax_streams = divider.append_axes("right", size="3%", pad=0.7)
+            cbar_streams = plt.colorbar(cont, cax=cax_streams)
+            #cbar_streams.ax.set_ylabel('Velocity', rotation=270)
+            cbar_streams.set_label('velocity', labelpad=5)
+
+
+        if(plot_quiver):
+            print('Plot quiver...')
+
+            global quiver_pos_x
+            global quiver_pos_y
+            if quiver_pos_x is None:
+                # skip=slice(None,None,25)
+                skip=slice(None,None,75)
+                quiver_pos_x = cell_centers_x[skip]
+                quiver_pos_y = cell_centers_y[skip]
+
+            finder = triangulation.get_trifinder()
+            tri_indices = finder(quiver_pos_x, quiver_pos_y)
+            mu = u[tri_indices].copy()
+            mv = v[tri_indices].copy()
+            mask = np.sqrt(mu**2 + mv**2) < 1e-12
+            mx = np.ma.masked_where(mask, quiver_pos_x)
+            my = np.ma.masked_where(mask, quiver_pos_y)
+            mu = np.ma.masked_where(mask, mu)
+            mv = np.ma.masked_where(mask, mv)
+
+            # # colormapped vectors:
+            # col = np.sqrt(u**2 + v**2)
+            # Q = ax.quiver(cell_centers_x[skip], cell_centers_y[skip], u[skip], v[skip], col[skip], units='width', cmap=plt.get_cmap('plasma'), alpha=0.8)
+            # uniform color:
+            Q = ax.quiver(mx, my, mu, mv,  units='width', color=quivercolor, alpha=0.66)
+
+            # qk = plt.quiverkey(Q, 0.9, 0.95, 2, r'$2 \frac{m}{s}$',
+            #        labelpos='E',
+            #        coordinates='figure',
+            #        fontproperties={'weight': 'bold'})
+
+            # cax_streams = divider.append_axes("right", size="3%", pad=0.7)
+            # cbar_streams = plt.colorbar(Q, cax=cax_streams)
+            # #cbar_streams.ax.set_ylabel('Velocity', rotation=270)
+            # cbar_streams.set_label('velocity', labelpad=5)
+
+    filename_base = basename(filename)
+    filename_base = os.path.splitext(filename_base)[0]
+    #filename_base = os.path.splitext(filename)[0]
+
+    if not os.path.exists('plots'):
+        print('--- Make directory...')
+        os.makedirs('plots', exist_ok=True)
+    output_dir = 'plots/{}'.format(output_folder)
+    if not os.path.exists(output_dir):
+        print('--- Make directory...')
+        os.makedirs(output_dir, exist_ok=True)
+    # if not os.path.exists('plots/png'):
+    #     print('--- Make directory...')
+    #     os.makedirs('plots/png', exist_ok=True)
+    # if not os.path.exists('plots/pdf'):
+    #     print('--- Make directory...')
+    #     os.makedirs('plots/pdf', exist_ok=True)
+
+
+    print('--- Save plot...')
+    plt.tight_layout()
+
+    # get index from filename
+    # plot_index = map(int, re.findall('\d+', filename_base))
+    # plot_index = plot_index[0]
+    # plt.savefig('plots/jpg/frame_%09d.jpg' % plot_index,
+    #             bbox_inches='tight', dpi=200)  # dpi=300 200
+    plt.savefig(os.path.join(output_dir, '{}.jpg'.format(filename_base)),
+                bbox_inches='tight', dpi=300)  # dpi=300 200
+
+    # regex = re.compile(r'\d+')
+    # file_number = regex.search(filename).group(0)
+    # # regex.findall(filename_base)
+    # plt.savefig(os.path.join(output_dir, '{}_{}.jpg'.format(scalar_variable, file_number)),
+    #             bbox_inches='tight', dpi=300)  # dpi=300 200
+
+    #print('Save pdf...')
+    #plt.savefig('pdfs/'+filename_base+'.pdf', bbox_inches='tight', dpi=300)
+    # plt.show()
+
+    # pp = PdfPages('plots/pdf/' + filename_base + '.pdf')
+    # pp.savefig()
+    # pp.close()
+
+    plt.close()
+
+    print("--- Done!")
+    return 1
+
+############################################
+
+
+if __name__ == "__main__":
+    fn = './data/global_mesh000000.vtu'
+    # fn = './data/domain_marker000000.vtu'
+    output_folder = 'marker'
+    plot_vtk_file(fn, scalar_variable='f', output_folder=output_folder)
+
+
+    # parser = OptionParser()
+    # parser.add_option("-f", "--file", dest="filename",
+    #                   default=False, help="open FILE", metavar="FILE")
+    # parser.add_option("-q", "--quiet",
+    #                   action="store_false", dest="verbose", default=False,
+    #                   help="don't print status messages to stdout")
+    #
+    # (options, args) = parser.parse_args()
+    #
+    # # If no option is given
+    # if (options.filename == False):
+    #     # If no option is given, plot the whole
+    #     # content of the input-directory:
+    #     input_directory = '../build/bin/out'
+    #     # input_directory = '../build_bc_fix/bin/out'
+    #     # input_directory = '../build_2/bin/out'
+    #     # input_directory = '../build_3/bin/out'
+    #     print("Input directory: %s" % input_directory)
+    #
+    #     filenames = sorted(glob(input_directory + '/solution*.vt*'))
+    #     interface_filenames = sorted(glob(input_directory + '/interface*.vt*'))
+    #     print("Found the following files:")
+    #     print(filenames)
+    #
+    #     cmap_binedges = None
+    #     datarange = None
+    #
+    #     # output = 'velocity'
+    #     # output = 'temperature'
+    #     output = 'density'
+    #     plot_velocity = True
+    #     plot_velocity = False
+    #     # plot_schlieren = True
+    #
+    #     cmap = plt.cm.Spectral_r
+    #     # cmap = plt.cm.cividis
+    #     # cmap = plt.cm.RdBu_r
+    #     # cmap = plt.cm.Blues
+    #
+    #     print("get data range...")
+    #     datarange = get_data_range(filenames, output)
+    #     # datarange = get_data_range(filenames, output, separat_ranges=True)
+    #     # cmap_binedges = get_data_bins(filenames, output)
+    #
+    #     print("Start plotting...")
+    #     output_interval = 50
+    #     output_interval = 1
+    #     # for fn in filenames[::output_interval]:
+    #     for fn, fni in zip(filenames[::output_interval], interface_filenames[::output_interval]):
+    #         print("Plot file {}, {}".format(fn, fni))
+    #         plot_vtk_file(fn, interface_filename=fni, scalar_variable=output, cmap_binedges=cmap_binedges, datarange=datarange, plot_velocity=plot_velocity, cmap=cmap)
+    #         # schlieren_plot_vtk_file(fn, interface_filename=fni, scalar_variable=output, cmap_binedges=cmap_binedges, datarange=datarange, plot_velocity=plot_velocity, cmap=cmap)
+    #
+    #     video_length = 8.0
+    #     nframes = len(filenames[::output_interval])
+    #     framerate = nframes / video_length
+    #     # framerate = 60
+    #     print('framerate = {}'.format(framerate))
+    #
+    #     command = ' ffmpeg -framerate {} -pattern_type glob -i "./plots/{}/solution_*.jpg" -filter:v "crop=2*floor(in_w/2):2*floor(in_h/2)" -pix_fmt yuv420p -vcodec h264 -r 60 ./plots/{}.mp4 '.format(framerate, output, output, output)
+    #     # command = ' ffmpeg -framerate {} -pattern_type glob -i "./plots/{}/*.jpg" -filter:v "crop=2*floor(in_w/2):2*floor(in_h/2)" -pix_fmt yuv420p -vcodec h264 -r 60 ./plots/{}.mp4 '.format(framerate, output, output)
+    #     subprocess.call(command, shell=True)
+    #
+    # else:
+    #     # If the file is given as a option, plot it
+    #     print("Plot file %s" % options.filename)
+    #     plot_vtk_file(options.filename)
diff --git a/Plotskripte_Jim/plot3d_vtu.py b/Plotskripte_Jim/plot3d_vtu.py
new file mode 100644
index 0000000..683a186
--- /dev/null
+++ b/Plotskripte_Jim/plot3d_vtu.py
@@ -0,0 +1,510 @@
+from glob import glob
+import os
+import sys
+# import fileinput
+# import re
+import shutil
+
+import numpy as np
+# import math
+# import random
+
+# import matplotlib.pyplot as plt
+# from mpl_toolkits.mplot3d import Axes3D
+import matplotlib
+import matplotlib.pyplot as plt
+from matplotlib.pyplot import cm
+
+# from vtk import *
+import vtk
+from vtk.util.numpy_support import vtk_to_numpy
+from vtk import vtkUnstructuredGridReader
+from vtk import vtkXMLUnstructuredGridReader
+
+os.environ["PV_ALLOW_BATCH_INTERACTION"] = "1"
+
+
+
+# ------------------------------------------------------------------
+
+def cm2inch(*tupl):
+    # Convert centimeters to inch
+    inch = 2.54
+    if isinstance(tupl[0], tuple):
+        return tuple(i / inch for i in tupl[0])
+    else:
+        return tuple(i / inch for i in tupl)
+
+# ------------------------------------------------------------------
+
+
+def plot3d_vtu(filename, interface_filename=None, plotname=None, cmap=cm.Spectral_r):
+
+    scalar_variable = 'phase'
+    # interfacecolor = '#0F75FF'
+    interfacecolor = '#0050CC'
+    phasecolor = '#73a0e6'
+
+    # cmap = sns.light_palette(interfacecolor, as_cmap=True, reverse=False, n_colors=6)
+    # cmap = cmap(np.linspace(0.0, 0.2, 10))
+    # cmap = LinearSegmentedColormap.from_list('truncated', cmap)
+    # cmap = matplotlib.colors.ListedColormap(['white', phasecolor])
+    # alpha=0.3
+    # black_cmap = matplotlib.colors.ListedColormap(['white', 'black'])
+    cmap = cm.get_cmap('Spectral_r')
+
+
+    # Get the filetype
+    extension = os.path.splitext(filename)[1]
+
+    if (extension == '.vtu'):
+        reader = vtkXMLUnstructuredGridReader()
+    elif extension == '.vtk':
+        reader = vtkUnstructuredGridReader()
+
+    reader.SetFileName(filename)
+    reader.Update()
+    if interface_filename is not None:
+        # Get the filetype
+        interface_extension = os.path.splitext(interface_filename)[1]
+        if (interface_extension == '.vtu'):
+            interface_reader = vtkXMLUnstructuredGridReader()
+        elif interface_extension == '.vtk':
+            interface_reader = vtkUnstructuredGridReader()
+
+        interface_reader.SetFileName(interface_filename)
+        interface_reader.Update()
+        interface_output = interface_reader.GetOutput()
+
+    print('--- Read cell values...')
+    n_cell_arrays = reader.GetOutput().GetCellData().GetNumberOfArrays()
+    n_point_arrays = reader.GetOutput().GetPointData().GetNumberOfArrays()
+
+    # time_field = reader.GetOutput().GetFieldData().GetArray(0)
+    # simulation_time = vtk_to_numpy(time_field)[0]
+
+    modified_idx = -1
+    for i in range(n_cell_arrays):
+        var_name = reader.GetOutput().GetCellData().GetArrayName(i)
+        print(var_name)
+        if (var_name == scalar_variable):
+            scalar_idx = i
+        if (var_name == 'modified'):
+            modified_idx = i
+
+
+    # scalar_vtk_array = reader.GetOutput().GetCellData().GetArray(scalar_idx)
+    # reader.GetOutput().GetCellData().SetActiveScalars(scalar_variable)
+
+    # reader.GetOutput().GetCellData().SetActiveScalars('phase')
+
+    # output = interface_reader.GetOutput()
+    output = reader.GetOutput()
+
+
+    bounds = reader.GetOutput().GetBounds()
+    center = reader.GetOutput().GetCenter()
+    scalar_range = output.GetScalarRange()
+    print(scalar_range)
+
+    # ----------------------------------------
+
+    # cmap_inst.set_bad(color='white')
+    vmin, vmax = scalar_range
+    lut = vtk.vtkLookupTable()
+    lut.SetTableRange(vmin, vmax)
+    # lut.SetNanColor(.1, .5, .99, 1.0)
+    n = 256
+    lut.SetNumberOfTableValues(n)
+    lut.Build()
+    # for i in range(n):
+    for i, val in enumerate(np.linspace(0, 1, n)):
+        # R,G,B = colors.colorMap(i, 'jet', 0, n)
+        rgb_val = cmap(val)
+        lut.SetTableValue(i, rgb_val[0], rgb_val[1], rgb_val[2], 1)
+
+    # ----------------------------------------
+
+    mapper = vtk.vtkDataSetMapper()
+    mapper.SetInputData(output)
+    # mapper.SetInputConnection(output_port)
+    mapper.SetScalarRange(scalar_range)
+    mapper.SetScalarModeToUseCellData()
+    # mapper.SelectColorArray('phase')
+    mapper.SetLookupTable(lut)
+
+    # ----------------------------------------
+    #
+    # interface_mapper = vtk.vtkDataSetMapper()
+    # interface_mapper.SetInputData(interface_output)
+    # interface_actor = vtk.vtkActor()
+    # interface_actor.SetMapper(interface_mapper)
+    # interface_actor.GetProperty().EdgeVisibilityOn()
+    # interface_actor.GetProperty().SetLineWidth(1.0)
+    # interface_actor.GetProperty().SetColor(matplotlib.colors.to_rgb(interfacecolor))
+    # # interface_actor.GetProperty().SetSpecular(0)
+    # # interface_actor.GetProperty().SetSpecularPower(0)
+    # interface_actor.GetProperty().SetOpacity(0.2)
+    # interface_actor.GetProperty().SetRepresentationToWireframe()
+    # # interface_actor.GetProperty().SetAmbient(0.1)
+    # # interface_actor.GetProperty().SetDiffuse(1)
+
+
+    # ----------------------------------------
+    #
+    # plane = vtk.vtkPlane()
+    # plane.SetOrigin(output.GetCenter())
+    # plane.SetNormal(0, 0, 1)
+    # plane_cut = vtk.vtkCutter()
+    # plane_cut.GenerateTrianglesOff()
+    # plane_cut.SetInputData(output)
+    # plane_cut.SetCutFunction(plane)
+    # cut_mapper = vtk.vtkDataSetMapper()
+    # cut_mapper.SetInputConnection(plane_cut.GetOutputPort())
+    # cut_mapper.SetScalarRange(scalar_range)
+    # cut_mapper.ScalarVisibilityOn()
+    # cut_mapper.SetScalarModeToUseCellData()
+    # cut_mapper.SetLookupTable(lut)
+
+    # ----------------------------------------
+    #
+    # extracter = vtk.vtkExtractGeometry()
+    #
+    # extracter.SetInputData(output);
+    # extracter.SetImplicitFunction(plane);
+    # extracter.ExtractInsideOn();
+    # extracter.ExtractBoundaryCellsOn();
+    #
+    # extract_mapper = vtk.vtkDataSetMapper()
+    # extract_mapper.SetInputConnection(extracter.GetOutputPort())
+    # extract_mapper.SetScalarRange(scalar_range)
+    # extract_mapper.ScalarVisibilityOn()
+    # extract_mapper.SetScalarModeToUseCellData()
+    # extract_mapper.SetLookupTable(lut)
+
+    # # ----------------------------------------
+    #
+    # clipper = vtk.vtkClipDataSet()
+    # clipper.SetClipFunction(plane)
+    # clipper.SetInputData(reader.GetOutput())
+    # clipper.InsideOutOn()
+    # clipper.SetValue(0.0)
+    # clipper.GenerateClippedOutputOn()
+    # clipper.Update()
+    #
+    # clip_mapper = vtk.vtkDataSetMapper()
+    # clip_mapper.SetInputData(clipper.GetOutput())
+    #
+    # # ----------------------------------------
+
+    main_mapper = mapper
+    # Create the Actor
+    actor = vtk.vtkActor()
+    actor.SetMapper(main_mapper)
+    # actor.SetMapper(clip_mapper)
+    # actor.SetMapper(extract_mapper)
+    # actor.SetMapper(cut_mapper)
+    actor.GetProperty().EdgeVisibilityOn()
+    actor.GetProperty().SetLineWidth(1.0)
+    actor.GetProperty().SetSpecular(0)
+    actor.GetProperty().SetSpecularPower(0)
+    # actor.GetProperty().SetOpacity(0.5)
+    # actor.GetProperty().SetRepresentationToWireframe()
+    actor.GetProperty().SetAmbient(1)
+    actor.GetProperty().SetDiffuse(0)
+
+    # ----------------------------------------
+
+    colorbar = vtk.vtkScalarBarActor()
+    colorbar.SetLookupTable(main_mapper.GetLookupTable())
+    colorbar.SetTitle(scalar_variable)
+    # colorbar.SetOrientationToHorizontal()
+    colorbar.SetOrientationToVertical()
+    # colorbar.GetLabelTextProperty().SetColor(0,0,0)
+    # colorbar.GetTitleTextProperty().SetColor(0,0,0)
+    tprop = vtk.vtkTextProperty()
+    # tprop.SetFontFamily(vtk.VTK_FONT_FILE)
+    # # tprop.SetFontFamilyAsString('Linux Biolinum O')
+    # # tprop.SetFontFile('./LibertinusSans-Bold.otf')
+    # tprop.SetFontFile('./LibertinusSans-Regular.otf')
+    tprop.SetColor(0,0,0)
+    # tprop.BoldOn()
+    colorbar.SetLabelTextProperty(tprop)
+    colorbar.SetTitleTextProperty(tprop)
+    # position it in window
+    coord = colorbar.GetPositionCoordinate()
+    coord.SetCoordinateSystemToNormalizedViewport()
+    coord.SetValue(0.9,0.2)
+    colorbar.SetWidth(.075)
+    colorbar.SetHeight(.5)
+
+    # ----------------------------------------
+
+    # backface = vtk.vtkProperty()
+    colors = vtk.vtkNamedColors()
+    # backface.SetColor(colors.GetColor3d("tomato"))
+    # actor.SetBackfaceProperty(backface)
+
+    # ----------------------------------------
+    # camera = renderer.GetActiveCamera()
+    camera = vtk.vtkCamera()
+    # camera.SetPosition(0, 0, 10)
+    # camera.SetPosition(-1, -1, 5)
+    # camera.SetFocalPoint(0, 0, 0)
+    camera.SetViewUp(0, 1, 0)
+    camera.ParallelProjectionOn()
+    # camera.SetParallelScale(7)
+    # renderer.ResetCamera()
+
+    spacing = [1.1,1.1]
+    xc = center[0] + 0.5*(bounds[0] + bounds[1])*spacing[0]
+    yc = center[1] + 0.5*(bounds[2] + bounds[3])*spacing[1]
+    xd = (bounds[1] - bounds[0] + 1)*spacing[0]
+    yd = (bounds[3] - bounds[2] + 1)*spacing[1]
+    d = camera.GetDistance()
+    camera.SetParallelScale(0.5*yd)
+    camera.SetFocalPoint(xc,yc,0.0)
+    camera.SetPosition(xc,yc,+d)
+
+    # ----------------------------------------
+
+    # renderer = vtk.vtkRenderer()
+
+       #
+    light = vtk.vtkLight()
+    light.SetFocalPoint(0,0,0)
+    light.SetPosition(-1, -1, 10)
+
+    # ----------------------------------------
+
+    axis_actor = vtk.vtkCubeAxesActor2D()
+    axis_actor.SetBounds(output.GetBounds())
+    axis_actor.SetLabelFormat("%6.4g")
+    axis_actor.SetFlyModeToOuterEdges()
+    axis_actor.SetFontFactor(1.5)
+    # tprop = vtk.vtkTextProperty()
+    # tprop.SetFontFamilyAsString('Linux Biolinum O')
+    # tprop.SetColor(0,0,0)
+    # tprop.ShadowOn()
+    axis_actor.ScalingOn()
+    axis_actor.SetAxisTitleTextProperty(tprop)
+    axis_actor.SetAxisLabelTextProperty(tprop)
+    axis_actor.SetFlyModeToOuterEdges()
+    # axis_actor.SetFlyModeToClosestTriad()
+    # axis_actor.YAxisVisibilityOn()
+    axis_actor.YAxisVisibilityOff()
+
+
+    # axis_actor = vtk.vtkCubeAxesActor()
+    # axis_actor.SetUseTextActor3D(1)
+    # axis_actor.SetBounds(output.GetBounds())
+    # axis_actor.SetCamera(renderer.GetActiveCamera())
+    # axis_actor.GetTitleTextProperty(0).SetColor(0,0,0)
+    # axis_actor.GetTitleTextProperty(0).SetFontSize(48)
+    # axis_actor.GetLabelTextProperty(0).SetColor(0,0,0)
+    # axis_actor.GetTitleTextProperty(1).SetColor(0,0,0)
+    # axis_actor.GetLabelTextProperty(1).SetColor(0,0,0)
+    # axis_actor.GetTitleTextProperty(2).SetColor(0,0,0)
+    # axis_actor.GetLabelTextProperty(2).SetColor(0,0,0)
+    # axis_actor.DrawXGridlinesOn()
+    # axis_actor.DrawYGridlinesOn()
+    # axis_actor.DrawZGridlinesOn()
+    # axis_actor.SetGridLineLocation(axis_actor.VTK_GRID_LINES_FURTHEST)
+    # axis_actor.XAxisMinorTickVisibilityOff()
+    # axis_actor.YAxisMinorTickVisibilityOff()
+    # axis_actor.ZAxisMinorTickVisibilityOff()
+    # axis_actor.SetFlyModeToStaticEdges()
+
+    # ----------------------------------------
+
+    # # Create a text actor.
+    # txt_actor = vtk.vtkTextActor()
+    # txt_actor.SetInput('t = {:.1e}'.format(simulation_time))
+    # txt_actor.SetTextProperty(tprop)
+    # txt_actor.GetTextProperty().SetFontSize(16)
+    # # txt_actor.SetTextScaleModeToProp()
+    # txt_actor.SetTextScaleModeToViewport()
+    # # txt_actor.SetTextScaleModeToNone()
+    # txt_actor.GetTextProperty().SetJustificationToRight()
+    # txt_actor.GetTextProperty().SetVerticalJustificationToTop()
+    # txt_actor.GetTextProperty().UseTightBoundingBoxOn()
+    #
+    # # txt_actor.SetDisplayPosition(20, 30)
+    # coord = txt_actor.GetPositionCoordinate()
+    # coord.SetCoordinateSystemToNormalizedViewport()
+    # coord.SetValue(0.975,0.975)
+    # # txt_actor.SetWidth(1)
+    # # txt_actor.SetHeight(1)
+
+    # ----------------------------------------
+
+    renderer = vtk.vtkRenderer()
+    renderer.SetActiveCamera(camera)
+    # renderer.GetActiveCamera().Zoom(3.33)
+    # renderer.GetActiveCamera().Dolly(1.0)
+    renderer.AddActor(actor)
+    # renderer.AddActor(colorbar)
+    axis_actor.SetCamera(renderer.GetActiveCamera())
+    renderer.AddActor(axis_actor)
+    # renderer.AddActor(txt_actor)
+    renderer.AddLight(light)
+    # renderer.AddActor(interface_actor)
+    renderer.SetBackground(1, 1, 1)
+    # renderer.SetBackground(colors.GetColor3d("Wheat"))
+    renderer.UseHiddenLineRemovalOn()
+    # renderer.ResetCameraClippingRange()
+
+    # renderer.SetUseDepthPeeling(1)
+    # renderer.SetOcclusionRatio(0.1)
+    # renderer.SetMaximumNumberOfPeels(100)
+    # render_window.SetMultiSamples(0)
+    # render_window.SetAlphaBitPlanes(1)
+
+    # Create the RendererWindow
+    renderer_window = vtk.vtkRenderWindow()
+    renderer_window.SetOffScreenRendering(1)
+    renderer_window.AddRenderer(renderer)
+    # renderer_window.Render()
+    # create a renderwindowinteractor
+    renderer_window.Render()
+    # renderer_window.SetWindowName('Plot')
+    # renderer_window.SetSize(600, 300)
+    # renderer_window.SetSize(500, 500)
+    basesize = 1000
+    renderer_window.SetSize(int(1.1*basesize), basesize)
+    # renderer_window.SetSize(2000, 2000)
+    # win_size = renderer_window.GetSize()
+    # factor = 2
+    # renderer_window.SetSize(factor * win_size[0], factor * win_size[1])
+
+    # ----------------------------------------
+
+    # write plot
+    renderer_window.Render()
+    windowto_image_filter = vtk.vtkWindowToImageFilter()
+    windowto_image_filter.SetInput(renderer_window)
+    rgba = False
+    if rgba:
+        windowto_image_filter.SetInputBufferTypeToRGBA()
+        windowto_image_filter.Update()
+    else:
+        windowto_image_filter.SetInputBufferTypeToRGB()
+        windowto_image_filter.ReadFrontBufferOff()
+        windowto_image_filter.Update()
+
+    # ----------------------------------------
+
+    filename = plotname + '.jpg'
+    writer = vtk.vtkJPEGWriter()
+    # writer = vtk.vtkPNGWriter()
+    writer.SetFileName(filename)
+    writer.SetInputConnection(windowto_image_filter.GetOutputPort())
+    writer.Write()
+
+    return 1
+
+    # Create the RendererWindowInteractor and display the vtk_file
+    # interactor = vtk.vtkRenderWindowInteractor()
+    # interactor.SetRenderWindow(renderer_window)
+    interactor.Initialize()
+    interactor.Start()
+
+    print("--- Done!")
+    return 1
+
+# ------------------------------------------------------------------------------
+
+
+if __name__ == "__main__":
+
+    fn = './data/global_mesh000000.vtu'
+    fn = './data/domain_marker000000.vtu'
+
+    plot3d_vtu(fn, interface_filename=None, plotname='./plots/vtu_plot')
+
+
+    # basepath = '../build/bin'
+    # # basepath = '../build_benchmarks/bin'
+    # # basepath = '.'
+    # paths = os.listdir(basepath)
+    # # print(paths)
+    # folders = [os.path.basename(p) for p in paths if os.path.isdir(
+    #     os.path.join(basepath, p))]
+    # # print(folders)
+    # for idx, f in enumerate(folders):
+    #     print("{}: {}".format(idx, f))
+    # idx = int(input("Choose dataset directory: "))
+    # # idx = 1
+    # main_dir = folders[idx]
+    # benchmark_name = folders[idx]
+    # # main_dir = os.path.join('.', folders[idx])
+    # print("Directory: {}".format(main_dir))
+    # main_dir = os.path.join(basepath, main_dir)
+    # print("Directory: {}".format(main_dir))
+    #
+    # all_paths = []
+    # all_paths.extend(glob(main_dir + '/benchmark*'))
+    # paths = [p for p in all_paths if os.path.isdir(p)]
+    # result = []
+    # for p in paths:
+    #     if any(fname.startswith('statistics') for fname in os.listdir(p)):
+    #         result.append(p)
+    # paths = result
+    # paths = sorted(paths)
+    #
+    # # print(paths)
+    # npaths = len(paths)
+    # print('There are %d benchmark folders.' % npaths)
+    # # sys.exit()
+    #
+    #
+    # # ==========================================================================
+    # # plot grids
+    #
+    # paths = [p for p in all_paths if os.path.isdir(p)]
+    # paths = sorted(paths)
+    # for idx, f in enumerate(paths):
+    #     print("{}: {}".format(idx, os.path.basename(f) ))
+    # idx = int(input("Choose data to plot: "))
+    # # idx = 0
+    # plot_path = paths[idx]
+    # input_directory = os.path.join(plot_path, 'vtk')
+    #
+    # print("Input directory: %s" % input_directory)
+    #
+    # filenames = sorted(glob(input_directory + '/solution*.vt*'))
+    # interface_filenames = sorted(glob(input_directory + '/interface*.vt*'))
+    # print("Found the following files:")
+    # print(filenames)
+    #
+    # cmap = plt.cm.Spectral_r
+    # # cmap = plt.cm.cividis
+    # # cmap = plt.cm.RdBu_r
+    # # cmap = plt.cm.Blues
+    #
+    # plot_output_dir = os.path.join(main_dir, 'plots')
+    # os.makedirs(plot_output_dir, exist_ok=True)
+    # shutil.copy(plot_path + '/log.json', plot_output_dir + '/log.json')
+    #
+    # print("Start plotting...")
+    # output_interval = 1
+    # # output_interval = 50
+    #
+    # file_list = list(zip(filenames[::output_interval], interface_filenames[::output_interval]))
+    #
+    # np.random.shuffle(file_list)
+    #
+    # # for fn in filenames[::output_interval]:
+    # # for fn, fni in zip(filenames[::output_interval], interface_filenames[::output_interval]):
+    #
+    # # file_list = [random.sample(file_list, k=1)]
+    #
+    # for fn, fni in file_list:
+    #     print("Plot file {}, {}".format(fn, fni))
+    #     plotname = os.path.basename(fn)
+    #     plotname = os.path.splitext(plotname)[0]
+    #     plotname = os.path.join(plot_output_dir, plotname)
+    #     # if os.path.isfile(plotname):
+    #     #     continue
+    #     print(plotname)
+    #     plot3d_vtu(fn, interface_filename=fni, plotname=plotname, cmap=cmap)
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic-different-permeabilties.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic-different-permeabilties.py
index 60bb341..6993dbc 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic-different-permeabilties.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic-different-permeabilties.py
@@ -54,10 +54,10 @@ resolutions = {
 # starttimes gives a list of starttimes to run the simulation from.
 # The list is looped over and a simulation is run with t_0 as initial time
 #  for each element t_0 in starttimes.
-starttimes = {0: 0.0}
-# starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9}
+# starttimes = {0: 0.0}
+starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9}
 timestep_size = 0.001
-number_of_timesteps = 1500
+number_of_timesteps = 1
 
 # LDD scheme parameters  ######################################################
 Lw1 = 0.01  # /timestep_size
@@ -69,41 +69,41 @@ Lnw2 = Lw2
 Lw3 = 0.002 # /timestep_size
 Lnw3 = 0.002
 
-Lw4 = 0.002 # /timestep_size
-Lnw4 = 0.002
+Lw4 = 0.0001 # /timestep_size
+Lnw4 = 0.0001
 
-Lw5 = 0.002 # /timestep_size
-Lnw5 = 0.002
+Lw5 = 0.0001 # /timestep_size
+Lnw5 = 0.0001
 
-Lw6 = 0.002  # /timestep_size
-Lnw6 = 0.002
+Lw6 = 0.0001  # /timestep_size
+Lnw6 = 0.0001
 
 lambda12_w = 1
-lambda12_nw = 0.4
+lambda12_nw = 1
 
 lambda23_w = 1
-lambda23_nw = 0.4
+lambda23_nw = 1
 
 lambda24_w = 1
-lambda24_nw= 0.4
+lambda24_nw= 1
 
 lambda25_w= 1
-lambda25_nw= 0.4
+lambda25_nw= 1
 
 lambda34_w = 1
-lambda34_nw = 0.4
+lambda34_nw = 1
 
 lambda36_w = 1
-lambda36_nw = 0.4
+lambda36_nw = 1
 
 lambda45_w = 1
-lambda45_nw = 0.4
+lambda45_nw = 1
 
 lambda46_w = 1
-lambda46_nw = 0.4
+lambda46_nw = 1
 
 lambda56_w = 1
-lambda56_nw = 0.4
+lambda56_nw = 1
 
 include_gravity = False
 debugflag = False
@@ -113,10 +113,10 @@ analyse_condition = False
 # when number_of_timesteps is high, it might take a long time to write all
 # timesteps to disk. Therefore, you can choose to only write data of every
 # plot_timestep_every timestep to disk.
-plot_timestep_every = 3
+plot_timestep_every = 1
 # Decide how many timesteps you want analysed. Analysed means, that
 # subsequent errors of the L-iteration within the timestep are written out.
-number_of_timesteps_to_analyse = 5
+number_of_timesteps_to_analyse = 1
 
 # fine grained control over data to be written to disk in the mesh study case
 # as well as for a regular simuation for a fixed grid.
@@ -225,9 +225,9 @@ porosity = {
     1: 0.57,  #0.2,  # Clayey gravels, clayey sandy gravels
     2: 0.57,  #0.22, # Silty gravels, silty sandy gravels
     3: 0.005,  #0.57, # Clayey sands
-    4: 0.005,  #0.2 # Silty or sandy clay
-    5: 0.005,  #
-    6: 0.005,  #
+    4: 0.0005,  #0.2 # Silty or sandy clay
+    5: 0.0001,  #
+    6: 0.001,  #
 }
 
 # subdom_num : subdomain L for L-scheme
@@ -292,8 +292,8 @@ intrinsic_permeability = {
     1: 0.1,  # sand
     2: 0.1,  # sand, there is a range
     3: 0.001,  #10e-2,  # clay has a range
-    4: 0.001,  #10e-3
-    5: 0.001,  #10e-2,  # clay has a range
+    4: 0.0005,  #10e-3
+    5: 0.0001,  #10e-2,  # clay has a range
     6: 0.001,  #10e-3
 }
 
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic-pure-dd.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic-pure-dd.py
new file mode 100755
index 0000000..52cad94
--- /dev/null
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic-pure-dd.py
@@ -0,0 +1,1015 @@
+#!/usr/bin/python3
+"""Layered soil simulation with inner patch.
+
+This program sets up an LDD simulation
+"""
+
+import dolfin as df
+import sympy as sym
+import functools as ft
+import LDDsimulation as ldd
+import helpers as hlp
+import datetime
+import os
+import pandas as pd
+
+# check if output directory exists
+if not os.path.exists('./output'):
+    os.mkdir('./output')
+    print("Directory ", './output',  " created ")
+else:
+    print("Directory ", './output',  " already exists. Will use as output \
+    directory")
+
+date = datetime.datetime.now()
+datestr = date.strftime("%Y-%m-%d")
+
+
+# init sympy session
+sym.init_printing()
+# solver_tol = 6E-7
+use_case = "TP-R-layered-soil-with-inner-patch-realistic"
+# name of this very file. Needed for log output.
+thisfile = "TP-R-layered_soil_with_inner_patch-realistic-pure-dd.py"
+max_iter_num = 700
+FEM_Lagrange_degree = 1
+mesh_study = False
+resolutions = {
+                # 1: 2e-6,  # h=2
+                # 2: 2e-6,  # h=1.1180
+                # 4: 2e-6,  # h=0.5590
+                # 8: 2e-6,  # h=0.2814
+                16: 8e-6, # h=0.1412
+                # 32: 2e-6,
+                # 64: 2e-6,
+                # 128: 2e-6
+                }
+
+# GRID #######################
+# mesh_resolution = 20
+timestep_size = 0.001
+number_of_timesteps = 5
+plot_timestep_every = 1
+# decide how many timesteps you want analysed. Analysed means, that we write
+# out subsequent errors of the L-iteration within the timestep.
+number_of_timesteps_to_analyse = 5
+starttimes = [0.0]
+
+Lw1 = 0.025  # /timestep_size
+Lnw1 = Lw1
+
+Lw2 = 0.025  # /timestep_size
+Lnw2 = Lw2
+
+Lw3 = 0.025  # /timestep_size
+Lnw3 = Lw3
+
+Lw4 = 0.025  # /timestep_size
+Lnw4 = 0.025
+
+Lw5 = 0.025  # /timestep_size
+Lnw5 = 0.025
+
+Lw6 = 0.025  # /timestep_size
+Lnw6 = Lw6
+
+lambda12_w = 4.0
+lambda12_nw = 4.0
+
+lambda23_w = 4.0
+lambda23_nw = 4.0
+
+lambda24_w = 4.0
+lambda24_nw= 4.0
+
+lambda25_w= 4.0
+lambda25_nw= 4.0
+
+lambda34_w = 4.0
+lambda34_nw = 4.0
+
+lambda36_w = 4.0
+lambda36_nw = 4.0
+
+lambda45_w = 4.0
+lambda45_nw = 4.0
+
+lambda46_w = 4.0
+lambda46_nw = 4.0
+
+lambda56_w = 4.0
+lambda56_nw = 4.0
+
+include_gravity = False
+debugflag = True
+analyse_condition = False
+
+
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+    datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+    )
+
+
+# toggle what should be written to files
+if mesh_study:
+    write_to_file = {
+        'space_errornorms': True,
+        'meshes_and_markers': True,
+        'L_iterations_per_timestep': False,
+        'solutions': True,
+        'absolute_differences': True,
+        'condition_numbers': analyse_condition,
+        'subsequent_errors': True
+    }
+else:
+    write_to_file = {
+        'space_errornorms': True,
+        'meshes_and_markers': True,
+        'L_iterations_per_timestep': False,
+        'solutions': True,
+        'absolute_differences': True,
+        'condition_numbers': analyse_condition,
+        'subsequent_errors': True
+    }
+
+
+# global domain
+subdomain0_vertices = [df.Point(-1.0,-1.0), #
+                        df.Point(1.0,-1.0),#
+                        df.Point(1.0,1.0),#
+                        df.Point(-1.0,1.0)]
+
+interface12_vertices = [df.Point(-1.0, 0.8),
+                        df.Point(0.3, 0.8),
+                        df.Point(0.5, 0.9),
+                        df.Point(0.8, 0.7),
+                        df.Point(1.0, 0.65)]
+
+
+                        # interface23
+interface23_vertices = [df.Point(-1.0, 0.0),
+                        df.Point(-0.35, 0.0),
+                        # df.Point(6.5, 4.5),
+                        df.Point(0.0, 0.0)]
+
+interface24_vertices = [interface23_vertices[2],
+                        df.Point(0.6, 0.0),
+                        ]
+
+interface25_vertices = [interface24_vertices[1],
+                        df.Point(1.0, 0.0)
+                        ]
+
+
+interface32_vertices = [interface23_vertices[2],
+                        interface23_vertices[1],
+                        interface23_vertices[0]]
+
+
+interface36_vertices = [df.Point(-1.0, -0.6),
+                        df.Point(-0.6, -0.45)]
+
+
+interface46_vertices = [interface36_vertices[1],
+                        df.Point(0.3, -0.25)]
+
+interface56_vertices = [interface46_vertices[1],
+                        df.Point(0.65, -0.6),
+                        df.Point(1.0, -0.7)]
+
+
+
+
+interface34_vertices = [interface36_vertices[1],
+                        interface23_vertices[2]]
+# interface36
+
+interface45_vertices = [interface56_vertices[0],
+                        df.Point(0.7, -0.2),
+                        interface25_vertices[0]
+                        ]
+
+
+# interface_vertices introduces a global numbering of interfaces.
+interface_def_points = [interface12_vertices,
+                        interface23_vertices,
+                        interface24_vertices,
+                        interface25_vertices,
+                        interface34_vertices,
+                        interface36_vertices,
+                        interface45_vertices,
+                        interface46_vertices,
+                        interface56_vertices,
+                        ]
+
+adjacent_subdomains = [[1,2],
+                       [2,3],
+                       [2,4],
+                       [2,5],
+                       [3,4],
+                       [3,6],
+                       [4,5],
+                       [4,6],
+                       [5,6]
+                       ]
+
+# subdomain1.
+subdomain1_vertices = [interface12_vertices[0],
+                        interface12_vertices[1],
+                        interface12_vertices[2],
+                        interface12_vertices[3],
+                        interface12_vertices[4],  # southern boundary, 12 interface
+                        subdomain0_vertices[2],  # eastern boundary, outer boundary
+                        subdomain0_vertices[3]]  # northern boundary, outer on_boundary
+
+# vertex coordinates of the outer boundaries. If it can not be specified as a
+# polygon, use an entry per boundary polygon.
+# This information is used for defining
+# the Dirichlet boundary conditions. If a domain is completely internal, the
+# dictionary entry should be 0: None
+subdomain1_outer_boundary_verts = {
+    0: [subdomain1_vertices[4],
+        subdomain1_vertices[5],  # eastern boundary, outer boundary
+        subdomain1_vertices[6],
+        subdomain1_vertices[0]]
+}
+
+#subdomain1
+subdomain2_vertices = [interface23_vertices[0],
+                        interface23_vertices[1],
+                        interface23_vertices[2],
+                        interface24_vertices[1],
+                        interface25_vertices[1],  # southern boundary, 23 interface
+                        subdomain1_vertices[4],  # eastern boundary, outer boundary
+                        subdomain1_vertices[3],
+                        subdomain1_vertices[2],
+                        subdomain1_vertices[1],
+                        subdomain1_vertices[0] ]  # northern boundary, 12 interface
+
+subdomain2_outer_boundary_verts = {
+    0: [subdomain2_vertices[9],
+        subdomain2_vertices[0]],
+    1: [subdomain2_vertices[4],
+        subdomain2_vertices[5]]
+}
+
+
+subdomain3_vertices = [interface36_vertices[0],
+                       interface36_vertices[1],
+                       # interface34_vertices[0],
+                       interface34_vertices[1],
+                       # interface32_vertices[0],
+                       interface32_vertices[1],
+                       interface32_vertices[2]
+                       ]
+
+subdomain3_outer_boundary_verts = {
+    0: [subdomain3_vertices[4],
+        subdomain3_vertices[0]]
+}
+
+
+# subdomain3
+subdomain4_vertices = [interface46_vertices[0],
+                       interface46_vertices[1],
+                       interface45_vertices[1],
+                       interface24_vertices[1],
+                       interface24_vertices[0],
+                       interface34_vertices[1]
+                       ]
+
+subdomain4_outer_boundary_verts = None
+
+subdomain5_vertices = [interface56_vertices[0],
+                       interface56_vertices[1],
+                       interface56_vertices[2],
+                       interface25_vertices[1],
+                       interface25_vertices[0],
+                       interface45_vertices[1],
+                       interface45_vertices[0]
+]
+
+subdomain5_outer_boundary_verts = {
+    0: [subdomain5_vertices[2],
+        subdomain5_vertices[3]]
+}
+
+
+
+subdomain6_vertices = [subdomain0_vertices[0],
+                       subdomain0_vertices[1], # southern boundary, outer boundary
+                       interface56_vertices[2],
+                       interface56_vertices[1],
+                       interface56_vertices[0],
+                       interface36_vertices[1],
+                       interface36_vertices[0]
+                       ]
+
+subdomain6_outer_boundary_verts = {
+    0: [subdomain6_vertices[6],
+        subdomain6_vertices[0],
+        subdomain6_vertices[1],
+        subdomain6_vertices[2]]
+}
+
+
+subdomain_def_points = [subdomain0_vertices,#
+                      subdomain1_vertices,#
+                      subdomain2_vertices,#
+                      subdomain3_vertices,#
+                      subdomain4_vertices,
+                      subdomain5_vertices,
+                      subdomain6_vertices
+                      ]
+
+
+# if a subdomain has no outer boundary write None instead, i.e.
+# i: None
+# if i is the index of the inner subdomain.
+outer_boundary_def_points = {
+    # subdomain number
+    1: subdomain1_outer_boundary_verts,
+    2: subdomain2_outer_boundary_verts,
+    3: subdomain3_outer_boundary_verts,
+    4: subdomain4_outer_boundary_verts,
+    5: subdomain5_outer_boundary_verts,
+    6: subdomain6_outer_boundary_verts
+}
+
+
+isRichards = {
+    1: True,
+    2: True,
+    3: False,
+    4: False,
+    5: False,
+    6: False
+    }
+
+# isRichards = {
+#     1: True,
+#     2: True,
+#     3: True,
+#     4: True,
+#     5: True,
+#     6: True
+#     }
+
+# Dict of the form: { subdom_num : viscosity }
+viscosity = {
+    1: {'wetting' :1,
+         'nonwetting': 1/50},
+    2: {'wetting' :1,
+         'nonwetting': 1/50},
+    3: {'wetting' :1,
+         'nonwetting': 1/50},
+    4: {'wetting' :1,
+         'nonwetting': 1/50},
+    5: {'wetting' :1,
+         'nonwetting': 1/50},
+    6: {'wetting' :1,
+         'nonwetting': 1/50},
+}
+
+# Dict of the form: { subdom_num : density }
+densities = {
+    1: {'wetting': 997.0,  #997
+         'nonwetting': 1.225},  #1},  #1.225},
+    2: {'wetting': 997.0,  #997
+         'nonwetting': 1.225},  #1.225},
+    3: {'wetting': 997.0,  #997
+         'nonwetting': 1.225},  #1.225},
+    4: {'wetting': 997.0,  #997
+         'nonwetting': 1.225},  #1.225}
+    5: {'wetting': 997.0,  #997
+         'nonwetting': 1.225},  #1.225},
+    6: {'wetting': 997.0,  #997
+         'nonwetting': 1.225}  #1.225}
+}
+
+gravity_acceleration = 9.81
+# porosities taken from
+# https://www.geotechdata.info/parameter/soil-porosity.html
+# Dict of the form: { subdom_num : porosity }
+porosity = {
+    1: 0.2,  #0.2,  # Clayey gravels, clayey sandy gravels
+    2: 0.2,  #0.22, # Silty gravels, silty sandy gravels
+    3: 0.2,  #0.37, # Clayey sands
+    4: 0.2,  #0.2 # Silty or sandy clay
+    5: 0.2,  #
+    6: 0.2,  #
+}
+
+# subdom_num : subdomain L for L-scheme
+L = {
+    1: {'wetting' :Lw1,
+         'nonwetting': Lnw1},
+    2: {'wetting' :Lw2,
+         'nonwetting': Lnw2},
+    3: {'wetting' :Lw3,
+         'nonwetting': Lnw3},
+    4: {'wetting' :Lw4,
+         'nonwetting': Lnw4},
+    5: {'wetting' :Lw5,
+         'nonwetting': Lnw5},
+    6: {'wetting' :Lw6,
+         'nonwetting': Lnw6}
+}
+
+
+# interface_num : lambda parameter for the L-scheme on that interface.
+# Note that interfaces are numbered starting from 0, because
+# adjacent_subdomains is a list and not a dict. Historic fuckup, I know
+# We have defined above as interfaces
+# # interface_vertices introduces a global numbering of interfaces.
+# interface_def_points = [interface12_vertices,
+#                         interface23_vertices,
+#                         interface24_vertices,
+#                         interface25_vertices,
+#                         interface34_vertices,
+#                         interface36_vertices,
+#                         interface45_vertices,
+#                         interface46_vertices,
+#                         interface56_vertices,
+#                         ]
+lambda_param = {
+    0: {'wetting': lambda12_w,
+         'nonwetting': lambda12_nw},#
+    1: {'wetting': lambda23_w,
+         'nonwetting': lambda23_nw},#
+    2: {'wetting': lambda24_w,
+         'nonwetting': lambda24_nw},#
+    3: {'wetting': lambda25_w,
+         'nonwetting': lambda25_nw},#
+    4: {'wetting': lambda34_w,
+         'nonwetting': lambda34_nw},#
+    5: {'wetting': lambda36_w,
+         'nonwetting': lambda36_nw},#
+    6: {'wetting': lambda45_w,
+         'nonwetting': lambda45_nw},#
+    7: {'wetting': lambda46_w,
+         'nonwetting': lambda46_nw},#
+    8: {'wetting': lambda56_w,
+         'nonwetting': lambda56_nw},#
+}
+
+
+# after Lewis, see pdf file
+intrinsic_permeability = {
+    1: 0.01,  # sand
+    2: 0.01,  # sand, there is a range
+    3: 0.01,  #10e-2,  # clay has a range
+    4: 0.01,  #10e-3
+    5: 0.01,  #10e-2,  # clay has a range
+    6: 0.01,  #10e-3
+}
+
+
+# relative permeabilty functions on subdomain 1
+def rel_perm1w(s):
+    # relative permeabilty wetting on subdomain1
+    return intrinsic_permeability[1]*s**2
+
+
+def rel_perm1nw(s):
+    # relative permeabilty nonwetting on subdomain1
+    return intrinsic_permeability[1]*(1-s)**2
+
+
+# relative permeabilty functions on subdomain 2
+def rel_perm2w(s):
+    # relative permeabilty wetting on subdomain2
+    return intrinsic_permeability[2]*s**2
+
+
+def rel_perm2nw(s):
+    # relative permeabilty nonwetting on subdomain2
+    return intrinsic_permeability[2]*(1-s)**2
+
+
+# relative permeabilty functions on subdomain 3
+def rel_perm3w(s):
+    # relative permeabilty wetting on subdomain3
+    return intrinsic_permeability[3]*s**3
+
+
+def rel_perm3nw(s):
+    # relative permeabilty nonwetting on subdomain3
+    return intrinsic_permeability[3]*(1-s)**3
+
+
+# relative permeabilty functions on subdomain 4
+def rel_perm4w(s):
+    # relative permeabilty wetting on subdomain4
+    return intrinsic_permeability[4]*s**3
+
+
+def rel_perm4nw(s):
+    # relative permeabilty nonwetting on subdomain4
+    return intrinsic_permeability[4]*(1-s)**3
+
+
+# relative permeabilty functions on subdomain 5
+def rel_perm5w(s):
+    # relative permeabilty wetting on subdomain5
+    return intrinsic_permeability[5]*s**3
+
+
+def rel_perm5nw(s):
+    # relative permeabilty nonwetting on subdomain5
+    return intrinsic_permeability[5]*(1-s)**3
+
+
+# relative permeabilty functions on subdomain 6
+def rel_perm6w(s):
+    # relative permeabilty wetting on subdomain6
+    return intrinsic_permeability[6]*s**3
+
+
+def rel_perm6nw(s):
+    # relative permeabilty nonwetting on subdomain6
+    return intrinsic_permeability[6]*(1-s)**3
+
+
+_rel_perm1w = ft.partial(rel_perm1w)
+_rel_perm1nw = ft.partial(rel_perm1nw)
+
+_rel_perm2w = ft.partial(rel_perm2w)
+_rel_perm2nw = ft.partial(rel_perm2nw)
+
+_rel_perm3w = ft.partial(rel_perm3w)
+_rel_perm3nw = ft.partial(rel_perm3nw)
+
+_rel_perm4w = ft.partial(rel_perm4w)
+_rel_perm4nw = ft.partial(rel_perm4nw)
+
+_rel_perm5w = ft.partial(rel_perm5w)
+_rel_perm5nw = ft.partial(rel_perm5nw)
+
+_rel_perm6w = ft.partial(rel_perm6w)
+_rel_perm6nw = ft.partial(rel_perm6nw)
+
+subdomain1_rel_perm = {
+    'wetting': _rel_perm1w,
+    'nonwetting': _rel_perm1nw
+}
+
+subdomain2_rel_perm = {
+    'wetting': _rel_perm2w,
+    'nonwetting': _rel_perm2nw
+}
+
+subdomain3_rel_perm = {
+    'wetting': _rel_perm3w,
+    'nonwetting': _rel_perm3nw
+}
+
+subdomain4_rel_perm = {
+    'wetting': _rel_perm4w,
+    'nonwetting': _rel_perm4nw
+}
+
+subdomain5_rel_perm = {
+    'wetting': _rel_perm5w,
+    'nonwetting': _rel_perm5nw
+}
+
+subdomain6_rel_perm = {
+    'wetting': _rel_perm6w,
+    'nonwetting': _rel_perm6nw
+}
+
+# dictionary of relative permeabilties on all domains.
+# relative_permeability = {
+#     1: subdomain1_rel_perm,
+#     2: subdomain2_rel_perm,
+#     3: subdomain3_rel_perm,
+#     4: subdomain4_rel_perm,
+#     5: subdomain5_rel_perm,
+#     6: subdomain6_rel_perm
+# }
+relative_permeability = {
+    1: subdomain1_rel_perm,
+    2: subdomain1_rel_perm,
+    3: subdomain1_rel_perm,
+    4: subdomain1_rel_perm,
+    5: subdomain1_rel_perm,
+    6: subdomain1_rel_perm
+}
+
+
+# definition of the derivatives of the relative permeabilities
+# relative permeabilty functions on subdomain 1
+def rel_perm1w_prime(s):
+    # relative permeabilty on subdomain1
+    return intrinsic_permeability[1]*2*s
+
+
+def rel_perm1nw_prime(s):
+    # relative permeabilty on subdomain1
+    return -1.0*intrinsic_permeability[1]*2*(1-s)
+
+
+def rel_perm2w_prime(s):
+    # relative permeabilty on subdomain2
+    return intrinsic_permeability[2]*2*s
+
+
+def rel_perm2nw_prime(s):
+    # relative permeabilty on subdomain2
+    return -1.0*intrinsic_permeability[2]*2*(1-s)
+
+
+# definition of the derivatives of the relative permeabilities
+# relative permeabilty functions on subdomain 3
+def rel_perm3w_prime(s):
+    # relative permeabilty on subdomain3
+    return intrinsic_permeability[3]*3*s**2
+
+
+def rel_perm3nw_prime(s):
+    # relative permeabilty on subdomain3
+    return -1.0*intrinsic_permeability[3]*3*(1-s)**2
+
+
+# definition of the derivatives of the relative permeabilities
+# relative permeabilty functions on subdomain 4
+def rel_perm4w_prime(s):
+    # relative permeabilty on subdomain4
+    return intrinsic_permeability[4]*3*s**2
+
+
+def rel_perm4nw_prime(s):
+    # relative permeabilty on subdomain4
+    return -1.0*intrinsic_permeability[4]*3*(1-s)**2
+
+
+# definition of the derivatives of the relative permeabilities
+# relative permeabilty functions on subdomain 5
+def rel_perm5w_prime(s):
+    # relative permeabilty on subdomain5
+    return intrinsic_permeability[5]*3*s**2
+
+
+def rel_perm5nw_prime(s):
+    # relative permeabilty on subdomain5
+    return -1.0*intrinsic_permeability[5]*3*(1-s)**2
+
+
+# definition of the derivatives of the relative permeabilities
+# relative permeabilty functions on subdomain 6
+def rel_perm6w_prime(s):
+    # relative permeabilty on subdomain6
+    return intrinsic_permeability[6]*3*s**2
+
+
+def rel_perm6nw_prime(s):
+    # relative permeabilty on subdomain6
+    return -1.0*intrinsic_permeability[6]*3*(1-s)**2
+
+
+_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
+_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
+_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
+_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
+_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
+_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
+_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
+_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
+_rel_perm5w_prime = ft.partial(rel_perm5w_prime)
+_rel_perm5nw_prime = ft.partial(rel_perm5nw_prime)
+_rel_perm6w_prime = ft.partial(rel_perm6w_prime)
+_rel_perm6nw_prime = ft.partial(rel_perm6nw_prime)
+
+subdomain1_rel_perm_prime = {
+    'wetting': _rel_perm1w_prime,
+    'nonwetting': _rel_perm1nw_prime
+}
+
+
+subdomain2_rel_perm_prime = {
+    'wetting': _rel_perm2w_prime,
+    'nonwetting': _rel_perm2nw_prime
+}
+
+subdomain3_rel_perm_prime = {
+    'wetting': _rel_perm3w_prime,
+    'nonwetting': _rel_perm3nw_prime
+}
+
+
+subdomain4_rel_perm_prime = {
+    'wetting': _rel_perm4w_prime,
+    'nonwetting': _rel_perm4nw_prime
+}
+
+subdomain5_rel_perm_prime = {
+    'wetting': _rel_perm5w_prime,
+    'nonwetting': _rel_perm5nw_prime
+}
+
+subdomain6_rel_perm_prime = {
+    'wetting': _rel_perm6w_prime,
+    'nonwetting': _rel_perm6nw_prime
+}
+
+
+# dictionary of relative permeabilties on all domains.
+# ka_prime = {
+#     1: subdomain1_rel_perm_prime,
+#     2: subdomain2_rel_perm_prime,
+#     3: subdomain3_rel_perm_prime,
+#     4: subdomain4_rel_perm_prime,
+#     5: subdomain5_rel_perm_prime,
+#     6: subdomain6_rel_perm_prime,
+# }
+ka_prime = {
+    1: subdomain1_rel_perm_prime,
+    2: subdomain1_rel_perm_prime,
+    3: subdomain1_rel_perm_prime,
+    4: subdomain1_rel_perm_prime,
+    5: subdomain1_rel_perm_prime,
+    6: subdomain1_rel_perm_prime,
+}
+
+
+
+# S-pc-relation ship. We use the van Genuchten approach, i.e.
+# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume
+# m = 1-1/n (see Helmig) and assume that residual saturation is Sw
+# this function needs to be monotonically decreasing in the capillary pressure
+# pc.
+# Since in the richards case pc=-pw, this becomes as a function of pw a mono
+# tonically INCREASING function like in our Richards-Richards paper. However
+# since we unify the treatment in the code for Richards and two-phase, we need
+# the same requierment
+# for both cases, two-phase and Richards.
+# def saturation(pc, n_index, alpha):
+#     # inverse capillary pressure-saturation-relationship
+#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
+#
+# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
+# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
+# def saturation_sym(pc, n_index, alpha):
+#     # inverse capillary pressure-saturation-relationship
+#     #df.conditional(pc > 0,
+#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
+#
+#
+# # derivative of S-pc relationship with respect to pc. This is needed for the
+# # construction of a analytic solution.
+# def saturation_sym_prime(pc, n_index, alpha):
+#     # inverse capillary pressure-saturation-relationship
+#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
+##
+# # note that the conditional definition of S-pc in the nonsymbolic part will be
+# # incorporated in the construction of the exact solution below.
+# S_pc_sym = {
+#     1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+#     2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+#     3: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+#     4: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+#     5: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+#     6: ft.partial(saturation_sym, n_index=3, alpha=0.001)
+# }
+#
+# S_pc_sym_prime = {
+#     1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+#     2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+#     3: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+#     4: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+#     5: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+#     6: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001)
+# }
+#
+# sat_pressure_relationship = {
+#     1: ft.partial(saturation, n_index=3, alpha=0.001),
+#     2: ft.partial(saturation, n_index=3, alpha=0.001),
+#     3: ft.partial(saturation, n_index=3, alpha=0.001),
+#     4: ft.partial(saturation, n_index=3, alpha=0.001),
+#     5: ft.partial(saturation, n_index=3, alpha=0.001),
+#     6: ft.partial(saturation, n_index=3, alpha=0.001)
+# }
+
+def saturation(pc, n_index):
+    # inverse capillary pressure-saturation-relationship
+    return df.conditional(pc > 0, 1/((1 + pc)**(1/(n_index + 1))), 1)
+
+
+def saturation_sym(pc, n_index):
+    # inverse capillary pressure-saturation-relationship
+    return 1/((1 + pc)**(1/(n_index + 1)))
+
+
+def saturation_sym_prime(pc, n_index):
+    # inverse capillary pressure-saturation-relationship
+    return -1/((n_index+1)*(1 + pc)**((n_index+2)/(n_index+1)))
+
+
+S_pc_sym = {
+    1: ft.partial(saturation_sym, n_index=1),
+    2: ft.partial(saturation_sym, n_index=1),
+    3: ft.partial(saturation_sym, n_index=1),
+    4: ft.partial(saturation_sym, n_index=1),
+    5: ft.partial(saturation_sym, n_index=1),
+    6: ft.partial(saturation_sym, n_index=1)
+}
+
+S_pc_sym_prime = {
+    1: ft.partial(saturation_sym_prime, n_index=1),
+    2: ft.partial(saturation_sym_prime, n_index=1),
+    3: ft.partial(saturation_sym_prime, n_index=1),
+    4: ft.partial(saturation_sym_prime, n_index=1),
+    5: ft.partial(saturation_sym_prime, n_index=1),
+    6: ft.partial(saturation_sym_prime, n_index=1)
+}
+
+sat_pressure_relationship = {
+    1: ft.partial(saturation, n_index=1),
+    2: ft.partial(saturation, n_index=1),
+    3: ft.partial(saturation, n_index=1),
+    4: ft.partial(saturation, n_index=1),
+    5: ft.partial(saturation, n_index=1),
+    6: ft.partial(saturation, n_index=1)
+}
+
+#############################################
+# Manufacture source expressions with sympy #
+#############################################
+x, y = sym.symbols('x[0], x[1]')  # needed by UFL
+t = sym.symbols('t', positive=True)
+
+
+p_e_sym = {
+    1: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+        'nonwetting': 0*t },
+    2: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+        'nonwetting': 0*t },
+    3: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+        'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
+    4: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+        'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
+    5: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+        'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
+    6: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+        'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
+}
+
+pc_e_sym = dict()
+for subdomain, isR in isRichards.items():
+    if isR:
+        pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
+    else:
+        pc_e_sym.update(
+            {subdomain: p_e_sym[subdomain]['nonwetting']
+                - p_e_sym[subdomain]['wetting']}
+            )
+
+
+symbols = {"x": x,
+           "y": y,
+           "t": t}
+# turn above symbolic code into exact solution for dolphin and
+# construct the rhs that matches the above exact solution.
+exact_solution_example = hlp.generate_exact_solution_expressions(
+                        symbols=symbols,
+                        isRichards=isRichards,
+                        symbolic_pressure=p_e_sym,
+                        symbolic_capillary_pressure=pc_e_sym,
+                        saturation_pressure_relationship=S_pc_sym,
+                        saturation_pressure_relationship_prime=S_pc_sym_prime,
+                        viscosity=viscosity,
+                        porosity=porosity,
+                        relative_permeability=relative_permeability,
+                        relative_permeability_prime=ka_prime,
+                        densities=densities,
+                        gravity_acceleration=gravity_acceleration,
+                        include_gravity=include_gravity,
+                        )
+source_expression = exact_solution_example['source']
+exact_solution = exact_solution_example['exact_solution']
+initial_condition = exact_solution_example['initial_condition']
+
+# Dictionary of dirichlet boundary conditions.
+dirichletBC = dict()
+# similarly to the outer boundary dictionary, if a patch has no outer boundary
+# None should be written instead of an expression.
+# This is a bit of a brainfuck:
+# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind.
+# Since a domain patch can have several disjoint outer boundary parts, the
+# expressions need to get an enumaration index which starts at 0.
+# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of
+# subdomain ind and boundary part j.
+# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting']
+# return the actual expression needed for the dirichlet condition for both
+# phases if present.
+
+# subdomain index: {outer boudary part index: {phase: expression}}
+for subdomain in isRichards.keys():
+    # if subdomain has no outer boundary, outer_boundary_def_points[subdomain]
+    # is None
+    if outer_boundary_def_points[subdomain] is None:
+        dirichletBC.update({subdomain: None})
+    else:
+        dirichletBC.update({subdomain: dict()})
+        # set the dirichlet conditions to be the same code as exact solution on
+        # the subdomain.
+        for outer_boundary_ind in outer_boundary_def_points[subdomain].keys():
+            dirichletBC[subdomain].update(
+                {outer_boundary_ind: exact_solution[subdomain]}
+                )
+
+
+# LOG FILE OUTPUT #############################################################
+# read this file and print it to std out. This way the simulation can produce a
+# log file with ./TP-R-layered_soil.py | tee simulation.log
+f = open(thisfile, 'r')
+print(f.read())
+f.close()
+
+
+# RUN #########################################################################
+for starttime in starttimes:
+    for mesh_resolution, solver_tol in resolutions.items():
+        # initialise LDD simulation class
+        simulation = ldd.LDDsimulation(
+            tol=1E-14,
+            LDDsolver_tol=solver_tol,
+            debug=debugflag,
+            max_iter_num=max_iter_num,
+            FEM_Lagrange_degree=FEM_Lagrange_degree,
+            mesh_study=mesh_study
+            )
+
+        simulation.set_parameters(
+            use_case=use_case,
+            output_dir=output_string,
+            subdomain_def_points=subdomain_def_points,
+            isRichards=isRichards,
+            interface_def_points=interface_def_points,
+            outer_boundary_def_points=outer_boundary_def_points,
+            adjacent_subdomains=adjacent_subdomains,
+            mesh_resolution=mesh_resolution,
+            viscosity=viscosity,
+            porosity=porosity,
+            L=L,
+            lambda_param=lambda_param,
+            relative_permeability=relative_permeability,
+            saturation=sat_pressure_relationship,
+            starttime=starttime,
+            number_of_timesteps=number_of_timesteps,
+            number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
+            plot_timestep_every=plot_timestep_every,
+            timestep_size=timestep_size,
+            sources=source_expression,
+            initial_conditions=initial_condition,
+            dirichletBC_expression_strings=dirichletBC,
+            exact_solution=exact_solution,
+            densities=densities,
+            include_gravity=include_gravity,
+            gravity_acceleration=gravity_acceleration,
+            write2file=write_to_file,
+            )
+
+        simulation.initialise()
+        output_dir = simulation.output_dir
+        # simulation.write_exact_solution_to_xdmf()
+        output = simulation.run(analyse_condition=analyse_condition)
+        for subdomain_index, subdomain_output in output.items():
+            mesh_h = subdomain_output['mesh_size']
+            for phase, error_dict in subdomain_output['errornorm'].items():
+                filename = output_dir \
+                    + "subdomain{}".format(subdomain_index)\
+                    + "-space-time-errornorm-{}-phase.csv".format(phase)
+                # for errortype, errornorm in error_dict.items():
+
+                # eocfile = open("eoc_filename", "a")
+                # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
+                # eocfile.close()
+                # if subdomain.isRichards:mesh_h
+                data_dict = {
+                    'mesh_parameter': mesh_resolution,
+                    'mesh_h': mesh_h,
+                }
+                for norm_type, errornorm in error_dict.items():
+                    data_dict.update(
+                        {norm_type: errornorm}
+                    )
+                errors = pd.DataFrame(data_dict, index=[mesh_resolution])
+                # check if file exists
+                if os.path.isfile(filename) is True:
+                    with open(filename, 'a') as f:
+                        errors.to_csv(
+                            f,
+                            header=False,
+                            sep='\t',
+                            encoding='utf-8',
+                            index=False
+                            )
+                else:
+                    errors.to_csv(
+                        filename,
+                        sep='\t',
+                        encoding='utf-8',
+                        index=False
+                        )
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
index 7a512be..6b336df 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
@@ -61,7 +61,6 @@ timestep_size = 0.001
 number_of_timesteps = 1000
 
 # LDD scheme parameters  ######################################################
-
 Lw1 = 0.01  # /timestep_size
 Lnw1 = Lw1
 
-- 
GitLab