diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index b630fa5f49c2f5a30248794dd5f4c29967d0611f..e75dccbe4f4d03f8799aa0eb411c3c74b417e095 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -5,6 +5,7 @@ LDD simulation class
 import os
 import dolfin as df
 import numpy as np
+import numpy.linalg as la
 import typing as tp
 import mshr
 import domainPatch as dp
@@ -293,7 +294,9 @@ class LDDsimulation(object):
         self._initialised = True
 
 
-    def run(self, analyse_solver_at_timesteps: tp.List[float] = None):
+    def run(self,
+            analyse_solver_at_timesteps: tp.List[float] = None,
+            analyse_condition: bool = False):
         """ run the simulation
 
         This method implements the time stepping, calls the acutal solvers and writes
@@ -335,9 +338,13 @@ class LDDsimulation(object):
             # check if the solver should beanalised or not.
             if analyse_solver_at_timesteps is not None and np.isin(self.timestep_num, analyse_solver_at_timesteps, assume_unique=True):
                 print("analyising timestep ...")
-                self.LDDsolver(debug = self.debug_solver, analyse_timestep = True)
+                self.LDDsolver(debug=self.debug_solver,
+                                analyse_timestep=True,
+                                analyse_condition=analyse_condition)
             else:
-                self.LDDsolver(debug = self.debug_solver, analyse_timestep = False)
+                self.LDDsolver(debug=self.debug_solver,
+                               analyse_timestep=False,
+                               analyse_condition=False)
 
             self.t += self.timestep_size
             self.timestep_num += 1
@@ -391,7 +398,8 @@ class LDDsimulation(object):
 
     def LDDsolver(self, time: float = None, #
                   debug: bool = False, #
-                  analyse_timestep: bool = False):
+                  analyse_timestep: bool = False,
+                  analyse_condition: bool = False):
         """ calulate the solution on all patches for the next timestep
 
         The heart of the simulation. The function LDDsolver implements the LDD
@@ -433,8 +441,10 @@ class LDDsimulation(object):
                 # set subdomain iteration to current iteration
                 subdomain.iteration_number = iteration
                 # solve the problem on subdomain
-                self.Lsolver_step(subdomain_index = sd_index,#
-                                  debug = debug,
+                condition_nums = self.Lsolver_step(
+                                  subdomain_index=sd_index,#
+                                  debug=debug,
+                                  analyse_condition=analyse_condition
                                   )
                 # bypass subsequent error calculation for the first iteration, becauses
                 # it is zero anyways.
@@ -464,8 +474,20 @@ class LDDsimulation(object):
                     self.write_subsequent_errors_to_csv(
                         filename = subsequent_error_filename, #
                         subdomain_index = sd_index,
-                        errors = subsequent_iter_error,
-                        time = time)
+                        errors = subsequent_iter_error
+                        )
+                    if analyse_condition:
+                        # save the calculated condition numbers of the assembled
+                        # matrices a separate file for monitoring
+                        condition_number_filename = self.output_dir\
+                            +self.output_filename_parameter_part[sd_index]\
+                            +"condition_numbers" +"_at_timestep"+\
+                            "{number}".format(number=self.timestep_num) +".csv"
+                        self.write_condition_numbers_to_csv(
+                            filename = condition_number_filename, #
+                            subdomain_index = sd_index,
+                            condition_number = condition_nums
+                            )
 
                 # prepare next iteration
                 # write the newly calculated solution to the inteface dictionaries
@@ -496,7 +518,7 @@ class LDDsimulation(object):
 
             total_subsequent_error = np.amax(max_subsequent_error)
             if debug:
-                print(f"time = {time}: max subsequent error of both subdomain in iteration {iteration} = ", total_subsequent_error)
+                print(f"time = {time} (at timestep {self.timestep_num}):: max subsequent error of both subdomain in iteration {iteration} = ", total_subsequent_error)
             if total_subsequent_error < error_criterion:
                 # if debug:
                 print(f"all phases on all subdomains reached a subsequent error = {total_subsequent_error} < tol = {error_criterion} after {iteration} iterations.")
