diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index ac93c5e5aebcea9a3bf1ebfe2fa87ac936abacf6..bb9974426c93073626115d9523a1515ac2f6c449 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -142,6 +142,10 @@ 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()
+        # dictionary holding for each interface in self.has_interface the dof
+        # indices that are shared with another interface (possibly none). This info is needed for
+        # the calculation of the gli terms.
+        self._interface_has_common_dof_indices = dict()
         # List of objects of clas outerBoundary initialised by self._init_markers()
         # can be NONE if no outer boundary is present.
         self.outer_boundary = []
@@ -175,6 +179,7 @@ class DomainPatch(df.SubDomain):
         self._init_markers()
         self._init_measures()
         self._calc_dof_indices_of_interfaces()
+        self._calc_interface_has_common_dof_indices()
 
     # END constructor
 
@@ -198,22 +203,18 @@ class DomainPatch(df.SubDomain):
             from_func = self.pressure
 
         if self.isRichards:
-            for ind in self.has_interface:
-                interface[ind].read_pressure_dofs(from_function = from_func['wetting'], #
-                                        interface_dofs = self._dof_indices_of_interface[ind]['wetting'],#
-                                        dof_to_vert_map = self.dof2vertex['wetting'],#
+            subdomain_has_phases = ['wetting']
+        else:
+            subdomain_has_phases = ['wetting', 'nonwetting']
+
+        for ind in self.has_interface:
+            for phase in subdomain_has_phases:
+                self.interface[ind].read_pressure_dofs(from_function = from_func[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,#
-                                        phase = 'wetting',#
+                                        phase = phase,#
                                         subdomain_ind = self.subdomain_index)
-        else:
-            for ind in self.has_interface:
-                for phase in ['wetting', 'nonwetting']:
-                    interface[ind].read_pressure_dofs(from_function = from_func[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,#
-                                            phase = phase,#
-                                            subdomain_ind = self.subdomain_index)
 
 
     def govering_problem(self, phase: str) -> tp.Dict[str, fl.Form]:
@@ -248,8 +249,8 @@ class DomainPatch(df.SubDomain):
             # we need to have all interfaces in the form
             interface_forms = []
             for interface in self.has_interface:
-                interface_forms.append((dt*Lambda*pw*phi_w)*ds(interface)
-            form1 = (Lw*pw*phi_w)*dx
+                interface_forms.append((dt*Lambda*pw*phi_w)*ds(interface))
+            form1 = Lw*pw*phi_w*dx
             form2 = dt/mu_w*df.dot(kw(S(pw_prev_iter))*df.grad(pw), df.grad(v))*dx
             form = form1 + form2 + df.sum(interface_forms)
             # # assemble rhs
@@ -348,7 +349,7 @@ class DomainPatch(df.SubDomain):
         subdomain = self.subdomain_index
         dt = self.timestep_size
         Lambda = self.lambda_param
-        
+
         for ind in self.has_interface:
             # self._calc_gli_term_on(interface, iteration)
             interface = self.interface[ind]
@@ -358,6 +359,7 @@ class DomainPatch(df.SubDomain):
             ds = self.ds(ind)
             # needed for the read_gli_dofs() functions
             interface_dofs = self._dof_indices_of_interface[ind],#
+
             if iteration == 0:
                 n = df.FacetNormal(self.mesh)
                 pw_prev_iter = self.pressure_prev_timestep
@@ -710,6 +712,12 @@ class DomainPatch(df.SubDomain):
                      'nonwetting': self.interface[ind].dofs_on_interface(V['nonwetting'], marker, ind)}
                 )
 
+    def _calc_interface_has_common_dof_indices(self):
+        """ populate dictionary self._interface_has_common_dof_indices
+
+        Determin for each interface in self.has_interface the dof indices
+        """
+
 
 # END OF CLASS DomainPatch