diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 16e7d366ee803bd55f7c6fa4c47a752d69024770..a0f597649f67254520f466268ea07ec4225fd5f7 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -90,8 +90,27 @@ class LDDsimulation(object):
         self._parameters_set = False
         # dictionary of objects of class DomainPatch initialised by self._init_subdomains()
         self.subdomain = dict()
-        # dictionary to store the dirichletBC in.
+        # dictionary to hold the input Expressions for the external boundary. These are
+        # turned into df.Expressions objects by the method update_DirichletBC_dictionaries()
+        # this is mostly in place when examples with exact solutions are calculated.
+        self.dirichletBC_Expressions = None
+        # dictionary to hold the df.Expressions for the external boundary created
+        # from the above self.dirichletBC_Expressions. These are
+        # turned into df.dirichletBC objects by the method update_DirichletBC_dictionaries()
+        self.dirichletBC_dfExpression = dict()
+        # dictionary to store the df.dirichletBC objects in. These are used by the
+        # solver to set the boundary conditions
         self.outerBC = dict()
+        # Flag to toggle between a case with exact solution and a purely numerical
+        # simulation. This changes the behaviour of self.update_DirichletBC_dictionaries()
+        self.exact_solution_case = False
+        # dictionary holding the source expressions for each subdomain. These are
+        # saved to the respective subdomains by the self._init_subdomains() method.
+        self.sources = None#
+        # dictionary holding expressions for the initial conditions for each subdomain.
+        # These are interpolated and stored to each respective subdomain by
+        # self._init_initial_values()
+        self.initial_conditions = None
 
     def set_parameters(self,#
                        output_dir: str,#
@@ -111,7 +130,8 @@ class LDDsimulation(object):
                        timestep_size: tp.Dict[int, float],#
                        sources: tp.Dict[int, tp.Dict[str, str]],#
                        initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
-                       outer_dirichletBC: tp.Dict[int, tp.Dict[str, str]],#
+                       dirichletBC_Expressions: tp.Dict[int, tp.Dict[str, str]],#
+                       exact_solution_case: bool = False#
                        )-> None:
 
         """ set parameters of an instance of class LDDsimulation"""
@@ -131,7 +151,8 @@ class LDDsimulation(object):
         self.timestep_size = timestep_size
         self.sources = sources
         self.initial_conditions = initial_conditions
-        self.outer_dirichletBC = outer_dirichletBC
+        self.dirichletBC_Expressions = dirichletBC_Expressions
+        self.exact_solution_case = exact_solution_case
         self._parameters_set = True
 
     def initialise(self) -> None:
@@ -153,24 +174,27 @@ class LDDsimulation(object):
         iteration, that is solves the problem given the initial data at timestep
         time.
         """
-        iteration = 0
         # dictionary holding bool variables for each subdomain indicating wether
         # minimal error has been achieved, i.e. the calculation can be considered
         # done or not.
         subdomain_calc_isdone = dict()
         # before the iteration gets started all iteration numbers must be nulled
-        # and
+        # and the subdomain_calc_isdone dictionary must be set to False
         for ind, subdomain in self.subdomain.items():
             subdomain_calc_isdone.update({ind, False})
+            # reset all interface[has_interface].current_iteration[ind] = 0, for
+            # index has_interface in subdomain.has_interface
+            subdomain.null_all_interface_iteration_numbers()
             # update the time of the source terms before starting the iteration.
             # The sources and sinks are the same throughout the iteration.
             self._eval_sources(interpolation_degree = 2,
                                time = time,
                                subdomain_index = ind)
+            self.update_DirichletBC_dictionaries()
             # the following only needs to be done when time is not equal 0 since
             # our initialisation functions already took care of setting what follows
-            # for the initial iteration step.
-            if time > 0:
+            # for the initial iteration step at time = 0.
+            if time > self.tol:
                 # this is the beginning of the solver, so iteration must be set
                 # to 0.
                 subdomain.iteration_number = 0
@@ -182,8 +206,9 @@ class LDDsimulation(object):
 
         ### actual iteration starts here
         # gobal stopping criterion for the iteration.
+        iteration = 0
         all_subdomains_are_done = False
-        while iteration < self._max_iter_num and not all_subdomains_are_done:
+        while iteration <= self._max_iter_num and not all_subdomains_are_done:
             # we need to loop over the subdomains and solve an L-scheme type step
             # on each subdomain. here it should be possible to parallelise in
             # the future.
@@ -191,9 +216,10 @@ class LDDsimulation(object):
                 # check if the calculation on subdomain has already be marked as
                 # finished.
                 if not subdomain_calc_isdone[sd_index]:
+                    # set subdomain iteration to current iteration
+                    subdomain.iteration_number = iteration
                     # solve the problem on subdomain
                     self.Lsolver_step(subdomain_index = sd_index,#
-                                      iteration = iteration,#
                                       debug = True
                     )
                     subsequent_iterations_err = self.calc_iteration_error(
@@ -229,8 +255,7 @@ class LDDsimulation(object):
             iteration += 1
             # end iteration while loop.
 
-    def Lsolver_step(subdomain_index: int],#
-                     iteration: int,#
+    def Lsolver_step(self, subdomain_index: int],#
                      debug: bool = False) -> None:
         """ L-scheme solver iteration step for an object of class subdomain
 
