diff --git a/domainPatch.py b/domainPatch.py
index 2f78134f41839280f6106b66dab9be8cb05a6f1e..3bb5286e61d32b44e8bba94774ef4958da87addd 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -332,6 +332,8 @@ class DomainPatch(df.SubDomain):
             neighbour_iter_num = interface.current_iteration[neighbour]
             dt = self.timestep_size
             ds = self.ds(ind)
+            # needed for the write_gli_dofs() functions
+            interface_dofs = self._dof_indices_of_interface[ind],#
 
             Lambda = self.lambda_param
             if iteration == 0:
@@ -365,8 +367,15 @@ class DomainPatch(df.SubDomain):
                     self.gli_assembled = {
                         'wetting': df.assemble(dt*gl0*phi_w*ds)
                         }
-                    ### TODO write dofs to interface dict
-                    #interface.gli_term[subdomain].update({'wetting': gl0})
+                    # write glo to the current iteration dictionary of the interface
+                    interface.write_gli_dofs(from_array = self.gli_assembled['wetting'], #
+                                   interface_dofs = interface_dofs['wetting'],#
+                                   dof_to_vert_map = self.dof2vertex['wetting'],#
+                                   local_to_parent_vertex_map = self.parent_mesh_index,#
+                                   phase = 'wetting',#
+                                   subdomain_ind = subdomain,#
+                                   previous_iter = False
+                                   )
                 else:
                     phi_w  = self.testfunction['wetting']
                     mu_w = mu['wetting']
@@ -378,7 +387,7 @@ class DomainPatch(df.SubDomain):
                     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
+                        # nonwetting phase and assemble the terms as zero form.
                         fluxw = -1/mu_w*kw(S(pc_prev))*df.grad(pw_prev)
                         gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
                         self.gli_assembled = {
@@ -386,7 +395,6 @@ class DomainPatch(df.SubDomain):
                             # 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']
@@ -401,10 +409,25 @@ class DomainPatch(df.SubDomain):
                             '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})
-
+                    # write glo to the current iteration dictionary of the interface
+                    # wetting
+                    interface.write_gli_dofs(from_array = self.gli_assembled['wetting'], #
+                                   interface_dofs = interface_dofs['wetting'],#
+                                   dof_to_vert_map = self.dof2vertex['wetting'],#
+                                   local_to_parent_vertex_map = self.parent_mesh_index,#
+                                   phase = 'wetting',#
+                                   subdomain_ind = subdomain,#
+                                   previous_iter = False
+                                   )
+                    # nonwetting
+                    interface.write_gli_dofs(from_array = self.gli_assembled['nonwetting'], #
+                                   interface_dofs = interface_dofs['nonwetting'],#
+                                   dof_to_vert_map = self.dof2vertex['nonwetting'],#
+                                   local_to_parent_vertex_map = self.parent_mesh_index,#
+                                   phase = 'nonwetting',#
+                                   subdomain_ind = subdomain,#
+                                   previous_iter = False
+                                   )
                 interface.current_iteration[subdomain] += 1
 
             # else: # iteration > 0
@@ -940,7 +963,7 @@ class Interface(BoundaryPart):
             print("have written previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
 
 
-    def write_gli_dofs(self, from_array: df.Function, #
+    def write_gli_dofs(self, from_array: np.array, #
                    interface_dofs: np.array,#
                    dof_to_vert_map: np.array,#
                    local_to_parent_vertex_map: np.array,#