diff --git a/LDDsimulation.py b/LDDsimulation.py
index b86e73a3c9bc3270654bdc10dbebcfba493d4572..b2e2d112a8854a0a62a00786d7f9963bdc92d471 100644
--- a/LDDsimulation.py
+++ b/LDDsimulation.py
@@ -328,6 +328,27 @@ class LDDsimulation(object):
                 subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting']),#
                                                     'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
 
+    def _init_DirichletBC(self, interpolation_degree: int = 2):
+        """ set initial values
+        """
+        for num, subdomain in self.subdomain.items():
+            mesh = subdomain.mesh
+            V = subdomain.function_space
+            # p0 contains both pw_0 and pnw_0 for subdomain[num]
+            p0 = self.initial_conditions[num]
+            if subdomain.isRichards:
+                # note that the default interpolation degree is 2
+                pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
+                subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting'])}
+                subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting'])}
+            else:
+                pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
+                pnw0 = df.Expression(p0['nonwetting'], domain = mesh, degree = interpolation_degree)
+                subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting']),#
+                                                'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
+                subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting']),#
+                                                    'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
+
     def _eval_sources(self, interpolation_degree: int = 2, time: float = 0):
         """ evaluate time dependent source terms or initialise them if time == 0
         """
diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
index bb90d7d423475810870fa10ec2a946c48e6e32f6..bf9c3f79fb994379d756e1fc2c8e53b758031e0d 100755
--- a/RR-2-patch-test.py
+++ b/RR-2-patch-test.py
@@ -125,6 +125,13 @@ initial_condition = {
     2: {'wetting': '-x[1]*x[1]'}
 }
 
+exact_solution = {
+    1: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'},#
+    2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}
+}
+
+dirichletBC = exact_solution
+
 # def saturation(pressure, subdomain_index):
 #     # inverse capillary pressure-saturation-relationship
 #     return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
diff --git a/domainPatch.py b/domainPatch.py
index ad563c4025fda86a6dc232614b3fc562f519cd62..5e298c6a543aab94d31f07078b7f6ef4301a6bad 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -139,7 +139,7 @@ class DomainPatch(df.SubDomain):
         ### Class variables set by methods
         # 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()
+        self._dof_indices_of_interface = dict()
         # measures are set by self._init_measures()
         self.iteration_number = 0
         self.dx = None
@@ -165,7 +165,7 @@ class DomainPatch(df.SubDomain):
     def write_pressure_to_interfaces(self):
         """ save the interface values of self.pressure to all neighbouring interfaces"""
         if self.isRichards:
-            for ind in has_interface:
+            for ind in self.has_interface:
                 interface[ind].write_dofs(from_function = self.pressure['wetting'], #
                                         interface_dofs = self._dof_indices_of_interface[ind]['wetting'],#
                                         dof_to_vert_map = self.dof2vertex['wetting'],#
@@ -173,7 +173,7 @@ class DomainPatch(df.SubDomain):
                                         phase = 'wetting',#
                                         subdomain_ind = self.subdomain_index)
         else:
-            for ind in has_interface:
+            for ind in self.has_interface:
                 for phase in ['wetting', 'nonwetting']:
                     interface[ind].write_dofs(from_function = self.pressure[phase], #
                                             interface_dofs = self._dof_indices_of_interface[ind][phase],#
@@ -311,7 +311,7 @@ class DomainPatch(df.SubDomain):
         """ calculate dof indices of each interface """
         V = self.function_space
         marker = self.boundary_marker
-        for ind in has_interface:
+        for ind in self.has_interface:
             self._dof_indices_of_interface.update(#
                 {ind: dict()}
             )