@@ -246,6 +271,26 @@ class LDDsimulation(object):
         debug:              bool                        flag to toggle weather or not
                                                         to write out each iteration.
         """
+        subdomain = self.subdomain[subdomain_index]
+        iteration = subdomain.iteration_number
+        if subdomain.isRichards:
+            # extract L-scheme form and rhs (without gli term) from subdomain.
+            governing_problem = subdomain.governing_problem(phase = 'wetting')
+            form = governing_problem[form]
+            rhs_without_gli = governing_problem[rhs_without_gli]
+            # assemble the form and rhs
+            form_assembled = df.assemble(form)
+            rhs_without_gli_assembled = df.assemble(rhs_without_gli)
+            # calculate the gli term
+            gli_term_assembled = subdomain.calc_gli_term(iteration = iteration)
+            # subdomain.calc_gli_term() asslembles gli but on the left hand side
+            # so gli_term_assembled needs to actually be subtracted from the rhs.
+            rhs_assembled = rhs_without_gli_assembled - gli_term_assembled
+            # apply outer Dirichlet boundary conditions if present.
+
+
+        else:
+            print("implement two phase assembly")
 
 
 
@@ -446,32 +491,55 @@ class LDDsimulation(object):
             # populate interface dictionaries with pressure values.
             subdomain.write_pressure_to_interfaces(initial_values = True)
 
-    def set_DirichletBC(self, interpolation_degree: int = 2, #
+    def update_DirichletBC_dictionary(self,
+                        interpolation_degree: int = 2, #
                         time: float = 0,#
                         # exact_solution: bool = False,#
-                        first_init: bool = False):
-        """ set dirichlet boundary values at time time.
-        """
-        for num, subdomain in self.subdomain.items():
-            mesh = subdomain.mesh
-            V = subdomain.function_space
-            # if first_init==True we have not initialised the dictionary
-            if first_init:
-                self.outerBC.update({num: dict()})
+                        subdomain_index: int):
+        """ update time of the dirichlet boundary object for subdomain with
+            index subdomain_index.
 
-            pDC = self.outer_dirichletBC[num]
-            boundary_marker = subdomain.outer_boundary_marker
-            if subdomain.isRichards:
-                # note that the default interpolation degree is 2
+        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.
+        """
+        num = subdomain_index
+        subdomain = self.subdomain[num]
+        mesh = subdomain.mesh
+        V = subdomain.function_space
+        # Here the dictionary has to be created first because t = 0.
+        if np.abs(time) < self.tol :
+            self.outerBC.update({num: dict()})
+            self.dirichletBC_dfExpression.update({num: dict()})
+
+        pDC = self.dirichletBC_Expressions[num]
+        boundary_marker = subdomain.outer_boundary_marker
+        if subdomain.isRichards:
+            # note that the default interpolation degree is 2
+            if np.abs(time) < self.tol :
+                # time = 0 so the Dolfin Expression has to be created first.
                 pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
-                self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)})
+                self.dirichletBC_dfExpression[num].update({'wetting': pDCw})
             else:
+                pDCw = self.dirichletBC_dfExpression[num]['wetting'].t = time
+
+            self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)})
+        else:
+            if np.abs(time) < self.tol :
+                # time = 0 so the Dolfin Expression has to be created first.
                 pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
                 pDCnw = df.Expression(pDC['nonwetting'], domain = mesh, degree = interpolation_degree, t=time)
-                self.outerBC[num].update(#
-                    {'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)},#
-                    {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 1)},#
-                )
+                self.dirichletBC_dfExpression[num].update(
+                    {'wetting': pDCw},#
+                    {'nonwetting': pDCnw},#
+                    )
+            else:
+                pDCw = self.dirichletBC_dfExpression[num]['wetting'].t = time
+                pDCnw = self.dirichletBC_dfExpression[num]['nonwetting'].t = time
+
+            self.outerBC[num].update(#
+                {'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)},#
+                {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 1)},#
+            )
 
     def _eval_sources(self, interpolation_degree: int = 2, #
                             time: float = 0,