diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 1de1635cc8b7e867fcea926e4bd13bd5508e8a1c..8eed97a5091b7918aa76ff86d3bc80646b4c2626 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -132,7 +132,8 @@ class DomainPatch(df.SubDomain):
                  include_gravity: bool = False,#
                  gravity_acceleration: float = 9.81,
                  interpolation_degree: int = 10,
-                 tol = None,#
+                 tol: float = None,#
+                 degree: int = 1,
                  ):
         # because we declare our own __init__ method for the derived class BoundaryPart,
         # we overwrite the __init__-method of the parent class df.SubDomain. However,
@@ -163,6 +164,7 @@ class DomainPatch(df.SubDomain):
         # timestep size, tau in the paper
         self.timestep_size = timestep_size
         self.interpolation_degree = interpolation_degree
+        self.FEM_Lagrange_degree = degree
         ### Class variables set by methods
         # sometimes we need to loop over the phases and the length of the loop is
         # determined by the model we are assuming
@@ -674,16 +676,21 @@ class DomainPatch(df.SubDomain):
             pass
 
     #### PRIVATE METHODS
-    def _init_function_space(self) -> None:
+    def _init_function_space(self, degree: int = None) -> None:
         """ create function space for solution and trial functions
 
-        Note that P1 FEM is hard coded here, as it is assumed in other methods
-        aswell. This method sets the class variable
+        This method sets the class variable
             self.function_space["pressure"]
             self.trialfunction
             self.pressure
             self.testfunction
+
+        INPUT
+        degree      int     degree of the Lagrange ansatz space
         """
+        if degree is None:
+            degree = self.FEM_Lagrange_degree
+
         function_name = ["pressure", "gli", "flux"]
 
         self.function_space = dict()
@@ -692,16 +699,20 @@ class DomainPatch(df.SubDomain):
             if function == "pressure":
                 self.trialfunction = dict()
                 self.testfunction = dict()
-            # self.neighbouring_interface_pressure = dict()
+
+            # we want the same space for both primary variables (pw
+            # and pnw) to be able to calculate pnw-pw without issues.
+            pressure_space = df.FunctionSpace(self.mesh, 'P', degree)
             for phase in self.has_phases:
                 if function == "pressure":
-                    self.function_space[function].update({phase: df.FunctionSpace(self.mesh, 'P', 1)})
+                    self.function_space[function].update({phase: pressure_space})
                     self.trialfunction.update({phase: df.TrialFunction(self.function_space["pressure"][phase])})
                     self.testfunction.update({phase: df.TestFunction(self.function_space["pressure"][phase])})
                 elif function == "gli":
                     degree = self.function_space["pressure"][phase].ufl_element().degree()
                     self.function_space[function].update({phase: df.FunctionSpace(self.mesh, 'DG', degree)})
                 else:
+                    # here function = flux
                     degree = self.function_space["pressure"][phase].ufl_element().degree()
                     self.function_space[function].update({phase: df.VectorFunctionSpace(self.mesh, 'DG', degree)})
 
@@ -811,7 +822,7 @@ class DomainPatch(df.SubDomain):
         self.ds = df.Measure('ds', domain = self.mesh, subdomain_data = self.interface_marker)
 
 
-    def _calc_interface_dof_indices_and_coordinates(self, debug: bool = False):
+    def _calc_interface_dof_indices_and_coordinates(self, debug: bool = True):
         """ calculate dictionaries containing for each local facet index of
         interfaces a dictionary containing {interface_dof_index: dof_coordinates}
         """
@@ -847,11 +858,11 @@ class DomainPatch(df.SubDomain):
                         if debug:
                             print(f"\ndofs on interface for function space {function} for phase {phase}:")
                         self.interface_dof_indices_and_coordinates[function][interface_index].update(#
-                            {phase: interface.dofs_and_coordinates(function_space[phase], marker, marker_value, debug=False)}
+                            {phase: interface.dofs_and_coordinates(function_space[phase], marker, marker_value, debug=debug)}
                         )
 
 
-    def _calc_corresponding_dof_indices(self, debug=True):
+    def _calc_corresponding_dof_indices(self, debug=False):
         """ calculate dictionary which for each interface and each phase holds
         for each facet index, the dof indices of the pressures and flux components
         corresponding to a given dof index of the gli function.
@@ -1053,9 +1064,10 @@ class DomainPatch(df.SubDomain):
             if self.isRichards:
                 capillary_pressure.assign(-self.pressure["wetting"])
             else:
-                pc_temp = self.pressure["nonwetting"].vector()[:] - self.pressure["wetting"].vector()[:]
-                capillary_pressure.vector().set_local(pc_temp)
-                # capillary_pressure.assign(pc_temp)
+                # pc_temp = self.pressure["nonwetting"].vector()[:] - self.pressure["wetting"].vector()[:]
+                # capillary_pressure.vector().set_local(pc_temp)
+                pc_temp = self.pressure["nonwetting"] - self.pressure["wetting"]
+                capillary_pressure.assign(pc_temp)
             capillary_pressure.rename("pc_num", "pc_num")
             file.write(capillary_pressure, t)
         else: