diff --git a/LDDsimulation.py b/LDDsimulation.py
index bc3064f21b706e0dcdf327c51d1b8fa6ca5716dd..f137716346bf1d2142d4cedcb8ff3ca2076b1735 100644
--- a/LDDsimulation.py
+++ b/LDDsimulation.py
@@ -93,7 +93,7 @@ class LDDsimulation(object):
                        output_dir: str,#
                        subdomain_def_points: tp.List[tp.List[df.Point]],#
                        # this is actually an array of bools
-                       isRichards: bool,#
+                       isRichards: tp.Dict[int, bool],#
                        interface_def_points: tp.List[tp.List[df.Point]],#
                        outer_boundary_def_points: tp.Dict[int, tp.Dict[int, tp.List[df.Point]]],#
                        adjacent_subdomains: tp.List[np.ndarray],#
@@ -267,11 +267,10 @@ class LDDsimulation(object):
     def _init_subdomains(self) -> None:
         """ generate list of opjects of class DomainPatch.
         """
-        is_richards = self.isRichards
 
         # The list is_richards has the same length as the nuber of subdomains.
         # Therefor it is used in the subdomain initialisation loop.
-        for subdom_num, isR in enumerate(is_richards, 1):
+        for subdom_num, isR in self.isRichards.items():
             interface_list = self._get_interfaces(subdom_num)
             self.subdomain.update(#
                     {subdom_num : dp.DomainPatch(#
diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
index 8fbb509932c075a4af6afb1f2804466fbec14865..165394a3bf542373cad8496bf8367d13b6f84286 100755
--- a/RR-2-patch-test.py
+++ b/RR-2-patch-test.py
@@ -76,7 +76,10 @@ outer_boundary_def_points = {
 # 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 = [True, True]
+isRichards = {
+    1: True, #
+    2: True
+    }
 
 Tmax = 2
 dt1 = 0.1
@@ -172,6 +175,7 @@ dirichletBC = exact_solution
 #
 # sa
 
+mesh_resolution = 3
 
 # initialise LDD simulation class
 simulation = ldd.LDDsimulation()
@@ -181,7 +185,7 @@ simulation.set_parameters(output_dir = "",#
     interface_def_points = interface_def_points,#
     outer_boundary_def_points = outer_boundary_def_points,#
     adjacent_subdomains = adjacent_subdomains,#
-    mesh_resolution = 2,#
+    mesh_resolution = mesh_resolution,#
     viscosity = viscosity,#
     porosity = porosity,#
     L = L,#
@@ -251,7 +255,7 @@ 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)
+# 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)
@@ -263,7 +267,7 @@ dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
 # 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")
+# 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)
 
@@ -273,7 +277,7 @@ 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)
+# print("Dofs on Interface dom1", dofs_on_interface2)
 
 # markermesh = interface_marker2.mesh()
 # functionSpaceMesh = V2.mesh()
@@ -282,12 +286,12 @@ print("Dofs on Interface dom1", dofs_on_interface2)
 
 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)
+# 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)
+# 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)
@@ -303,7 +307,7 @@ f1nodal = f1test.vector().get_local()
 # for i in range(len(f1nodal)):
 #     f1nodal[i] = i
 
-print("nodal values of f1nodal\n", f1nodal)
+# print("nodal values of f1nodal\n", f1nodal)
 
 #
 # print('dofs on interface \n', f1vec[dofs_on_interface])
@@ -311,24 +315,24 @@ print("nodal values of f1nodal\n", f1nodal)
 # 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())
+# 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("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)
@@ -394,8 +398,8 @@ 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, df.Constant(0.5), interface_marker1, 1)
-p_gamma2 = df.DirichletBC(V2, df.Constant(0.1), interface_marker2, 1)
+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)
@@ -419,26 +423,39 @@ lambda2 = lambda1
 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):
