diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 4e77d417531fc95e24a4116d364b633b45e33685..95138cf5700e63d4989f492106faac9458460a39 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -116,9 +116,9 @@ class LDDsimulation(object):
         # dictionary to hold the input Expressions for the external boundary. These are
         # turned into df.Expressions objects by the method update_DirichletBC_dictionary()
         # this is mostly in place when examples with exact solutions are calculated.
-        self.dirichletBC_Expressions = None
+        self.dirichletBC_expression_strings = None
         # dictionary to hold the df.Expressions for the external boundary created
-        # from the above self.dirichletBC_Expressions. These are
+        # from the above self.dirichletBC_expression_strings. These are
         # turned into df.dirichletBC objects by the method update_DirichletBC_dictionary()
         self.dirichletBC_dfExpression = dict()
         # dictionary to store the df.dirichletBC objects in. These are used by the
@@ -223,7 +223,7 @@ class LDDsimulation(object):
         self.timestep_size = timestep_size
         self.sources = sources
         self.initial_conditions = initial_conditions
-        self.dirichletBC_Expressions = dirichletBC_Expressions
+        self.dirichletBC_expression_strings = dirichletBC_Expressions
         self.exact_solution = exact_solution
         self.write2file = write2file
         self._parameters_set = True
@@ -946,13 +946,14 @@ class LDDsimulation(object):
             interpolation_degree = self.interpolation_degree
 
         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]
             for phase in subdomain.has_phases:
                 # note that the default interpolation degree is 2
-                pa0 = df.Expression(p0[phase], domain = mesh, degree = interpolation_degree, t=self.t)
+                pa0 = df.Expression(p0[phase], domain = subdomain.mesh, #
+                                               degree = interpolation_degree,#
+                                               t=self.starttime)
                 pa0_interpolated = df.interpolate(pa0, V[phase])
                 # alocate memory
                 subdomain.pressure.update(
@@ -984,7 +985,7 @@ class LDDsimulation(object):
             index subdomain_index.
 
         In case we don't have a case with an exact solution, we should have 0
-        expressions in self.dirichletBC_Expressions, so nothing should happen.
+        expressions in self.dirichletBC_expression_strings, so nothing should happen.
         """
         if interpolation_degree is None:
             interpolation_degree = self.interpolation_degree
@@ -994,7 +995,6 @@ class LDDsimulation(object):
         # we only need to set Dirichlet boundary conditions if we have an outer
         # boundary.
         if subdomain.outer_boundary is not None:
-            mesh = subdomain.mesh
             V = subdomain.function_space
             # loop over all outer boundary parts. This is a bit of a brainfuck:
             # subdomain.outer_boundary gives a dictionary of the outer boundaries of the subdomain.
@@ -1004,6 +1004,7 @@ class LDDsimulation(object):
             # the wetting and the nonwetting phase, so subdomain.outer_boundary[j]['wetting'] and
             # subdomain.outer_boundary[j]['nonwetting'] return
             # the actual expression needed for the dirichlet condition for both phases if present.
+            boundary_marker = subdomain.outer_boundary_marker
             for index, boundary in enumerate(subdomain.outer_boundary):
                 if np.abs(time-self.starttime) < self.tol:
                     # Here the dictionary has to be created first because t = 0.
@@ -1014,9 +1015,8 @@ class LDDsimulation(object):
                     # part of the subdomain.
                     self.outerBC[num].update({index: dict()})
                     self.dirichletBC_dfExpression[num].update({index: dict()})
-
-                pDC = self.dirichletBC_Expressions[num][index]
-                boundary_marker = subdomain.outer_boundary_marker
+                # this contains the Expression strings input by the user.
+                pDC = self.dirichletBC_expression_strings[num][index]
                 # as with the interfaces, since numbering of the outer boundary
                 # parts of each subdomain starts at 0 and subdomain.outer_boundary_marker
                 # gets set to 0 on all facets of the mesh, we need to up the value for
@@ -1026,11 +1026,13 @@ class LDDsimulation(object):
                     # note that the default interpolation degree is 2
                     if np.abs(time - self.starttime) < self.tol :
                         # time = 0 so the Dolfin Expression has to be created first.
-                        pDCa = df.Expression(pDC[phase], domain = mesh, degree = interpolation_degree, t=time)
+                        pDCa = df.Expression(pDC[phase], domain = subdomain.mesh, #
+                                                         degree = interpolation_degree, #
+                                                         t=self.starttime)
                         self.dirichletBC_dfExpression[num][index].update({phase: pDCa})
                     else:
-                        pDCa = self.dirichletBC_dfExpression[num][index][phase].t = time
-
+                        pDCa = self.dirichletBC_dfExpression[num][index][phase]
+                        pDCa.t = time
                     self.outerBC[num][index].update(
                         {phase: df.DirichletBC(V[phase], pDCa, boundary_marker, boundary_marker_value)}
                         )