diff --git a/domainPatch.py b/domainPatch.py
index 1abfc25388a6af5b3429b207855f3ee23dd3097b..2f78134f41839280f6106b66dab9be8cb05a6f1e 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -176,7 +176,7 @@ class DomainPatch(df.SubDomain):
         """ save the interface values of self.pressure to all neighbouring interfaces"""
         if self.isRichards:
             for ind in self.has_interface:
-                interface[ind].write_dofs(from_function = self.pressure['wetting'], #
+                interface[ind].write_pressure_dofs(from_function = self.pressure['wetting'], #
                                         interface_dofs = self._dof_indices_of_interface[ind]['wetting'],#
                                         dof_to_vert_map = self.dof2vertex['wetting'],#
                                         local_to_parent_vertex_map = self.parent_mesh_index,#
@@ -185,7 +185,7 @@ class DomainPatch(df.SubDomain):
         else:
             for ind in self.has_interface:
                 for phase in ['wetting', 'nonwetting']:
-                    interface[ind].write_dofs(from_function = self.pressure[phase], #
+                    interface[ind].write_pressure_dofs(from_function = self.pressure[phase], #
                                             interface_dofs = self._dof_indices_of_interface[ind][phase],#
                                             dof_to_vert_map = self.dof2vertex[phase],#
                                             local_to_parent_vertex_map = self.parent_mesh_index,#
@@ -891,14 +891,12 @@ class Interface(BoundaryPart):
                     self.gli_term_prev[subdom_ind][phase].update({vert_ind: 0.0})
         # end init_interface_values
 
-    def write_dofs(self, from_function: df.Function, #
+    def write_pressure_dofs(self, from_function: df.Function, #
                    interface_dofs: np.array,#
                    dof_to_vert_map: np.array,#
                    local_to_parent_vertex_map: np.array,#
                    phase: str,#
                    subdomain_ind: int,#
-                   pressure: bool = False,#
-                   gl: bool = False,#
                    previous_iter: bool = False,#
                    ) -> None:
         """ write dofs of from_function with indices interface_dofs to self.pressure_values
@@ -908,11 +906,58 @@ class Interface(BoundaryPart):
         of the interface
             self.pressure_values[subdomain_ind][phase], (if pressure == True)
             self.pressure_values_prev[subdomain_ind][phase], (if previous_iter == True)
+
+        # Parameters
+        from_function:              #type: df.Function, function to save the dofs of
+        interface_dofs:             #type: np.array,    index array of dofs to be saved
+        dof_to_vert_map:            #type: np.array,    dof to vert map of the function
+                                                        space on the submesh.
+        local_to_parent_vertex_map: #type: np.array,
+        phase                       #type: str          is either 'wetting' or
+                                                        'nonwetting' and determins
+                                                        the dict to be written to.
+        subdomain_ind               #type: int          subdomain index
+        previous_iter               #type: bool         determins wether a previous
+                                                        iterate is being written or not
+        """
+        if self.pressure_values == None:
+            raise RuntimeError("self.pressure_values not initiated. You forgot to \
+                                run self.init_interface_values")
+
+        parent = local_to_parent_vertex_map
+        d2v = dof_to_vert_map
+        if not previous_iter:
+            for dof_index in interface_dofs:
+                # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                self.pressure_values[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
+            print("have written pressure interface values\n", self.pressure_values[subdomain_ind][phase])
+        else:
+            for dof_index in interface_dofs:
+                # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
+                # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                self.pressure_values_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
+            print("have written previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
+
+
+    def write_gli_dofs(self, from_array: df.Function, #
+                   interface_dofs: np.array,#
+                   dof_to_vert_map: np.array,#
+                   local_to_parent_vertex_map: np.array,#
+                   phase: str,#
+                   subdomain_ind: int,#
+                   previous_iter: bool = False,#
+                   ) -> None:
+        """ write dofs of from_array with indices interface_dofs to self.gli_term
+
+        Write the degrees of freedom of the array from_array with indices
+        interface_dofs, i.e. the dofs of the interface, to one of the dictionaries
+        of the interface
             self.gli_term[subdomain_ind][phase], (if gl == True)
             self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
 
         # Parameters
-        from_function:              #type: df.Function, function to save the dofs to
+        from_array:                 #type: np.array,    array to save the dofs of
         interface_dofs:             #type: np.array,    index array of dofs to be saved
         dof_to_vert_map:            #type: np.array,    dof to vert map of the function
                                                         space on the submesh.
@@ -921,10 +966,6 @@ class Interface(BoundaryPart):
                                                         'nonwetting' and determins
                                                         the dict to be written to.
         subdomain_ind               #type: int          subdomain index
-        pressure                    #type: bool         determins wether pressure
-                                                        values are being written or not
-        gl                          #type: bool         determins wether gl
-                                                        values are being written or not
         previous_iter               #type: bool         determins wether a previous
                                                         iterate is being written or not
         """
@@ -934,33 +975,19 @@ class Interface(BoundaryPart):
 
         parent = local_to_parent_vertex_map
         d2v = dof_to_vert_map
-        # chose where to write the dofs to
-        if pressure:
-            if not previous_iter:
-                for dof_index in interface_dofs:
-                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
-                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
-                    self.pressure_values[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
-                print("have written pressure interface values\n", self.pressure_values[subdomain_ind][phase])
-            else:
-                for dof_index in interface_dofs:
-                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
-                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
-                    self.pressure_values_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
-                print("have written previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
-        if gl:
-            if not previous_iter:
-                for dof_index in interface_dofs:
-                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
-                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
-                    self.gli_term[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
-                print("have written gl interface values\n", self.gli_term[subdomain_ind][phase])
-            else:
-                for dof_index in interface_dofs:
-                    # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
-                    # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
-                    self.gli_term_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
-                print("have written previous gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
+        if not previous_iter:
+            for dof_index in interface_dofs:
+                # from_array is orderd by dof and not by vertex numbering of the mesh.
+                # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                self.gli_term[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
+            print("have written gl interface values\n", self.gli_term[subdomain_ind][phase])
+        else:
+            for dof_index in interface_dofs:
+                # from_array is orderd by dof and not by vertex numbering of the mesh.
+                # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
+                self.gli_term_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
+            print("have written previous gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
+
 
     def read_pressure_dofs(self, to_function: df.Function, #
                   interface_dofs: np.array,#