@@ -602,7 +624,8 @@ class LDDsimulation(object):
 
 
     def Lsolver_step(self, subdomain_index: int,#
-                     debug: bool = False) -> None:
+                     analyse_condition: bool = False,
+                     debug: bool = False) -> tp.Dict[str, float]:
         """ L-scheme solver iteration step for an object of class subdomain
 
         Lsolver_step implements L-scheme solver iteration step for an subdomain
@@ -619,7 +642,9 @@ class LDDsimulation(object):
         """
         subdomain = self.subdomain[subdomain_index]
         iteration = subdomain.iteration_number
-
+        condition_number = None
+        if analyse_condition:
+            condition_number = dict()
         for phase in subdomain.has_phases:
             # extract L-scheme form and rhs (without gli term) from subdomain.
             governing_problem = subdomain.governing_problem(phase = phase)
@@ -629,6 +654,11 @@ class LDDsimulation(object):
             form_assembled = df.assemble(form)
             rhs_assembled = df.assemble(rhs)
 
+            if analyse_condition:
+                condition_number.update(
+                    {phase: la.cond(form_assembled.array())}
+                )
+
             # apply outer Dirichlet boundary conditions if present.
             if subdomain.outer_boundary is not None:
                 # apply boundary dirichlet conditions for all outer boundaries.
@@ -650,6 +680,8 @@ class LDDsimulation(object):
             if debug and subdomain.mesh.num_cells() < 36:
                 print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
 
