diff --git a/LDDsimulation/domainSubstructuring.py b/LDDsimulation/domainSubstructuring.py
index 09f5e72680fdc967be1c39af595949e310c935bc..4048ef4957c74f0925c383a70008d4150fd5cbf7 100644
--- a/LDDsimulation/domainSubstructuring.py
+++ b/LDDsimulation/domainSubstructuring.py
@@ -34,6 +34,99 @@ class domainSubstructuring(object):
         """Set self._outer_boundary_def_points."""
         raise(NotImplementedError())
 
+class twoSoilLayers(domainSubstructuring):
+    """layered soil substructuring with inner patch."""
+
+    def __init__(self):
+        """Layered soil case with inner patch."""
+        super().__init__()
+        hlp.print_once("\n Layered Soil with inner Patch:\n")
+        # global domain
+        self.__subdomain0_vertices = [
+            df.Point(-1.0, -1.0),
+            df.Point(1.0, -1.0),
+            df.Point(1.0, 1.0),
+            df.Point(-1.0, 1.0)
+            ]
+
+        self.__interface_def_points()
+        self.__adjacent_subdomains()
+        self.__subdomain_def_points()
+        self.__outer_boundary_def_points()
+
+    def __interface_def_points(self):
+        """Set self._interface_def_points."""
+
+        self.__interface12_vertices = [
+            df.Point(-1.0, 0.0),
+            df.Point(1.0, 0.0)
+            ]
+
+        # interface_vertices introduces a global numbering of interfaces.
+        self.interface_def_points = [
+            self.__interface12_vertices,
+            ]
+
+    def __adjacent_subdomains(self):
+        """Set self._adjacent_subdomains."""
+        self.adjacent_subdomains = [
+            [1, 2],
+            ]
+
+    def __subdomain_def_points(self):
+        """Set self._subdomain_def_points."""
+        # subdomain1.
+        self.__subdomain1_vertices = [
+            self.__interface12_vertices[0],
+            self.__interface12_vertices[1],
+            self.__sub_domain0_vertices[2],
+            self.__sub_domain0_vertices[3]
+            ]
+
+        # subdomain2
+        self.__subdomain2_vertices = [
+            self.__sub_domain0_vertices[0],
+            self.__sub_domain0_vertices[1],
+            self.__interface12_vertices[1],
+            self.__interface12_vertices[0]]
+
+        self.subdomain_def_points = [
+            self.__subdomain0_vertices,
+            self.__subdomain1_vertices,
+            self.__subdomain2_vertices,
+            ]
+
+    def __outer_boundary_def_points(self):
+        """Set self._outer_boundary_def_points."""
+        # vertex coordinates of the outer boundaries. If it can not be
+        # specified as a polygon, use an entry per boundary polygon.
+        # This information is used for defining the Dirichlet boundary
+        # conditions. If a domain is completely internal, the
+        # dictionary entry should be 0: None
+        self.__subdomain1_outer_boundary_verts = {
+            0: [self.__interface12_vertices[1],
+                self.__sub_domain0_vertices[2],
+                self.__sub_domain0_vertices[3],
+                self.__interface12_vertices[0]]
+        }
+
+        self.__subdomain2_outer_boundary_verts = {
+            0: [self.__interface12_vertices[0],
+                self.__sub_domain0_vertices[0],
+                self.__sub_domain0_vertices[1],
+                self.__interface12_vertices[1]]
+        }
+
+        # if a subdomain has no outer boundary write None instead, i.e.
+        # i: None
+        # if i is the index of the inner subdomain.
+        self.outer_boundary_def_points = {
+            # subdomain number
+            1: self.__subdomain1_outer_boundary_verts,
+            2: self.__subdomain2_outer_boundary_verts,
+        }
+
+
 
 class layeredSoilInnerPatch(domainSubstructuring):
     """layered soil substructuring with inner patch."""
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 4dd2660ab918798f2e03a3231751aa4f1d2bada3..006b47545bb2dda96e7e90fce516348788edc310 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
@@ -42,12 +42,12 @@ max_iter_num = 5
 FEM_Lagrange_degree = 1
 
 # GRID AND MESH STUDY SPECIFICATIONS  #########################################
-mesh_study = True
+mesh_study = False
 resolutions = {
                 # 1: 1e-5,
                 # 2: 1e-5,
-                4: 1e-5,
-                8: 1e-5,
+                # 4: 1e-5,
+                # 8: 1e-5,
                 16: 5e-6,
                 # 32: 5e-6,
                 # 64: 2e-6,
@@ -135,67 +135,11 @@ else:
 
 
 # DOMAIN AND INTERFACE  #######################################################
-# global simulation domain domain
-sub_domain0_vertices = [df.Point(-1.0, -1.0),
-                        df.Point(1.0, -1.0),
-                        df.Point(1.0, 1.0),
-                        df.Point(-1.0, 1.0)]
-# interface between subdomain1 and subdomain2
-interface12_vertices = [df.Point(-1.0, 0.0),
-                        df.Point(1.0, 0.0) ]
-# subdomain1.
-sub_domain1_vertices = [interface12_vertices[0],
-                        interface12_vertices[1],
-                        sub_domain0_vertices[2],
-                        sub_domain0_vertices[3]]
-
-# vertex coordinates of the outer boundaries. If it can not be specified as a
-# polygon, use an entry per boundary polygon. This information is used for defining
-# the Dirichlet boundary conditions. If a domain is completely internal, the
-# dictionary entry should be 0: None
-subdomain1_outer_boundary_verts = {
-    0: [interface12_vertices[1], #
-        sub_domain0_vertices[2],
-        sub_domain0_vertices[3], #
-        interface12_vertices[0]]
-}
-# subdomain2
-sub_domain2_vertices = [sub_domain0_vertices[0],
-                        sub_domain0_vertices[1],
-                        interface12_vertices[1],
-                        interface12_vertices[0] ]
-
-subdomain2_outer_boundary_verts = {
-    0: [interface12_vertices[0], #
-        sub_domain0_vertices[0],
-        sub_domain0_vertices[1],
-        interface12_vertices[1]]
-}
-
-# list of subdomains given by the boundary polygon vertices.
-# Subdomains are given as a list of dolfin points forming
-# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
-# to create the subdomain. subdomain_def_points[0] contains the
-# vertices of the global simulation domain and subdomain_def_points[i] contains the
-# vertices of the subdomain i.
-subdomain_def_points = [sub_domain0_vertices,#
-                      sub_domain1_vertices,#
-                      sub_domain2_vertices]
-# in the below list, index 0 corresponds to the 12 interface which has index 1
-interface_def_points = [interface12_vertices]
-
-# if a subdomain has no outer boundary write None instead, i.e.
-# i: None
-# if i is the index of the inner subdomain.
-outer_boundary_def_points = {
-    # subdomain number
-    1 : subdomain1_outer_boundary_verts,
-    2 : subdomain2_outer_boundary_verts
-}
-
-# adjacent_subdomains[i] contains the indices of the subdomains sharing the
-# interface i (i.e. given by interface_def_points[i]).
-adjacent_subdomains = [[1,2]]
+substructuring = dss.twoSoilLayers()
+interface_def_points = substructuring.interface_def_points
+adjacent_subdomains = substructuring.adjacent_subdomains
+subdomain_def_points = substructuring.subdomain_def_points
+outer_boundary_def_points = substructuring.outer_boundary_def_points
 
 # MODEL CONFIGURATION #########################################################
 isRichards = {