diff --git a/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
similarity index 100%
rename from LDDsimulation.py
rename to LDDsimulation/LDDsimulation.py
diff --git a/LDDsimulation/__init__.py b/LDDsimulation/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/domainPatch.py b/LDDsimulation/domainPatch.py
similarity index 99%
rename from domainPatch.py
rename to LDDsimulation/domainPatch.py
index 1f69ce34b2d40a73470a3ccaf328f79184e5678e..47d72f150860bb9b842e4eefef2ec368ff36a6fd 100644
--- a/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -142,6 +142,8 @@ class DomainPatch(df.SubDomain):
         # dictionary holding the dof indices corresponding to an interface of
         # given interface. see self._calc_dof_indices_of_interfaces()
         self._dof_indices_of_interface = dict()
+        # List of objects of clas outerBoundary initialised by self._init_markers()
+        # can be NONE if no outer boundary is present.
         self.outer_boundary = []
         # measures are set by self._init_measures()
         self.iteration_number = 0
diff --git a/helpers.py b/LDDsimulation/helpers.py
similarity index 100%
rename from helpers.py
rename to LDDsimulation/helpers.py
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
new file mode 100755
index 0000000000000000000000000000000000000000..12938f2c5825291b29f0ef72cae1b85cd2da7d81
--- /dev/null
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -0,0 +1,530 @@
+#!/usr/bin/python3
+import dolfin as df
+import mshr
+import numpy as np
+import sympy as sym
+import typing as tp
+import domainPatch as dp
+import LDDsimulation as ldd
+import functools as ft
+#import ufl as ufl
+
+##### Domain and Interface ####
+# global simulation domain domain
+sub_domain0_vertices = [df.Point(0.0,0.0), #
+                        df.Point(1.0,0.0),#
+                        df.Point(1.0,1.0),#
+                        df.Point(0.0,1.0)]
+# interface between subdomain1 and subdomain2
+interface12_vertices = [df.Point(0.0, 0.5),
+                        df.Point(1.0, 0.5) ]
+# subdomain1.
+sub_domain1_vertices = [interface12_vertices[0],
+                        interface12_vertices[1],
+                        df.Point(1.0,1.0),
+                        df.Point(0.0,1.0) ]
+
+# 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: [interface12_vertices[0], #
+        df.Point(0.0,1.0), #
+        df.Point(1.0,1.0), #
+        interface12_vertices[1]]
+}
+# subdomain2
+sub_domain2_vertices = [df.Point(0.0,0.0),
+                        df.Point(1.0,0.0),
+                        interface12_vertices[1],
+                        interface12_vertices[0] ]
+
+subdomain2_outer_boundary_verts = {
+    0: [interface12_vertices[1], #
+        df.Point(1.0,0.0), #
+        df.Point(0.0,0.0), #
+        interface12_vertices[0]]
+}
+# subdomain2_outer_boundary_verts = {
+#     0: [interface12_vertices[0], df.Point(0.0,0.0)],#
+#     1: [df.Point(0.0,0.0), df.Point(1.0,0.0)], #
+#     2: [df.Point(1.0,0.0), interface12_vertices[1]]
+# }
+# subdomain2_outer_boundary_verts = {
+#     0: None
+# }
+
+# list of subdomains given by the boundary polygon vertices.
+# Subdomains are given as a list of dolfin points forming
+# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
+# to create the subdomain. subdomain_def_points[0] contains the
+# vertices of the global simulation domain and subdomain_def_points[i] contains the
+# vertices of the subdomain i.
+subdomain_def_points = [sub_domain0_vertices,#
+                      sub_domain1_vertices,#
+                      sub_domain2_vertices]
+# in the below list, index 0 corresponds to the 12 interface which has index 1
+interface_def_points = [interface12_vertices]
+
+outer_boundary_def_points = {
+    # subdomain number
+    1 : subdomain1_outer_boundary_verts,
+    2 : subdomain2_outer_boundary_verts
+}
+
+# adjacent_subdomains[i] contains the indices of the subdomains sharing the
+# interface i (i.e. given by interface_def_points[i]).
+adjacent_subdomains = [[1,2]]
+isRichards = {
+    1: True, #
+    2: True
+    }
+
+Tmax = 2
+dt1 = 0.1
+dt2 = dt1
+# the timestep is taken as a dictionary to be able to run different timesteps on
+# different subdomains.
+timestep_size = {
+    1: dt1,
+    2: dt2
+}
+
+viscosity = {#
+# subdom_num : viscosity
+    1 : {'wetting' :1}, #
+    2 : {'wetting' :1}
+}
+
+porosity = {#
+# subdom_num : porosity
+    1 : 1,#
+    2 : 1
+}
+
+L = {#
+# subdom_num : subdomain L for L-scheme
+    1 : {'wetting' :0.25},#
+    2 : {'wetting' :0.25}
+}
+
+lambda_param = {#
+# subdom_num : lambda parameter for the L-scheme
+    1 : {'wetting' :4},#
+    2 : {'wetting' :4}
+}
+
+## relative permeabilty functions on subdomain 1
+def rel_perm1(s):
+    # relative permeabilty on subdomain1
+    return s**2
+
+_rel_perm1 = ft.partial(rel_perm1)
+
+subdomain1_rel_perm = {
+    'wetting': _rel_perm1,#
+    'nonwetting': None
+}
+## relative permeabilty functions on subdomain 2
+def rel_perm2(s):
+    # relative permeabilty on subdomain2
+    return s**3
+_rel_perm2 = ft.partial(rel_perm2)
+
+subdomain2_rel_perm = {
+    'wetting': _rel_perm2,#
+    'nonwetting': None
+}
+
+## dictionary of relative permeabilties on all domains.
+relative_permeability = {#
+    1: subdomain1_rel_perm,
+    2: subdomain2_rel_perm
+}
+
+def saturation(pressure, subdomain_index):
+    # inverse capillary pressure-saturation-relationship
+    return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+
+sat_pressure_relationship = {#
+    1: ft.partial(saturation, subdomain_index = 1),#
+    2: ft.partial(saturation, subdomain_index = 2)
+}
+
+source_expression = {
+    1: {'wetting': '4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )'},
+    2: {'wetting': '2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))'}
+}
+
+initial_condition = {
+    1: {'wetting': '-(x[0]*x[0] + x[1]*x[1])'},#
+    2: {'wetting': '-x[1]*x[1]'}
+}
+
+exact_solution = {
+    1: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'},#
+    2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}
+}
+
+dirichletBC = exact_solution
+
+# def saturation(pressure, subdomain_index):
+#     # inverse capillary pressure-saturation-relationship
+#     return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+#
+# sa
+
+mesh_resolution = 3
+
+# initialise LDD simulation class
+simulation = ldd.LDDsimulation()
+simulation.set_parameters(output_dir = "",#
+    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,#
+    timestep_size = timestep_size,#
+    sources = source_expression,#
+    initial_conditions = initial_condition,#
+    outer_dirichletBC = dirichletBC,#
+    )
+
+simulation.initialise()
+# simulation._init_meshes_and_markers(subdomain_def_points, mesh_resolution=2)
+# subdomain marker functions
+domain_marker = simulation.domain_marker
+mesh_subdomain = simulation.mesh_subdomain
+
+# mesh = mesh_subdomain[0]
+# interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
+# interface_marker.set_all(0)
+# interface = dp.Interface(vertices=interface_def_points[0], #
+#                                     tol=mesh.hmin()/100,#
+#                                     internal=True,#
+#                                     adjacent_subdomains = adjacent_subdomains[0])
+# print("interface vertices are: \n", interface.vertices)
+# interface.mark(interface_marker, 1)
+# simulation._init_interfaces(interface_def_points, adjacent_subdomains)
+
+interface = simulation.interface
+interface_marker = simulation.interface_marker
+
+subdoms = simulation.subdomain
+
+for iterations in range(1):
+    #
+
+    a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
+    rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
+    rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
+    extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
+    splitrhs = df.assemble(rhssplit1)
+    print("splitrhs: \n", splitrhs.get_local())
+    extrat_assembled = df.assemble(extratrem)
+    print("extratrem: \n", extrat_assembled.get_local())
+    added = splitrhs + extrat_assembled
+    print("both added together:\n", added.get_local())
+
+    A1 = df.assemble(a1)
+    dsform = (timestep*lambda1*u1*v1)*ds1(1)
+    dsform_vec = df.assemble(dsform)
+    b1 = df.assemble(rhs1)
+    #p_gamma1.apply(A1,b1)
+    # outerBC1.apply(A1,b1)
+    u  = df.Function(V1)
+    U1 = u.vector()
+    df.solve(A1,U1,b1)
+    # print("Form:\n", A1.array())
+    print("RHS:\n", b1.get_local())
+    #print("ds term:\n", dsform_vec.array())
+    # print('solution:\n', U1.get_local())
+    p1_i.vector()[:] = U1
+    interface[0].read_pressure_dofs(from_function = p1_i, #
+                            interface_dofs = dofs_on_interface1,#
+                            dof_to_vert_map = dom1_d2v,#
+                            local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+                            phase = 'wetting',#
+                            subdomain_ind = 1)
+#
+# Save mesh to file
+df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
+df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
+df.File('./test_global_interface_marker.pvd') << interface_marker
+df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
+df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
+df.File('./test_domain_marker.pvd') << domain_marker
+df.File('./test_domain_layered_soil_solution.pvd') << u
+df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
+df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
+df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
+df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
+
+# boundary_marker1 = simulation.subdomain[1].outer_boundary_marker
+# boundary_marker2 = simulation.subdomain[2].outer_boundary_marker
+# dx_1 = simulation.subdomain[1].dx
+# dx_2 = simulation.subdomain[2].dx
+# ds_1 = simulation.subdomain[1].ds
+# ds_2 = simulation.subdomain[2].ds
+#
+# interface_marker1 = simulation.subdomain[1].interface_marker#df.MeshFunction('size_t', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
+# interface_marker2 = simulation.subdomain[2].interface_marker#df.MeshFunction('size_t', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
+# # interface_marker1.set_all(0)
+# # interface_marker2.set_all(0)
+# # print(dir(interface_marker1))
+#
+# interface[0].mark(interface_marker1, 1)
+# interface[0].mark(interface_marker2, 1)
+# # interface_marker1 = simulation.subdomain[1].boundary_marker
+#
+# # domain measures
+#
+# # dx1 = simulation.subdomain[1].dx
+# dx1 = df.Measure('dx', domain = mesh_subdomain[1])
+# dx2 = df.Measure('dx', domain = mesh_subdomain[2])
+#
+# # boundary measures
+# ds1 = df.Measure('ds', domain = mesh_subdomain[1], subdomain_data = interface_marker1)
+# ds2 = df.Measure('ds', domain = mesh_subdomain[2], subdomain_data = interface_marker2)
+#
+# V1 = df.FunctionSpace(mesh_subdomain[1], 'P', 1)
+# p1_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2, t = 0)
+#
+# coordinates = mesh_subdomain[0].coordinates()
+# print("global mesh coordinates of domain: \n", coordinates)
+# coordinates1 = mesh_subdomain[1].coordinates()
+#
+# coordinates2 = mesh_subdomain[2].coordinates()
+# submesh1_data = mesh_subdomain[1].data()
+# dom1_parent_vertex_indices = submesh1_data.array('parent_vertex_indices',0)
+# # print("map of dom1 vertices to parent vertices on global domain:\n", dom1_parent_vertex_indices)
+# #
+# # coordinates2 = mesh_subdomain[2].coordinates()
+# # print("\nvertices of subdomain2: \n", coordinates2)
+# submesh2_data = mesh_subdomain[2].data()
+# dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
+# # print("'parent_vertex_indices',1 \n")
+# #
+# #
+# # print("on domain: \n")
+# # interface_coordinates = interface[0].coordinates(interface_marker, 1, print_coordinates = True)
+# # interface_coordinates = interface[0]._vertex_indices(interface_marker, 1, print_vertex_indices = True)
+# # print("\non subdomain2: \n")
+# # dom2_interface_coordinates = interface[0].coordinates(interface_marker2, 1, print_coordinates = True)
+# dom2_interface_def_points = interface[0]._vertex_indices(interface_marker2, 1, print_vertex_indices = True)
+#
+#
+#
+# dofs_on_interface1 = interface[0].dofs_on_interface(V1, interface_marker1, 1)
+# print("Dofs on Interface dom1", dofs_on_interface1)
+#
+# V2 = df.FunctionSpace(mesh_subdomain[2], 'P', 1)
+# dofs_on_interface2 = interface[0].dofs_on_interface(V2, interface_marker2, 1)
+# # print("Dofs on Interface dom1", dofs_on_interface2)
+#
+# # markermesh = interface_marker2.mesh()
+# # functionSpaceMesh = V2.mesh()
+# # print("Domain 2 mesh extracted from markerfunction:\n", markermesh.coordinates())
+# # print("Domain 2 mesh extracted from Function Space V2:\n", functionSpaceMesh.coordinates())
+#
+# dom1_d2v = df.dof_to_vertex_map(V1)
+# dom1_v2d = df.vertex_to_dof_map(V1)
+# # print("dof to vertex map for V1:\n", dom1_d2v)
+# # print("vertex to dof map for V1:\n", dom1_v2d)
+# dom2_d2v = df.dof_to_vertex_map(V2)
+# dom2_v2d = df.vertex_to_dof_map(V2)
+# # print("dof to vertex map for V2:\n", dom2_d2v)
+# # print("vertex to dof map for V2:\n", dom2_v2d)
+#
+# testf1 = df.Constant(41)
+# f1test = df.interpolate(testf1, V1)
+# testf2 = df.Constant(42)
+# f2test = df.interpolate(testf2, V2)
+#
+# # f1vec = f1test.compute_vertex_values()
+# # for i in range(len(f1vec)):
+# #     f1vec[i] = i
+#
+# # f1test.vector()[:] = f1vec
+# # f1nodal = f1test.vector().get_local()
+# # for i in range(len(f1nodal)):
+# #     f1nodal[i] = i
+#
+# # print("nodal values of f1nodal\n", f1nodal)
+#
+# #
+# # print('dofs on interface \n', f1vec[dofs_on_interface])
+# # f1test.vector()[:] = f1vec
+# # print("f1test.vector()[:]\n", f1test.vector()[:])
+# # print("f1test.vector()[dofs_on_interface]\n", f1test.vector()[dofs_on_interface])
+#
+# # print("constant function on dom1\n", f1vec)
+# #
+# # for i in dofs_on_interface1:
+# #     print("local dof ind %i, parent = %i", i,dom1_parent_vertex_indices[dom1_d2v[i]])
+# #     # f1test.vector() is enumerated by dof not by vertex
+# #     f1test.vector()[i] = i
+# # # this enumerates according to the dof numbering
+# # print("testfunction f1test with set interface values\n", f1test.vector()[:])
+# # print("testfunction f1test with get_local\n", f1test.vector().get_local())
+# # f1test.vector()[dom1_v2d[dom1_d2v[dofs_on_interface]]] = f1vec[dofs_on_interface]
+# # .vector()[:] and .vector().get_local() order values according to the vertex numbering
+# # print("testfunction f1test with compute_vertex_values\n", f1test.compute_vertex_values())
+# # print("testfunction f1test with compute_vertex_values mit d2vmap\n", f1test.compute_vertex_values()[dom1_d2v[:]])
+# # i=0
+# # for x in coordinates1:
+# #      print(f"vertex {i}: f1test[{dom1_v2d[i]}] = {f1test.vector()[dom1_v2d[i]]}\tu({x}) = {f1test(x)}")
+# #      i += 1
+# # print("\non subdomain1: \n")
+# print("vertices of subdomain1: \n", coordinates1)
+# # dom1_interface_coordinates = interface[0].coordinates(interface_marker1, 1, print_coordinates = True)
+# dom1_interface_def_points = interface[0]._vertex_indices(interface_marker1, 1, print_vertex_indices = True)
+#
+# interface[0].read_pressure_dofs(from_function = f1test, #
+#                         interface_dofs = dofs_on_interface1,#
+#                         dof_to_vert_map = dom1_d2v,#
+#                         local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+#                         phase = 'wetting',
+#                         subdomain_ind = 1)
+#
+# interface[0].write_pressure_dofs(to_function = f2test, #
+#                         interface_dofs = dofs_on_interface2,#
+#                         dof_to_vert_map = dom2_d2v,#
+#                         local_to_parent_vertex_map = dom2_parent_vertex_indices,#
+#                         phase = 'wetting',
+#                         subdomain_ind = 2)
+#
+# i=0
+# for x in coordinates2:
+#      print(f"vertex {i}: f2test[{dom2_v2d[i]}] = {f2test.vector()[dom2_v2d[i]]}\tu({x}) = {f2test(x)}")
+#      i += 1
+# print("\non subdomain2: \n")
+# print("vertices of subdomain2: \n", coordinates2)
+#
+# # dom2_interface_coordinates = interface[0].coordinates(interface_marker2, 1, print_coordinates = True)
+# dom2_interface_def_points = interface[0]._vertex_indices(interface_marker2, 1, print_vertex_indices = True)
+# print("Parent indices", dom2_parent_vertex_indices)
+#
+#
+# u1 = df.TrialFunction(V1)
+# v1 = df.TestFunction(V1)
+#
+# u2 = df.TrialFunction(V2)
+# v2 = df.TestFunction(V2)
+#
+# def relative_permeability(s, subdomain_index):
+#     if subdomain_index == 1:
+#         # relative permeabilty on subdomain1
+#         return s**2
+#     else:
+#         # relative permeabilty on subdomain2
+#         return s**3
+#
+# def saturation(pressure, subdomain_index):
+#     # inverse capillary pressure-saturation-relationship
+#     return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+#
+#
+# # exact solution
+# # cell = df.Cell("triangle", 2)
+# # x1 = df.SpatialCoordinate(mesh_subdomain[1])
+# # x2 = df.SpatialCoordinate(mesh_subdomain[2])
+#
+# p2_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])', domain = mesh_subdomain[2], degree = 2, t = 0)
+# p1_e = df.interpolate(p1_exact, V1)
+# p2_e = df.interpolate(p2_exact, V2)
+# #
+# # initial conditions
+# p1_initial = df.Expression('-(x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2)
+# p2_initial = df.Expression('-x[1]*x[1]', domain = mesh_subdomain[2], degree = 2)
+# p1_0 = df.interpolate(p1_initial, V1)
+# p2_0 = df.interpolate(p2_initial, V2)
+#
+# # Initial boundary value for the pressure on interface
+# p_gamma1 = df.DirichletBC(V1, p1_exact, interface_marker1, 1)
+# p_gamma2 = df.DirichletBC(V2, p2_exact, interface_marker2, 1)
+#
+# outerBC1 = df.DirichletBC(V1, p1_exact, boundary_marker1, 1)
+# outerBC2 = df.DirichletBC(V2, p2_exact, boundary_marker2, 1)
+#
+#
+# ### source terms
+# # source1 = df.Expression('4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )', domain = mesh_subdomain[1], degree = 2, t = 0)
+# # source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
+# source1 = df.Expression(source_expression[1]['wetting'], domain = mesh_subdomain[1], degree = 2, t = 0)
+# source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
+#
+# source1h = df.interpolate(source1, V1)
+# source2h = df.interpolate(source2, V2)
+#
+# L1 = 0.25
+# L2 = L1
+# timestep = 0.1
+# lambda1 = 4
+# lambda2 = lambda1
+# # initial iteartion
+# p1_i = p1_0
+# p2_i = p2_0
+# # initial g_i
+# n1 = df.FacetNormal(mesh_subdomain[1])
+# flux = -df.grad(p1_0)
+# gl0 = (df.dot(flux, n1) - lambda1*p1_0)*v1*ds1(1)
+# print('gl0 term assembliert', df.assemble(gl0).get_local())
+# g1_i = -lambda1*p1_i
+# g2_i = -lambda1*p2_i
+
+# for iterations in range(1):
+#     #
+#
+#     a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
+#     rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
+#     rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
+#     extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
+#     splitrhs = df.assemble(rhssplit1)
+#     print("splitrhs: \n", splitrhs.get_local())
+#     extrat_assembled = df.assemble(extratrem)
+#     print("extratrem: \n", extrat_assembled.get_local())
+#     added = splitrhs + extrat_assembled
+#     print("both added together:\n", added.get_local())
+#
+#     A1 = df.assemble(a1)
+#     dsform = (timestep*lambda1*u1*v1)*ds1(1)
+#     dsform_vec = df.assemble(dsform)
+#     b1 = df.assemble(rhs1)
+#     #p_gamma1.apply(A1,b1)
+#     # outerBC1.apply(A1,b1)
+#     u  = df.Function(V1)
+#     U1 = u.vector()
+#     df.solve(A1,U1,b1)
+#     # print("Form:\n", A1.array())
+#     print("RHS:\n", b1.get_local())
+#     #print("ds term:\n", dsform_vec.array())
+#     # print('solution:\n', U1.get_local())
+#     p1_i.vector()[:] = U1
+#     interface[0].read_pressure_dofs(from_function = p1_i, #
+#                             interface_dofs = dofs_on_interface1,#
+#                             dof_to_vert_map = dom1_d2v,#
+#                             local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+#                             phase = 'wetting',#
+#                             subdomain_ind = 1)
+# #
+# # Save mesh to file
+# df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
+# df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
+# df.File('./test_global_interface_marker.pvd') << interface_marker
+# df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
+# df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
+# df.File('./test_domain_marker.pvd') << domain_marker
+# df.File('./test_domain_layered_soil_solution.pvd') << u
+# df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
+# df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
+# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
+# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
deleted file mode 100755
index df0fb8065f30b0106ca23d31977b6e5e04d366d1..0000000000000000000000000000000000000000
--- a/RR-2-patch-test.py
+++ /dev/null
@@ -1,478 +0,0 @@
-#!/usr/bin/python3
-import dolfin as df
-import mshr
-import numpy as np
-import sympy as sym
-import typing as tp
-import domainPatch as dp
-import LDDsimulation as ldd
-import functools as ft
-#import ufl as ufl
-
-##### Domain and Interface ####
-# global simulation domain domain
-sub_domain0_vertices = [df.Point(0.0,0.0), #
-                        df.Point(1.0,0.0),#
-                        df.Point(1.0,1.0),#
-                        df.Point(0.0,1.0)]
-# interface between subdomain1 and subdomain2
-interface12_vertices = [df.Point(0.0, 0.5),
-                        df.Point(1.0, 0.5) ]
-# subdomain1.
-sub_domain1_vertices = [interface12_vertices[0],
-                        interface12_vertices[1],
-                        df.Point(1.0,1.0),
-                        df.Point(0.0,1.0) ]
-
-# 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: [interface12_vertices[0], #
-        df.Point(0.0,1.0), #
-        df.Point(1.0,1.0), #
-        interface12_vertices[1]]
-}
-# subdomain2
-sub_domain2_vertices = [df.Point(0.0,0.0),
-                        df.Point(1.0,0.0),
-                        interface12_vertices[1],
-                        interface12_vertices[0] ]
-
-subdomain2_outer_boundary_verts = {
-    0: [interface12_vertices[1], #
-        df.Point(1.0,0.0), #
-        df.Point(0.0,0.0), #
-        interface12_vertices[0]]
-}
-# subdomain2_outer_boundary_verts = {
-#     0: [interface12_vertices[0], df.Point(0.0,0.0)],#
-#     1: [df.Point(0.0,0.0), df.Point(1.0,0.0)], #
-#     2: [df.Point(1.0,0.0), interface12_vertices[1]]
-# }
-# subdomain2_outer_boundary_verts = {
-#     0: None
-# }
-
-# list of subdomains given by the boundary polygon vertices.
-# Subdomains are given as a list of dolfin points forming
-# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
-# to create the subdomain. subdomain_def_points[0] contains the
-# vertices of the global simulation domain and subdomain_def_points[i] contains the
-# vertices of the subdomain i.
-subdomain_def_points = [sub_domain0_vertices,#
-                      sub_domain1_vertices,#
-                      sub_domain2_vertices]
-# in the below list, index 0 corresponds to the 12 interface which has index 1
-interface_def_points = [interface12_vertices]
-
-outer_boundary_def_points = {
-    # subdomain number
-    1 : subdomain1_outer_boundary_verts,
-    2 : subdomain2_outer_boundary_verts
-}
-
-# adjacent_subdomains[i] contains the indices of the subdomains sharing the
-# interface i (i.e. given by interface_def_points[i]).
-adjacent_subdomains = [[1,2]]
-isRichards = {
-    1: True, #
-    2: True
-    }
-
-Tmax = 2
-dt1 = 0.1
-dt2 = dt1
-# the timestep is taken as a dictionary to be able to run different timesteps on
-# different subdomains.
-timestep_size = {
-    1: dt1,
-    2: dt2
-}
-
-viscosity = {#
-# subdom_num : viscosity
-    1 : {'wetting' :1}, #
-    2 : {'wetting' :1}
-}
-
-porosity = {#
-# subdom_num : porosity
-    1 : 1,#
-    2 : 1
-}
-
-L = {#
-# subdom_num : subdomain L for L-scheme
-    1 : {'wetting' :0.25},#
-    2 : {'wetting' :0.25}
-}
-
-lambda_param = {#
-# subdom_num : lambda parameter for the L-scheme
-    1 : {'wetting' :4},#
-    2 : {'wetting' :4}
-}
-
-## relative permeabilty functions on subdomain 1
-def rel_perm1(s):
-    # relative permeabilty on subdomain1
-    return s**2
-
-_rel_perm1 = ft.partial(rel_perm1)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1,#
-    'nonwetting': None
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2(s):
-    # relative permeabilty on subdomain2
-    return s**3
-_rel_perm2 = ft.partial(rel_perm2)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2,#
-    'nonwetting': None
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
-}
-
-def saturation(pressure, subdomain_index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-
-sat_pressure_relationship = {#
-    1: ft.partial(saturation, subdomain_index = 1),#
-    2: ft.partial(saturation, subdomain_index = 2)
-}
-
-source_expression = {
-    1: {'wetting': '4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )'},
-    2: {'wetting': '2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))'}
-}
-
-initial_condition = {
-    1: {'wetting': '-(x[0]*x[0] + x[1]*x[1])'},#
-    2: {'wetting': '-x[1]*x[1]'}
-}
-
-exact_solution = {
-    1: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'},#
-    2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}
-}
-
-dirichletBC = exact_solution
-
-# def saturation(pressure, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-#
-# sa
-
-mesh_resolution = 3
-
-# initialise LDD simulation class
-simulation = ldd.LDDsimulation()
-simulation.set_parameters(output_dir = "",#
-    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,#
-    timestep_size = timestep_size,#
-    sources = source_expression,#
-    initial_conditions = initial_condition,#
-    outer_dirichletBC = dirichletBC,#
-    )
-
-simulation.initialise()
-# simulation._init_meshes_and_markers(subdomain_def_points, mesh_resolution=2)
-# subdomain marker functions
-domain_marker = simulation.domain_marker
-mesh_subdomain = simulation.mesh_subdomain
-
-# mesh = mesh_subdomain[0]
-# interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
-# interface_marker.set_all(0)
-# interface = dp.Interface(vertices=interface_def_points[0], #
-#                                     tol=mesh.hmin()/100,#
-#                                     internal=True,#
-#                                     adjacent_subdomains = adjacent_subdomains[0])
-# print("interface vertices are: \n", interface.vertices)
-# interface.mark(interface_marker, 1)
-# simulation._init_interfaces(interface_def_points, adjacent_subdomains)
-
-interface = simulation.interface
-interface_marker = simulation.interface_marker
-
-boundary_marker1 = simulation.subdomain[1].outer_boundary_marker
-boundary_marker2 = simulation.subdomain[2].outer_boundary_marker
-dx_1 = simulation.subdomain[1].dx
-dx_2 = simulation.subdomain[2].dx
-ds_1 = simulation.subdomain[1].ds
-ds_2 = simulation.subdomain[2].ds
-
-interface_marker1 = simulation.subdomain[1].interface_marker#df.MeshFunction('size_t', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
-interface_marker2 = simulation.subdomain[2].interface_marker#df.MeshFunction('size_t', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
-# interface_marker1.set_all(0)
-# interface_marker2.set_all(0)
-# print(dir(interface_marker1))
-
-interface[0].mark(interface_marker1, 1)
-interface[0].mark(interface_marker2, 1)
-# interface_marker1 = simulation.subdomain[1].boundary_marker
-
-# domain measures
-
-# dx1 = simulation.subdomain[1].dx
-dx1 = df.Measure('dx', domain = mesh_subdomain[1])
-dx2 = df.Measure('dx', domain = mesh_subdomain[2])
-
-# boundary measures
-ds1 = df.Measure('ds', domain = mesh_subdomain[1], subdomain_data = interface_marker1)
-ds2 = df.Measure('ds', domain = mesh_subdomain[2], subdomain_data = interface_marker2)
-
-V1 = df.FunctionSpace(mesh_subdomain[1], 'P', 1)
-p1_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2, t = 0)
-
-coordinates = mesh_subdomain[0].coordinates()
-print("global mesh coordinates of domain: \n", coordinates)
-coordinates1 = mesh_subdomain[1].coordinates()
-
-coordinates2 = mesh_subdomain[2].coordinates()
-submesh1_data = mesh_subdomain[1].data()
-dom1_parent_vertex_indices = submesh1_data.array('parent_vertex_indices',0)
-# print("map of dom1 vertices to parent vertices on global domain:\n", dom1_parent_vertex_indices)
-#
-# coordinates2 = mesh_subdomain[2].coordinates()
-# print("\nvertices of subdomain2: \n", coordinates2)
-submesh2_data = mesh_subdomain[2].data()
-dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
-# print("'parent_vertex_indices',1 \n")
-#
-#
-# print("on domain: \n")
-# interface_coordinates = interface[0].coordinates(interface_marker, 1, print_coordinates = True)
-# interface_coordinates = interface[0]._vertex_indices(interface_marker, 1, print_vertex_indices = True)
-# print("\non subdomain2: \n")
-# dom2_interface_coordinates = interface[0].coordinates(interface_marker2, 1, print_coordinates = True)
-dom2_interface_def_points = interface[0]._vertex_indices(interface_marker2, 1, print_vertex_indices = True)
-
-
-dofs_on_interface1 = interface[0].dofs_on_interface(V1, interface_marker1, 1)
-print("Dofs on Interface dom1", dofs_on_interface1)
-
-V2 = df.FunctionSpace(mesh_subdomain[2], 'P', 1)
-dofs_on_interface2 = interface[0].dofs_on_interface(V2, interface_marker2, 1)
-# print("Dofs on Interface dom1", dofs_on_interface2)
-
-# markermesh = interface_marker2.mesh()
-# functionSpaceMesh = V2.mesh()
-# print("Domain 2 mesh extracted from markerfunction:\n", markermesh.coordinates())
-# print("Domain 2 mesh extracted from Function Space V2:\n", functionSpaceMesh.coordinates())
-
-dom1_d2v = df.dof_to_vertex_map(V1)
-dom1_v2d = df.vertex_to_dof_map(V1)
-# print("dof to vertex map for V1:\n", dom1_d2v)
-# print("vertex to dof map for V1:\n", dom1_v2d)
-dom2_d2v = df.dof_to_vertex_map(V2)
-dom2_v2d = df.vertex_to_dof_map(V2)
-# print("dof to vertex map for V2:\n", dom2_d2v)
-# print("vertex to dof map for V2:\n", dom2_v2d)
-
-testf1 = df.Constant(41)
-f1test = df.interpolate(testf1, V1)
-testf2 = df.Constant(42)
-f2test = df.interpolate(testf2, V2)
-
-f1vec = f1test.compute_vertex_values()
-for i in range(len(f1vec)):
-    f1vec[i] = i
-
-# f1test.vector()[:] = f1vec
-f1nodal = f1test.vector().get_local()
-# for i in range(len(f1nodal)):
-#     f1nodal[i] = i
-
-# print("nodal values of f1nodal\n", f1nodal)
-
-#
-# print('dofs on interface \n', f1vec[dofs_on_interface])
-# f1test.vector()[:] = f1vec
-# print("f1test.vector()[:]\n", f1test.vector()[:])
-# print("f1test.vector()[dofs_on_interface]\n", f1test.vector()[dofs_on_interface])
-
-# print("constant function on dom1\n", f1vec)
-#
-# for i in dofs_on_interface1:
-#     print("local dof ind %i, parent = %i", i,dom1_parent_vertex_indices[dom1_d2v[i]])
-#     # f1test.vector() is enumerated by dof not by vertex
-#     f1test.vector()[i] = i
-# # this enumerates according to the dof numbering
-# print("testfunction f1test with set interface values\n", f1test.vector()[:])
-# print("testfunction f1test with get_local\n", f1test.vector().get_local())
-# f1test.vector()[dom1_v2d[dom1_d2v[dofs_on_interface]]] = f1vec[dofs_on_interface]
-# .vector()[:] and .vector().get_local() order values according to the vertex numbering
-# print("testfunction f1test with compute_vertex_values\n", f1test.compute_vertex_values())
-# print("testfunction f1test with compute_vertex_values mit d2vmap\n", f1test.compute_vertex_values()[dom1_d2v[:]])
-# i=0
-# for x in coordinates1:
-#      print(f"vertex {i}: f1test[{dom1_v2d[i]}] = {f1test.vector()[dom1_v2d[i]]}\tu({x}) = {f1test(x)}")
-#      i += 1
-# print("\non subdomain1: \n")
-print("vertices of subdomain1: \n", coordinates1)
-# dom1_interface_coordinates = interface[0].coordinates(interface_marker1, 1, print_coordinates = True)
-dom1_interface_def_points = interface[0]._vertex_indices(interface_marker1, 1, print_vertex_indices = True)
-
-interface[0].write_pressure_dofs(from_function = f1test, #
-                        interface_dofs = dofs_on_interface1,#
-                        dof_to_vert_map = dom1_d2v,#
-                        local_to_parent_vertex_map = dom1_parent_vertex_indices,#
-                        phase = 'wetting',
-                        subdomain_ind = 1)
-
-interface[0].read_pressure_dofs(to_function = f2test, #
-                        interface_dofs = dofs_on_interface2,#
-                        dof_to_vert_map = dom2_d2v,#
-                        local_to_parent_vertex_map = dom2_parent_vertex_indices,#
-                        phase = 'wetting',
-                        subdomain_ind = 2)
-
-i=0
-for x in coordinates2:
-     print(f"vertex {i}: f2test[{dom2_v2d[i]}] = {f2test.vector()[dom2_v2d[i]]}\tu({x}) = {f2test(x)}")
-     i += 1
-print("\non subdomain2: \n")
-print("vertices of subdomain2: \n", coordinates2)
-
-# dom2_interface_coordinates = interface[0].coordinates(interface_marker2, 1, print_coordinates = True)
-dom2_interface_def_points = interface[0]._vertex_indices(interface_marker2, 1, print_vertex_indices = True)
-print("Parent indices", dom2_parent_vertex_indices)
-
-
-u1 = df.TrialFunction(V1)
-v1 = df.TestFunction(V1)
-
-u2 = df.TrialFunction(V2)
-v2 = df.TestFunction(V2)
-
-def relative_permeability(s, subdomain_index):
-    if subdomain_index == 1:
-        # relative permeabilty on subdomain1
-        return s**2
-    else:
-        # relative permeabilty on subdomain2
-        return s**3
-
-def saturation(pressure, subdomain_index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-
-
-# exact solution
-# cell = df.Cell("triangle", 2)
-# x1 = df.SpatialCoordinate(mesh_subdomain[1])
-# x2 = df.SpatialCoordinate(mesh_subdomain[2])
-
-p2_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])', domain = mesh_subdomain[2], degree = 2, t = 0)
-p1_e = df.interpolate(p1_exact, V1)
-p2_e = df.interpolate(p2_exact, V2)
-#
-# initial conditions
-p1_initial = df.Expression('-(x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2)
-p2_initial = df.Expression('-x[1]*x[1]', domain = mesh_subdomain[2], degree = 2)
-p1_0 = df.interpolate(p1_initial, V1)
-p2_0 = df.interpolate(p2_initial, V2)
-
-# Initial boundary value for the pressure on interface
-p_gamma1 = df.DirichletBC(V1, p1_exact, interface_marker1, 1)
-p_gamma2 = df.DirichletBC(V2, p2_exact, interface_marker2, 1)
-
-outerBC1 = df.DirichletBC(V1, p1_exact, boundary_marker1, 1)
-outerBC2 = df.DirichletBC(V2, p2_exact, boundary_marker2, 1)
-
-
-### source terms
-# source1 = df.Expression('4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )', domain = mesh_subdomain[1], degree = 2, t = 0)
-# source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
-source1 = df.Expression(source_expression[1]['wetting'], domain = mesh_subdomain[1], degree = 2, t = 0)
-source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
-
-source1h = df.interpolate(source1, V1)
-source2h = df.interpolate(source2, V2)
-
-L1 = 0.25
-L2 = L1
-timestep = 0.1
-lambda1 = 4
-lambda2 = lambda1
-# initial iteartion
-p1_i = p1_0
-p2_i = p2_0
-# initial g_i
-n1 = df.FacetNormal(mesh_subdomain[1])
-flux = -df.grad(p1_0)
-gl0 = (df.dot(flux, n1) - lambda1*p1_0)*v1*ds1(1)
-print('gl0 term assembliert', df.assemble(gl0).get_local())
-g1_i = -lambda1*p1_i
-g2_i = -lambda1*p2_i
-
-for iterations in range(1):
-    #   +  +
-    a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
-    rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
-    rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
-    extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
-    splitrhs = df.assemble(rhssplit1)
-    print("splitrhs: \n", splitrhs.get_local())
-    extrat_assembled = df.assemble(extratrem)
-    print("extratrem: \n", extrat_assembled.get_local())
-    added = splitrhs + extrat_assembled
-    print("both added together:\n", added.get_local())
-
-    A1 = df.assemble(a1)
-    dsform = (timestep*lambda1*u1*v1)*ds1(1)
-    dsform_vec = df.assemble(dsform)
-    b1 = df.assemble(rhs1)
-    #p_gamma1.apply(A1,b1)
-    # outerBC1.apply(A1,b1)
-    u  = df.Function(V1)
-    U1 = u.vector()
-    df.solve(A1,U1,b1)
-    # print("Form:\n", A1.array())
-    print("RHS:\n", b1.get_local())
-    #print("ds term:\n", dsform_vec.array())
-    # print('solution:\n', U1.get_local())
-    p1_i.vector()[:] = U1
-    interface[0].write_pressure_dofs(from_function = p1_i, #
-                            interface_dofs = dofs_on_interface1,#
-                            dof_to_vert_map = dom1_d2v,#
-                            local_to_parent_vertex_map = dom1_parent_vertex_indices,#
-                            phase = 'wetting',#
-                            subdomain_ind = 1)
-#
-# Save mesh to file
-df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
-df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
-df.File('./test_global_interface_marker.pvd') << interface_marker
-df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
-df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
-df.File('./test_domain_marker.pvd') << domain_marker
-df.File('./test_domain_layered_soil_solution.pvd') << u
-df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
-df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
-df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
-df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
diff --git a/Test_misc/taaaeeeeeeschd.py b/Test_misc/taaaeeeeeeschd.py
new file mode 100755
index 0000000000000000000000000000000000000000..3d07eb32324284f50817b26dac6654055684f2ee
--- /dev/null
+++ b/Test_misc/taaaeeeeeeschd.py
@@ -0,0 +1,45 @@
+#!/usr/bin/python3
+
+import dolfin as df
+
+mesh = df.UnitIntervalMesh(2)
+
+dx = df.dx
+ds = df.ds
+
+V = df.FunctionSpace(mesh, 'CG', 1)
+
+v = df.TestFunction(V)
+u = df.Function(V)
+
+# set k constant to 42
+# k = df.Constant(42)
+k = df.Expression('2*x[0]', domain = mesh, degree=1)
+k = df.interpolate(k, V)
+
+n = df.FacetNormal(mesh)
+print(k.vector().get_local())
+print(df.assemble(k*v*dx).get_local())
+print(df.assemble(df.dot(df.grad(k),n)*v*ds).get_local())
+
+
+# solve k*u'' = 0
+a = u.dx(0)*k*v.dx(0)*dx + u.dx(0)*k.dx(0)*v*dx
+
+print('Form: \n', df.assemble(a).get_local())
+bcL = df.DirichletBC(V,df.Constant(10.),"on_boundary && x[0] < 0.5")
+bcR = df.DirichletBC(V,df.Constant(1.),"on_boundary && x[0] > 0.5")
+bcs= [bcL,bcR]
+
+df.solve(a == 0,u,bcs=bcs)
+
+df.File('./taaeschd_u1.pvd') << u
+
+# change value of k at one boundary point
+k.vector()[0] =  1.e6
+print(k.vector().get_local())
+print(df.assemble(k*dx))
+print(df.assemble(k*ds))
+
+df.solve(a == 0,u,bcs=bcs)
+df.File('./taaeschd_u2.pvd') << u
diff --git a/layered_soil.py b/layered-soil-case/layered_soil.py
similarity index 100%
rename from layered_soil.py
rename to layered-soil-case/layered_soil.py