diff --git a/Usecases/README.md b/Usecases/README.md
index 0702c87f83d06cca44167001bd36683344146d0a..56ea5802b40ba67f298892baed9b5b4d425021b1 100644
--- a/Usecases/README.md
+++ b/Usecases/README.md
@@ -392,7 +392,13 @@ to have a central way to defining and reusing data functions.
 If you want to introduce other relative permeabilities, this is the place to do it.  
 - `Spc_on_subdomains`: The Spc get defined similarly to the relative permeablities by calling a a keyword specifying the type of function along with a dictionary of parameters. As in the case of the relative permeabilities
 the function `fts.generate_Spc_dicts(Spc_on_subdomains)` does the generation.
-Again, `../LDDsimulation/functions.py` is the place to start hacking to get other parametrisations going. 
+Again, `../LDDsimulation/functions.py` is the place to start hacking to get other parametrisations going.
+
+>>>
+**NOTE**
+Even if you are using TPR couplings, the coupling conditions need nonwetting information even on the Richards domain, see my thesis.
+Make sure to provide this information or the code will throw an error.
+>>>
 
 ### MANUFACTURED SOLUTION/SOURCE EXPRESSION
 ~~~python
@@ -403,11 +409,49 @@ x, y = sym.symbols('x[0], x[1]')  # needed by UFL
 t = sym.symbols('t', positive=True)
 
 p_e_sym = {
-    1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2},
-    2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2,
-        'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
-} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
+    1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y))},
+    2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
+        'nonwetting': (-2-t*(1.1+y + x**2))*y**2},
+}
+
+pc_e_sym = hlp.generate_exact_symbolic_pc(
+                isRichards=isRichards,
+                symbolic_pressure=p_e_sym
+            )
+
+symbols = {"x": x,
+           "y": y,
+           "t": t}
+# turn above symbolic code into exact solution for dolphin and
+# construct the rhs that matches the above exact solution.
+exact_solution_example = hlp.generate_exact_solution_expressions(
+                        symbols=symbols,
+                        isRichards=isRichards,
+                        symbolic_pressure=p_e_sym,
+                        symbolic_capillary_pressure=pc_e_sym,
+                        saturation_pressure_relationship=S_pc_sym,
+                        saturation_pressure_relationship_prime=S_pc_sym_prime,
+                        viscosity=viscosity,
+                        porosity=porosity,
+                        relative_permeability=relative_permeability,
+                        relative_permeability_prime=ka_prime,
+                        densities=densities,
+                        gravity_acceleration=gravity_acceleration,
+                        include_gravity=include_gravity,
+                        )
+source_expression = exact_solution_example['source']
+exact_solution = exact_solution_example['exact_solution']
+initial_condition = exact_solution_example['initial_condition']
 ~~~
+This block defines exact solution expressions and calculates source terms and
+initial conditions from it.
+If you dont want to have an exact solution, you can set
+`exact_solution = None`.
+Then you need to comment out the generation helper function specify
+`source_expression` and `initial_condition` directly similarly to `p_e_sym`.
+Have a look at the helper function `hlp.generate_exact_solution_expressions`
+specified in `../LDDsimulation/helpers.py`.
+
 
 ### BOUNDARY CONDITIONS
 ~~~python
@@ -424,3 +468,6 @@ dirichletBC = hlp.generate_exact_DirichletBC(
 # conditions. See the definiton of hlp.generate_exact_DirichletBC() to see
 # the structure.
 ~~~
+This block defines the boundary condition from the exact solution expression.
+Again, if no exact solution is assumed, `dirichletBC` needs to be specified
+manually, similar as explained above. 
diff --git a/Usecases/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/Usecases/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index 4f6d69518aa2b180264a8569d5c25b13a9085be9..a7943279ea450df9f1f9d353a0582a4bbe45df07 100755
--- a/Usecases/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/Usecases/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -231,10 +231,10 @@ x, y = sym.symbols('x[0], x[1]')  # needed by UFL
 t = sym.symbols('t', positive=True)
 
 p_e_sym = {
-    1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2},
-    2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2,
-        'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
-} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
+    1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y))},
+    2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
+        'nonwetting': (-2-t*(1.1+y + x**2))*y**2},
+}
 
 pc_e_sym = hlp.generate_exact_symbolic_pc(
                 isRichards=isRichards,