From b95c0f321952bbc6248bdf7349daf4cd38547c1a Mon Sep 17 00:00:00 2001
From: David <forenkram@gmx.de>
Date: Mon, 29 Jun 2020 15:49:29 +0200
Subject: [PATCH] clean up TP-R 2 patch examples using functions module

---
 .../TP-R-2-patch-pure-dd-realistic.py         | 188 ++--------------
 .../TP-R-2-patch-pure-dd.py                   | 195 +++--------------
 ...alistic-gravity-but-same-intrinsic-perm.py | 203 ++---------------
 ...c-no-gravity-but-varying-intrinsic-perm.py | 200 ++---------------
 ...R-2-patch-realistic-same-intrinsic-perm.py | 204 ++----------------
 .../TP-R-2-patch-realistic.py                 | 190 ++--------------
 ...2-patch-mesh-study-all-params-one-but-g.py | 188 ++--------------
 .../TP-R-2-patch-mesh-study-all-params-one.py | 188 ++--------------
 .../mesh_studies/TP-R-2-patch-mesh-study.py   | 188 ++--------------
 ...vity-but-same-intrinsic-perm-mesh-study.py | 200 ++---------------
 ...y-but-varying-intrinsic-perm-mesh-study.py | 200 ++---------------
 ...R-2-patch-realistic-same-intrinsic-perm.py | 200 ++---------------
 12 files changed, 249 insertions(+), 2095 deletions(-)

diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
index d0d7874..7b96ff4 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
@@ -2,15 +2,13 @@
 
 This program sets up an LDD simulation
 """
-
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -187,176 +185,28 @@ intrinsic_permeability = {
         "nonwetting": 1},
 }
 
-## relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]["wetting"]*s**2
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]["nonwetting"]*(1-s)**2
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,#
-    'nonwetting': _rel_perm1nw
-}
-# ## relative permeabilty functions on subdomain 2
-# def rel_perm2w(s):
-#     # relative permeabilty wetting on subdomain2
-#     return s**3
-# def rel_perm2nw(s):
-#     # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
-#     return (1-s)**3
-#
-# _rel_perm2w = ft.partial(rel_perm2w)
-# _rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm1w,#
-    'nonwetting': _rel_perm1nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]["wetting"]*2*s
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]["nonwetting"]*2*(1-s)
-
-# # definition of the derivatives of the relative permeabilities
-# # relative permeabilty functions on subdomain 1
-# def rel_perm2w_prime(s):
-#     # relative permeabilty on subdomain1
-#     return 3*s**2
-#
-# def rel_perm2nw_prime(s):
-#     # relative permeabilty on subdomain1
-#     return -3*(1-s)**2
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-# _rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-# _rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=1),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=1),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=1),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 1}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 #############################################
 # Manufacture source expressions with sympy #
@@ -397,6 +247,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -453,6 +304,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
index 7e26fbf..2b9d2cd 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
@@ -2,15 +2,13 @@
 
 This program sets up an LDD simulation
 """
-
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -166,9 +164,9 @@ date = datetime.datetime.now()
 datestr = date.strftime("%Y-%m-%d")
 
 # Name of the usecase that will be printed during simulation.
-use_case = "TPR-2-patch-realistic"
+use_case = "TPR-2-patch-params-one-pure-dd"
 # The name of this very file. Needed for creating log output.
-thisfile = "TP-R-2-patch-realistic.py"
+thisfile = "TP-R-2-patch-pure-dd.py"
 
 # GENERAL SOLVER CONFIG  ######################################################
 # maximal iteration per timestep
@@ -315,176 +313,35 @@ lambda_param = {#
          'nonwetting': lambda_nw}
 }
 
-## relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return s**2
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return (1-s)**2
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,#
-    'nonwetting': _rel_perm1nw
-}
-# ## relative permeabilty functions on subdomain 2
-# def rel_perm2w(s):
-#     # relative permeabilty wetting on subdomain2
-#     return s**3
-# def rel_perm2nw(s):
-#     # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
-#     return (1-s)**3
-#
-# _rel_perm2w = ft.partial(rel_perm2w)
-# _rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm1w,#
-    'nonwetting': _rel_perm1nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
+intrinsic_permeability = {
+    1: {"wetting": 1,
+        "nonwetting": 1},
+    2: {"wetting": 1,
+        "nonwetting": 1},
 }
 
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return 2*s
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -2*(1-s)
-
-# # definition of the derivatives of the relative permeabilities
-# # relative permeabilty functions on subdomain 1
-# def rel_perm2w_prime(s):
-#     # relative permeabilty on subdomain1
-#     return 3*s**2
-#
-# def rel_perm2nw_prime(s):
-#     # relative permeabilty on subdomain1
-#     return -3*(1-s)**2
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-# _rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-# _rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=1),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=1),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=1),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 1}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 #############################################
 # Manufacture source expressions with sympy #
