diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index ad5411327c16e7ae83e4ae230a714d62d1d875d6..1ccbbb1a5e814cd908ee2950d604d85713993375 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -627,37 +627,51 @@ class DomainPatch(df.SubDomain):
         Note that P1 FEM is hard coded here, as it is assumed in other methods
         aswell. This method sets the class variable
             self.function_space
+            self.trialfunction
             self.pressure
             self.testfunction
         """
-        if self.isRichards:
-            self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
-            self.trialfunction = {
-                'wetting': df.TrialFunction(self.function_space['wetting'])
-                }
-            self.testfunction = {
-                'wetting': df.TestFunction(self.function_space['wetting'])
-                }
-            self.neighbouring_interface_pressure = {
-                'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])
-                }
-        else:
-            self.function_space = {#
-                'wetting'    : df.FunctionSpace(self.mesh, 'P', 1),#
-                'nonwetting' : df.FunctionSpace(self.mesh, 'P', 1)#
-                }
-            self.trialfunction = {
-                'wetting'    : df.TrialFunction(self.function_space['wetting']),
-                'nonwetting' : df.TrialFunction(self.function_space['nonwetting'])
-                }
-            self.testfunction = {
-                'wetting'    : df.TestFunction(self.function_space['wetting']),
-                'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
-                }
-            self.neighbouring_interface_pressure = {#
-                'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting']),#
-                'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
-            }
+        self.function_space = dict()
+        self.trialfunction = dict()
+        self.testfunction = dict()
+        self.neighbouring_interface_pressure = dict()
+        for phase in self.has_phases:
+            self.function_space.update({phase: df.FunctionSpace(self.mesh, 'P', 1)})
+            self.trialfunction.update({phase: df.TrialFunction(self.function_space[phase])})
+            self.testfunction.update({phase: df.TestFunction(self.function_space[phase])})
+            self.neighbouring_interface_pressure.update(
+                {phase: df.interpolate(df.Constant(0), self.function_space[phase])}
+                )
+            print(f"Python id of self.neighbouring_interface_pressure[{phase}]", id(self.neighbouring_interface_pressure[phase]))
+
+        # if self.isRichards:
+        #     self.function_space = {'wetting': df.FunctionSpace(self.mesh, 'P', 1)}
+        #     self.trialfunction = {
+        #         'wetting': df.TrialFunction(self.function_space['wetting'])
+        #         }
+        #     self.testfunction = {
+        #         'wetting': df.TestFunction(self.function_space['wetting'])
+        #         }
+        #     self.neighbouring_interface_pressure = {
+        #         'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting'])
+        #         }
+        # else:
+        #     self.function_space = {#
+        #         'wetting'    : df.FunctionSpace(self.mesh, 'P', 1),#
+        #         'nonwetting' : df.FunctionSpace(self.mesh, 'P', 1)#
+        #         }
+        #     self.trialfunction = {
+        #         'wetting'    : df.TrialFunction(self.function_space['wetting']),
+        #         'nonwetting' : df.TrialFunction(self.function_space['nonwetting'])
+        #         }
+        #     self.testfunction = {
+        #         'wetting'    : df.TestFunction(self.function_space['wetting']),
+        #         'nonwetting' : df.TestFunction(self.function_space['nonwetting'])
+        #         }
+        #     self.neighbouring_interface_pressure = {#
+        #         'wetting' : df.interpolate(df.Constant(0), self.function_space['wetting']),#
+        #         'nonwetting' : df.interpolate(df.Constant(0), self.function_space['nonwetting'])
+        #     }
 
     def _init_dof_maps(self) -> None:
         """ calculate dof maps and the vertex to parent map.
@@ -745,9 +759,9 @@ class DomainPatch(df.SubDomain):
             )
             for phase in self.has_phases:
                 self._dof_indices_of_interface[ind].update(#
-                    {phase: self.interface[ind].dofs_on_interface(V[phase], marker, marker_value).copy()}
+                    {phase: self.interface[ind].dofs_on_interface(V[phase], marker, marker_value)}
                 )
-                print(f"subdomain{self.subdomain_index}: dofs on interface{ind}:", self._dof_indices_of_interface[ind][phase])
+                # print(f"subdomain{self.subdomain_index}: dofs on interface{ind}:", self._dof_indices_of_interface[ind][phase])
 
     def _calc_interface_has_common_dof_indices(self):
         """ populate dictionary self._interface_has_common_dof_indices
@@ -778,13 +792,16 @@ class DomainPatch(df.SubDomain):
             self._interface_has_common_dof_indices.update({interf_ind: dict()})
             if self.isRichards:
                 self._interface_has_common_dof_indices[interf_ind].update(
-                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w].copy()}
+                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w]}
                 )
             else:
                 self._interface_has_common_dof_indices[interf_ind].update(
-                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w].copy(),
-                    'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw].copy()}
+                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w],
+                    'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw]}
                 )
+                # for phase in self.has_phases:
+                #     print(f"self._interface_has_common_dof_indices[{interf_ind}][{phase}]", self._interface_has_common_dof_indices[interf_ind][phase])
+                #     print(f"python id of self._interface_has_common_dof_indices[{interf_ind}][{phase}]", id(self._interface_has_common_dof_indices[interf_ind][phase]))
 
 
     def write_solution_to_xdmf(self, file: tp.Type[LDDsimulation.SolutionFile], #
@@ -1045,7 +1062,6 @@ class Interface(BoundaryPart):
     """
     def __init__(self, #
                 global_index: int,#
-                # marker_value: int,#
                 adjacent_subdomains: np.array = None,#
                 # this is an array of bools
                 isRichards: np.array = None,#
@@ -1141,7 +1157,7 @@ class Interface(BoundaryPart):
         """
         vertex_indices = self._vertex_indices(interface_marker, #
                                               interface_marker_value,#
-                                              print_vertex_indices = True)
+                                              print_vertex_indices = False)
         self.pressure_values = dict()
         #self.pressure_values is a dictionay which contains for each index in
         # adjacent_subdomains a dictionary for wetting and nonwetting phase with