+        return condition_number
+
 
     def _init_solution_files(self):
         """ set up solution files for saving the solution of the LDD scheme for
@@ -714,12 +746,8 @@ class LDDsimulation(object):
                                 filename: str, #
                                 subdomain_index: int,#
                                 errors: tp.Dict[str, float],#
-                                time: float = None) -> None:
+                                ) -> None:
         """ write subsequent errors to csv file for plotting with pgfplots"""
-        # in case no time is given, use the time of the class.
-        subdomain = self.subdomain[subdomain_index]
-        if time is None:
-            time = self.t
 
         subdomain = self.subdomain[subdomain_index]
         iteration = subdomain.iteration_number
@@ -744,6 +772,37 @@ class LDDsimulation(object):
         else:
             self.errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
 
+
+    def write_condition_numbers_to_csv(self,#
+                                filename: str, #
+                                subdomain_index: int,#
+                                condition_number: tp.Dict[str, float],#
+                                ) -> None:
+        """ write subsequent errors to csv file for plotting with pgfplots"""
+        subdomain = self.subdomain[subdomain_index]
+        iteration = subdomain.iteration_number
+
+        if subdomain.isRichards:
+            condition = pd.DataFrame({ "iteration": iteration,
+                #"th. initial mass": df.assemble(rho*df.dx), \
+                #"energy": self.total_energy(), \
+                "wetting": condition_number["wetting"]}, index=[iteration])
+
+        else:
+            condition = pd.DataFrame({"iteration":iteration, #\
+                #"th. initial mass": df.assemble(rho*df.dx), \
+                #"energy": self.total_energy(), \
+                "wetting": condition_number["wetting"],
+                "nonwetting": condition_number["nonwetting"]}, index=[iteration])
+        # if self._rank == 0:
+        # check if file exists
+        if os.path.isfile(filename) == True:
+            with open(filename, 'a') as f:
+                condition.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
+        else:
+            condition.to_csv(filename, sep='\t', encoding='utf-8', index=False)
+
+
     def write_errornorms_to_csv(self,#
                                 filename: str, #
                                 subdomain_index: int,#
diff --git a/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index ef696e314ebb02f3793138c7ea766558dc4f01c0..99241158d9a982114ab23a6bb77dcbf4881628c3 100755
--- a/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -14,7 +14,7 @@ import helpers as hlp
 sym.init_printing()
 
 solver_tol = 5e-7
-############ GRID #######################ü
+######################## GRID #######################
 mesh_resolution = 30
 timestep_size = 0.0001
 number_of_timesteps = 50
diff --git a/TP-one-patch/TP-one-patch.py b/TP-one-patch/TP-one-patch.py
index c66fc33dca088bc975336061d86ee3a0a0452ae8..4df57fc61610efbc48c960cdd4023132b774a7ad 100755
--- a/TP-one-patch/TP-one-patch.py
+++ b/TP-one-patch/TP-one-patch.py
@@ -14,25 +14,28 @@ import helpers as hlp
 sym.init_printing()
 
 
-solver_tol = 5E-6
+solver_tol = 1E-7
 
 ############ GRID #######################ü
-mesh_resolution = 6
-timestep_size = 0.001
+mesh_resolution = 20
+timestep_size = 0.0001
 number_of_timesteps = 20
 # decide how many timesteps you want analysed. Analysed means, that we write out
 # subsequent errors of the L-iteration within the timestep.
 number_of_timesteps_to_analyse = 10
 starttime = 0
 
-Lw = 1/timestep_size
+Lw = 1  #/timestep_size
 Lnw=Lw
 
-l_param_w = 80
-l_param_nw = 80
+l_param_w = 40
+l_param_nw = 40
 
 include_gravity = False
-debugflag = False
+debugflag = True
+analyse_condition = False
+
+output_string = "./output/nondirichlet_number_of_timesteps{}_".format(number_of_timesteps)
 
 ##### Domain and Interface ####
 # global simulation domain domain
@@ -283,16 +286,16 @@ cutoff = gaussian/(gaussian + zero_on_shrinking)
 #     is_inside = ((1-sym.Abs(x) > epsilon) & (1-sym.Abs(y) > epsilon))
 #     sym.Piecewise((0, is_inside))
 
-# p_e_sym = {
-#     0: {'wetting': (-3 - (1+t*t)*(1 + x*x + y*y))*cutoff,
-#         'nonwetting': (-1 -t*(1+y + x**2)**2)*cutoff},
-# }
-
 p_e_sym = {
-    0: {'wetting': -(sym.cos(2*t-x - 2*y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2,
-        'nonwetting': -6*(sym.cos(t-x -y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2},
+    0: {'wetting': (-3 - (1+t*t)*(1 + x*x + y*y)),  #*cutoff,
+        'nonwetting': (-1 -t*(1+y + x**2)**2)},  #*cutoff},
 }
 
+# p_e_sym = {
+#     0: {'wetting': -(sym.cos(2*t-x - 2*y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2,
+#         'nonwetting': -6*(sym.cos(t-x -y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2},
+# }
+
 
 print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n")
 # # pw_sym_x*pw_sym_y
@@ -310,10 +313,10 @@ print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n")
 pc_e_sym = dict()
 for subdomain, isR in isRichards.items():
     if isR:
-        pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
+        pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
     else:
-        pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
-                                        - p_e_sym[subdomain]['wetting'].copy()})
+        pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting']
+                                        - p_e_sym[subdomain]['wetting']})
 
 symbols = {"x": x,
            "y": y,
@@ -379,8 +382,6 @@ write_to_file = {
     'L_iterations': True
 }
 
-output_string = "./output/testing_number_of_timesteps{}_".format(number_of_timesteps)
-
 # initialise LDD simulation class
 simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=debugflag)
 simulation.set_parameters(output_dir = output_string,#
@@ -411,4 +412,4 @@ simulation.set_parameters(output_dir = output_string,#
 
 simulation.initialise()
 # simulation.write_exact_solution_to_xdmf()
-simulation.run()
+simulation.run(analyse_condition=analyse_condition)