@@ -524,6 +381,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -580,6 +438,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py
index 1058674..fca9fe2 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py
@@ -3,19 +3,16 @@
 
 This program sets up an LDD simulation
 """
-
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
-
 # init sympy session
 sym.init_printing()
 
@@ -185,189 +182,28 @@ intrinsic_permeability = {
     2: 0.01,
 }
 
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,
-    'nonwetting': _rel_perm1nw
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,
-    'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -381,7 +217,6 @@ p_e_sym = {
         '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)
 
-
 pc_e_sym = hlp.generate_exact_symbolic_pc(
                 isRichards=isRichards,
                 symbolic_pressure=p_e_sym
@@ -401,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -457,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py
index f2f48ff..a2929e8 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
 """
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -183,189 +182,28 @@ intrinsic_permeability = {
     2: 0.0001,
 }
 
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,
-    'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,
-    'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
+# relative permeabilties
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -398,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -454,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py
index 92b123f..1b40691 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
 """
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -59,7 +58,7 @@ resolutions = {
 #  for each element t_0 in starttimes.
 starttimes = [0.0]
 timestep_size = 0.01
-number_of_timesteps = 800
+number_of_timesteps = 8
 
 # LDD scheme parameters  ######################################################
 Lw1 = 0.025
@@ -72,7 +71,7 @@ lambda_w = 4
 lambda_nw = 4
 
 include_gravity = False
-debugflag = False
+debugflag = True
 analyse_condition = False
 
 # I/O CONFIG  #################################################################
@@ -183,189 +182,28 @@ intrinsic_permeability = {
     2: 0.01,
 }
 
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,
-    'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,
-    'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
+# relative permeabilties
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -398,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -454,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py
index c95a5d1..6b16a1c 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py
@@ -3,15 +3,13 @@
 
 This program sets up an LDD simulation
 """
-
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -184,177 +182,28 @@ intrinsic_permeability = {
     2: 0.01,
 }
 
-## relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,#
-    'nonwetting': _rel_perm1nw
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,#
-    'nonwetting': _rel_perm2nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -368,7 +217,6 @@ p_e_sym = {
         '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)
 
-
 pc_e_sym = hlp.generate_exact_symbolic_pc(
                 isRichards=isRichards,
                 symbolic_pressure=p_e_sym
@@ -388,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -444,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py
index 113d0b8..d6bc8fe 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
 """
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -183,177 +182,28 @@ intrinsic_permeability = {
     2: 1.0,
 }
 
-## relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,#
-    'nonwetting': _rel_perm1nw
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,#
-    'nonwetting': _rel_perm2nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -386,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -442,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py
index cd3c859..c613e0b 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
 """
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -183,177 +182,28 @@ intrinsic_permeability = {
     2: 1.0,
 }
 
-## relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,#
-    'nonwetting': _rel_perm1nw
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,#
-    'nonwetting': _rel_perm2nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -386,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -442,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
index b69793d..ea48fc1 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
 """
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -183,177 +182,28 @@ intrinsic_permeability = {
     2: 1,
 }
 
-## relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,#
-    'nonwetting': _rel_perm1nw
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,#
-    'nonwetting': _rel_perm2nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -386,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -442,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py
index 2862c09..877dae6 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
 """
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -183,189 +182,28 @@ intrinsic_permeability = {
     2: 0.01,
 }
 
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,
-    'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,
-    'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -398,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -454,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py
index 0b305e1..249ab3f 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
 """
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -183,189 +182,28 @@ intrinsic_permeability = {
     2: 0.0001,
 }
 
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,
-    'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,
-    'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -398,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -454,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py
index d054ce4..5491080 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
 """
 import dolfin as df
 import sympy as sym
-import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
 import os
-import pandas as pd
 import multiprocessing as mp
 import domainSubstructuring as dss
 
@@ -183,189 +182,28 @@ intrinsic_permeability = {
     2: 0.01,
 }
 
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,
-    'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,
-    'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -3*intrinsic_permeability[2]*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-
-subdomain1_rel_perm_prime = {
-    'wetting': _rel_perm1w_prime,
-    'nonwetting': _rel_perm1nw_prime
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
 }
 
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
 
-subdomain2_rel_perm_prime = {
-    'wetting': _rel_perm2w_prime,
-    'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation1_sym(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_sym(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     #df.conditional(pc > 0,
-#     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
-#
-#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation1_sym_prime(pc, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-#     # inverse capillary pressure-saturation-relationship
-#     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
-# S_pc_sym = {
-#     1: ft.partial(saturation1_sym, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-#     1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-#     2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-#     1: ft.partial(saturation1, subdomain_index = 1),#,
-#     2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
-    # inverse capillary pressure-saturation-relationship
-    return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
-    1: ft.partial(saturation_sym, index=1),
-    2: ft.partial(saturation_sym, index=2),
-    # 3: ft.partial(saturation_sym, index=2),
-    # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=1),
-    2: ft.partial(saturation_sym_prime, index=2),
-    # 3: ft.partial(saturation_sym_prime, index=2),
-    # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
-    1: ft.partial(saturation, index=1),
-    2: ft.partial(saturation, index=2),
-    # 3: ft.partial(saturation, index=2),
-    # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+    1: {"testSpc": {"index": 1}},
+    2: {"testSpc": {"index": 2}},
 }
 
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
 
 ###############################################################################
 # Manufacture source expressions with sympy #
@@ -398,6 +236,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         saturation_pressure_relationship_prime=S_pc_sym_prime,
                         viscosity=viscosity,
                         porosity=porosity,
+                        intrinsic_permeability=intrinsic_permeability,
                         relative_permeability=relative_permeability,
                         relative_permeability_prime=ka_prime,
                         densities=densities,
@@ -454,6 +293,7 @@ if __name__ == '__main__':
         "L": L,
         "lambda_param": lambda_param,
         "relative_permeability": relative_permeability,
+        "intrinsic_permeability": intrinsic_permeability,
         "sat_pressure_relationship": sat_pressure_relationship,
         # "starttime": starttime,
         "number_of_timesteps": number_of_timesteps,
-- 
GitLab