-    #   + relative_permeability(saturation(p1_i, 1), 1)* +
-    a1 = (df.inner(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
+    #   +  +
+    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)
+    #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("Form:\n", A1.array())
     print("RHS:\n", b1.get_local())
-    print("ds term:\n", dsform_vec.array())
-    print('solution:\n', U1.get_local())
+    #print("ds term:\n", dsform_vec.array())
+    # print('solution:\n', U1.get_local())
     p1_i.vector()[:] = U1
     interface[0].write_dofs(from_function = p1_i, #
                             interface_dofs = dofs_on_interface1,#
diff --git a/domainPatch.py b/domainPatch.py
index 44b7964cb3602dc87dfcab9d3cb3d09bc9fb7f90..de36a41f76bcacedebc044bc4bf140dd99f8398d 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -150,15 +150,19 @@ class DomainPatch(df.SubDomain):
         # solution function(s) for the form set by _init_function_space
         self.pressure = None
         self.source = None
-        # dictionary holding the previous iteration.
+        # FEM function holding the previous iteration.
         self.pressure_prev_iter = None
-        # dictionary holding the pressures of the previous timestep t_n.
+        # FEM function holding the pressures of the previous timestep t_n.
         self.pressure_prev_timestep = None
         # test function(s) for the form set by _init_function_space
-        self.TestFunction = None
-        # # dictionary which is initiated by _init_function_space holding the
-        # # the gli term needed for the form assemply.
-        # self.gli = None
+        self.testfunction = None
+        # Dictionary of Arrays (one for 'wetting' and one for 'nonwetting')
+        # corresponding to the dofs of an FEM function which is initiated
+        # by calc_gli_term() holding the dofs of the assembled dt*gli*v*ds terms as
+        # a sparse vector (here, dt = self.timestep_size). This vector is read by the solver and added to the rhs
+        # which is given by self.govering_problem().
+        #
+        self.gli_assembled = None
         self._init_function_space()
         self._init_dof_and_vertex_maps()
         self._init_markers()
@@ -194,7 +198,7 @@ class DomainPatch(df.SubDomain):
             iter_num as a dictionary.
 
         return the governing form ant right hand side for phase phase and iteration
-        iter_num (without the gli terms) as a dictionary. 
+        iter_num (without the gli terms) as a dictionary.
         """
         # define measures
         dx = self.dx
@@ -227,7 +231,7 @@ class DomainPatch(df.SubDomain):
             # # assemble rhs for
             # rhs_interface_forms = []
             # self._calc_gli_term(iteration = iter_num)
-            # # gli = self.gli['wetting']
+            # # gli = self.gli_assembled['wetting']
             # for index in self.has_interface:
             #     gli = self.interface[index].gli_term[self.subdomain_index]['wetting']
             #     rhs_interface_forms.append((dt*gli*phi_w)*ds[index])
@@ -241,10 +245,7 @@ class DomainPatch(df.SubDomain):
             return form_and_rhs
         else:
             La = self.L['alpha']
-            # Lnw = self.L['nonwetting']
-            # this is pw_i/pnw_i in the paper
             pa = self.pressure['alpha']
-            # pnw = self.pressure['nonwetting']
             # this is pw_{i-1} in the paper
             pw_prev_iter = self.pressure_prev_iter['wetting']
             pnw_prev_iter = self.pressure_prev_iter['nonwetting']
@@ -253,17 +254,14 @@ class DomainPatch(df.SubDomain):
             pc_prev_iter = pnw_prev_iter - pw_prev_iter
             pc_prev_timestep = pnw_prev_timestep - pw_prev_timestep
             phi_a  = self.testfunction['alpha']
-            # phi_nw  = self.testfunction['nonwetting']
             Lambda_a = self.lambda_param['alpha']
-            # Lambda_nw = self.lambda_param['nonwetting']
             ka = self.relative_permeability['alpha']
-            # knw = self.relative_permeability['nonwetting']
             S = self.saturation
             mu_a = self.viscosity['alpha']
-            # mu_nw = self.viscosity['nonwetting']
             source_a = self.source['alpha']
-            # source_nw = self.source['nonwetting']
-
+            # the forms of wetting and nonwetting phase differ by a minus sign
+            # and by interface_form terms if the neighbour of the current patch
+            # assumes Richards model
             if phase == "wetting":
                 # gather interface forms
                 interface_forms = []
@@ -316,6 +314,132 @@ class DomainPatch(df.SubDomain):
                 'Expected either phase == wetting or phase == nonwetting ')
                 return None
 
+    def calc_gli_term(self, iteration: int) -> None: # tp.Dict[str, fl.Form]
+        """calculate the gl terms for the iteration number iteration
+
+        calculate the gl terms for the iteration number iteration and save them
+        to self.gli_assembled. The term gets saved as an already assembled sparse
+        ds form to be added to the rhs by the solver. This is due to the need of
+        having to communicate dofs over the interface. Additionally, the dofs are
+        being saved to an interface dictionary exactly as the pressures.
+        """
+        for ind in self.has_interface:
+            # self._calc_gli_term_on(interface, iteration)
+            interface = self.interface[ind]
+            subdomain = self.subdomain_index
+            # neighbour index
+            neighbour = interface.neighbour[subdomain]
+            neighbour_iter_num = interface.current_iteration[neighbour]
+            dt = self.timestep_size
+            ds = self.ds(ind)
+
+            Lambda = self.lambda_param
+            if iteration == 0:
+                n = df.FacetNormal(self.mesh)
+                pw_prev_iter = self.pressure_prev_timestep
+                k = self.relative_permeability
+                S = self.saturation
+                mu = self.viscosity
+
+                if not neighbour_iter_num == 0:
+                    # save the gl term from the neighbour first to use it in the
+                    # next iteration
+                    if self.isRichards:
+                        interface.gli_term_prev[subdomain].update(#
+                            {'wetting': interface.gli_term[neighbour]['wetting']}
+                            )
+                    else:
+                        interface.gli_term_prev[subdomain].update(#
+                            {'wetting': interface.gli_term[neighbour]['wetting']},#
+                            {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
+                            )
+
+                if self.isRichards:
+                    mu_w = mu['wetting']
+                    kw = k['wetting']
+                    pw_prev = pw_prev_iter['wetting']
+                    phi_w  = self.testfunction['wetting']
+                    #calc gl0 from the flux term
+                    flux = -1/mu_w*kw(S(pw_prev))*df.grad(pw_prev)
+                    gl0 = (df.dot(flux, n) - Lambda['wetting']*pw_prev)
+                    self.gli_assembled = {
+                        'wetting': df.assemble(dt*gl0*phi_w*ds)
+                        }
+                    ### TODO write dofs to interface dict
+                    #interface.gli_term[subdomain].update({'wetting': gl0})
+                else:
+                    phi_w  = self.testfunction['wetting']
+                    mu_w = mu['wetting']
+                    pw_prev = pw_prev_iter['wetting']
+                    pnw_prev = pw_prev_iter['nonwetting']
+                    pc_prev = pnw_prev - pw_prev
+                    kw = k['wetting']
+                    #calc gl0 from the flux term
+                    if interface.neighbour_isRichards[subdomain]:
+                        # if the neighbour of our subdomain (with index self.subdomain_index)
+                        # assumes a Richards model, we don't have gli terms for the
+                        # nonwetting phase and
+                        fluxw = -1/mu_w*kw(S(pc_prev))*df.grad(pw_prev)
+                        gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
+                        self.gli_assembled = {
+                            'wetting': df.assemble(dt*gwl0*phi_w*ds),
+                            # assemble the nonwetting phase gl0 term as zero functional
+                            'nonwetting': df.assemble(df.Constant(0)*ds)
+                            }
+                            ###TODO write gli_terms to interface
+                    else:
+                        # here we need to assemble a gli term for the nonwetting phase
+                        phi_nw  = self.testfunction['nonwetting']
+                        mu_nw = mu['nonwetting']
+                        knw = k['nonwetting']
+
+                        fluxw = -1/mu_w*kw(S(pc_prev))*df.grad(pw_prev)
+                        fluxnw = -1/mu_nw*knw(1-S(pc_prev))*df.grad(pnw_prev)
+                        gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
+                        gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
+                        self.gli_assembled = {
+                            'wetting': df.assemble(dt*gwl0*phi_w*ds),
+                            'nonwetting': df.assemble(dt*gnwl0*phi_nw*ds)
+                            }
+                        ###TODO write gli_terms to interface
+                        # interface.gli_term[subdomain].update({'wetting': gwl0})
+                        # interface.gli_term[subdomain].update({'nonwetting': gnwl0})
+
+                interface.current_iteration[subdomain] += 1
+
+            # else: # iteration > 0
+            #     if neighbour_iter_num == iteration:
+            #         if self.isRichards:
+            #             interface.gli_term_prev[subdomain].update(#
+            #                 {'wetting': interface.gli_term[neighbour]['wetting']}
+            #                 )
+            #
+            #             pw_prev = pw_prev_iter['wetting']
+            #             #calc gl0 from the flux term
+            #             flux = -1/mu_w*kw(S(pw_prev))*df.grad(pw_prev)
+            #             gl0 = (df.dot(flux, n) - Lambda['wetting']*pw_prev)
+            #             interface.gli_term[subdomain].update({'wetting': gl0})
+            #         else:
+            #             mu_w = mu['wetting']
+            #             mu_nw = mu['nonwetting']
+            #             pw_prev = pw_prev_iter['wetting']
+            #             pnw_prev = pw_prev_iter['nonwetting']
+            #             kw = k['wetting']
+            #             knw = k['nonwetting']
+            #             #calc gl0 from the flux term
+            #             fluxw = -1/mu_w*kw(S(pnw_prev - pw_prev))*df.grad(pw_prev)
+            #             fluxnw = -1/mu_nw*knw(1-S(pnw_prev - pw_prev))*df.grad(pnw_prev)
+            #             gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
+            #             gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
+            #             interface.gli_term[subdomain].update({'wetting': gwl0})
+            #             interface.gli_term[subdomain].update({'nonwetting': gnwl0})
+            #
+            #     else:
+            #         print('uiae')
+            #
+            #     interface.current_iteration[subdomain] += 1
+
+
     #### PRIVATE METHODS
     def _init_function_space(self) -> None:
         """ create function space for solution and trial functions
@@ -324,13 +448,13 @@ class DomainPatch(df.SubDomain):
         aswell. This method sets the class variable
             self.function_space
             self.pressure
-            self.TestFunction
+            self.testfunction
         """
         if self.isRichards:
             self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
             self.pressure = df.TrialFunction(self.function_space['wetting'])
             self.testfunction = df.TestFunction(self.function_space['wetting'])
-            # self.gli = {'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])}
+            # self.gli_assembled = {'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])}
         else:
             self.function_space = {#
                 'wetting'    : df.FunctionSpace(self.mesh, 'P', 1),#
@@ -344,7 +468,7 @@ class DomainPatch(df.SubDomain):
                 'wetting'    : df.TestFunction(self.function_space['wetting']),
                 'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
                 }
-            # self.gli = {#
+            # self.gli_assembled = {#
             #     'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting']),#
             #     'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
             # }
@@ -438,97 +562,6 @@ class DomainPatch(df.SubDomain):
                      'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, ind)}
                 )
 
-    def _calc_gli_term(self, iteration: int) -> None: # tp.Dict[str, fl.Form]
-        """calculate the gl terms for the iteration number iteration
-
-        calculate the gl terms for the iteration number iteration and save them
-        to self.gli.
-        """
-        for ind in self.has_interface:
-            # self._calc_gli_term_on(interface, iteration)
-            interface = self.interface[ind]
-            subdomain = self.subdomain_index
-            # neighbour index
-            neighbour = interface.neighbour[subdomain]
-            neighbour_iter_num = interface.current_iteration[neighbour]
-
-            Lambda = self.lambda_param
-            if iteration == 0:
-                n = df.FacetNormal(self.mesh)
-                pw_prev_iter = self.pressure_prev_timestep
-                k = self.relative_permeability
-                S = self.saturation
-                mu = self.viscosity
-
-                if not neighbour_iter_num == 0:
-                    # save the gl term from the neighbour first to use it in the
-                    # next iteration
-                    if self.isRichards:
-                        interface.gli_term_prev[subdomain].update(#
-                            {'wetting': interface.gli_term[neighbour]['wetting']}
-                            )
-                    else:
-                        interface.gli_term_prev[subdomain].update(#
-                            {'wetting': interface.gli_term[neighbour]['wetting']},#
-                            {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
-                            )
-
-                if self.isRichards:
-                    mu_w = mu['wetting']
-                    kw = k['wetting']
-                    pwn_prev = pw_prev_iter['wetting']
-                    #calc gl0 from the flux term
-                    flux = -1/mu_w*kw(S(pwn_prev))*df.grad(pwn_prev)
-                    gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
-                    interface.gli_term[subdomain].update({'wetting': gl0})
-                else:
-                    mu_w = mu['wetting']
-                    mu_nw = mu['nonwetting']
-                    pwn_prev = pw_prev_iter['wetting']
-                    pnwn_prev = pw_prev_iter['nonwetting']
-                    kw = k['wetting']
-                    knw = k['nonwetting']
-                    #calc gl0 from the flux term
-                    fluxw = -1/mu_w*kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
-                    fluxnw = -1/mu_nw*knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
-                    gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
-                    gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
-                    interface.gli_term[subdomain].update({'wetting': gwl0})
-                    interface.gli_term[subdomain].update({'nonwetting': gnwl0})
-
-                interface.current_iteration[subdomain] += 1
-
-            # else: # iteration > 0
-            #     if neighbour_iter_num == iteration:
-            #         if self.isRichards:
-            #             interface.gli_term_prev[subdomain].update(#
-            #                 {'wetting': interface.gli_term[neighbour]['wetting']}
-            #                 )
-            #
-            #             pwn_prev = pw_prev_iter['wetting']
-            #             #calc gl0 from the flux term
-            #             flux = -1/mu_w*kw(S(pwn_prev))*df.grad(pwn_prev)
-            #             gl0 = (df.dot(flux, n) - Lambda['wetting']*pwn_prev)
-            #             interface.gli_term[subdomain].update({'wetting': gl0})
-            #         else:
-            #             mu_w = mu['wetting']
-            #             mu_nw = mu['nonwetting']
-            #             pwn_prev = pw_prev_iter['wetting']
-            #             pnwn_prev = pw_prev_iter['nonwetting']
-            #             kw = k['wetting']
-            #             knw = k['nonwetting']
-            #             #calc gl0 from the flux term
-            #             fluxw = -1/mu_w*kw(S(pnwn_prev - pwn_prev))*df.grad(pwn_prev)
-            #             fluxnw = -1/mu_nw*knw(1-S(pnwn_prev - pwn_prev))*df.grad(pnwn_prev)
-            #             gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pwn_prev)
-            #             gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnwn_prev)
-            #             interface.gli_term[subdomain].update({'wetting': gwl0})
-            #             interface.gli_term[subdomain].update({'nonwetting': gnwl0})
-            #
-            #     else:
-            #         print('uiae')
-            #
-            #     interface.current_iteration[subdomain] += 1
 
 # END OF CLASS DomainPatch