diff --git a/LDDsimulation/domainSubstructuring.py b/LDDsimulation/domainSubstructuring.py
new file mode 100644
index 0000000000000000000000000000000000000000..09f5e72680fdc967be1c39af595949e310c935bc
--- /dev/null
+++ b/LDDsimulation/domainSubstructuring.py
@@ -0,0 +1,281 @@
+"""Class domainSubstructuring.
+
+Provides the base class and examples for defining substructurings for a
+the simulation domain
+"""
+import dolfin as df
+import helpers as hlp
+
+
+class domainSubstructuring(object):
+ """Base class for substructuring of a domain."""
+
+ def __init__(self):
+ """Set all variables."""
+ hlp.print_once("\n DD Substructuring:\n")
+ self.interface_def_points = None
+ self.adjacent_subdomains = None
+ self.subdomain_def_points = None
+ self.outer_boundary_def_points = None
+
+ def __interface_def_points(self):
+ """Set self._interface_def_points."""
+ raise(NotImplementedError())
+
+ def __adjacent_subdomains(self):
+ """Set self._adjacent_subdomains."""
+ raise(NotImplementedError())
+
+ def __subdomain_def_points(self):
+ """Set self._subdomain_def_points."""
+ raise(NotImplementedError())
+
+ def __outer_boundary_def_points(self):
+ """Set self._outer_boundary_def_points."""
+ raise(NotImplementedError())
+
+
+class layeredSoilInnerPatch(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")
+
+ 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."""
+ # 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.__interface12_vertices = [
+ df.Point(-1.0, 0.8),
+ df.Point(0.3, 0.8),
+ df.Point(0.5, 0.9),
+ df.Point(0.8, 0.7),
+ df.Point(1.0, 0.65)
+ ]
+
+ # interface23
+ self.__interface23_vertices = [
+ df.Point(-1.0, 0.0),
+ df.Point(-0.35, 0.0),
+ df.Point(0.0, 0.0)
+ ]
+
+ self.__interface24_vertices = [
+ self.__interface23_vertices[2],
+ df.Point(0.6, 0.0),
+ ]
+
+ self.__interface25_vertices = [
+ self.__interface24_vertices[1],
+ df.Point(1.0, 0.0)
+ ]
+
+ self.__interface32_vertices = [
+ self.__interface23_vertices[2],
+ self.__interface23_vertices[1],
+ self.__interface23_vertices[0]
+ ]
+
+ self.__interface36_vertices = [
+ df.Point(-1.0, -0.6),
+ df.Point(-0.6, -0.45)
+ ]
+
+ self.__interface46_vertices = [
+ self.__interface36_vertices[1],
+ df.Point(0.3, -0.25)
+ ]
+
+ self.__interface56_vertices = [
+ self.__interface46_vertices[1],
+ df.Point(0.65, -0.6),
+ df.Point(1.0, -0.7)
+ ]
+
+ self.__interface34_vertices = [
+ self.__interface36_vertices[1],
+ self.__interface23_vertices[2]
+ ]
+
+ # Interface 45 needs to be split, because of the shape. There can be
+ # triangles with two facets on the interface and this creates a rogue
+ # dof type error when integrating over that particular interface.
+ # Accordingly, the lambda_param dictionary has two entries for that
+ # interface.
+ self.__interface45_vertices_a = [
+ self.__interface56_vertices[0],
+ df.Point(0.7, -0.2),
+ ]
+
+ self.__interface45_vertices_b = [
+ df.Point(0.7, -0.2),
+ self.__interface25_vertices[0]
+ ]
+
+ # interface_vertices introduces a global numbering of interfaces.
+ self.interface_def_points = [
+ self.__interface12_vertices,
+ self.__interface23_vertices,
+ self.__interface24_vertices,
+ self.__interface25_vertices,
+ self.__interface34_vertices,
+ self.__interface36_vertices,
+ self.__interface45_vertices_a,
+ self.__interface45_vertices_b,
+ self.__interface46_vertices,
+ self.__interface56_vertices,
+ ]
+
+ def __adjacent_subdomains(self):
+ """Set self._adjacent_subdomains."""
+ self.adjacent_subdomains = [
+ [1, 2],
+ [2, 3],
+ [2, 4],
+ [2, 5],
+ [3, 4],
+ [3, 6],
+ [4, 5],
+ [4, 5],
+ [4, 6],
+ [5, 6]
+ ]
+
+ def __subdomain_def_points(self):
+ """Set self._subdomain_def_points."""
+ # subdomain1.
+ self.__subdomain1_vertices = [
+ self.__interface12_vertices[0],
+ self.__interface12_vertices[1],
+ self.__interface12_vertices[2],
+ self.__interface12_vertices[3],
+ self.__interface12_vertices[4],
+ self.__subdomain0_vertices[2],
+ self.__subdomain0_vertices[3]]
+
+ # subdomain1
+ self.__subdomain2_vertices = [
+ self.__interface23_vertices[0],
+ self.__interface23_vertices[1],
+ self.__interface23_vertices[2],
+ self.__interface24_vertices[1],
+ self.__interface25_vertices[1],
+ self.__subdomain1_vertices[4],
+ self.__subdomain1_vertices[3],
+ self.__subdomain1_vertices[2],
+ self.__subdomain1_vertices[1],
+ self.__subdomain1_vertices[0]]
+
+ self.__subdomain3_vertices = [
+ self.__interface36_vertices[0],
+ self.__interface36_vertices[1],
+ self.__interface34_vertices[1],
+ self.__interface32_vertices[1],
+ self.__interface32_vertices[2]
+ ]
+
+ # subdomain3
+ self.__subdomain4_vertices = [
+ self.__interface46_vertices[0],
+ self.__interface46_vertices[1],
+ self.__interface45_vertices_a[1],
+ self.__interface24_vertices[1],
+ self.__interface24_vertices[0],
+ self.__interface34_vertices[1]
+ ]
+
+ self.__subdomain5_vertices = [
+ self.__interface56_vertices[0],
+ self.__interface56_vertices[1],
+ self.__interface56_vertices[2],
+ self.__interface25_vertices[1],
+ self.__interface25_vertices[0],
+ self.__interface45_vertices_b[1],
+ self.__interface45_vertices_b[0]
+ ]
+
+ self.__subdomain6_vertices = [
+ self.__subdomain0_vertices[0],
+ self.__subdomain0_vertices[1],
+ self.__interface56_vertices[2],
+ self.__interface56_vertices[1],
+ self.__interface56_vertices[0],
+ self.__interface36_vertices[1],
+ self.__interface36_vertices[0]
+ ]
+
+ self.subdomain_def_points = [
+ self.__subdomain0_vertices,
+ self.__subdomain1_vertices,
+ self.__subdomain2_vertices,
+ self.__subdomain3_vertices,
+ self.__subdomain4_vertices,
+ self.__subdomain5_vertices,
+ self.__subdomain6_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.__subdomain1_vertices[4],
+ self.__subdomain1_vertices[5],
+ self.__subdomain1_vertices[6],
+ self.__subdomain1_vertices[0]]
+ }
+
+ self.__subdomain2_outer_boundary_verts = {
+ 0: [self.__subdomain2_vertices[9],
+ self.__subdomain2_vertices[0]],
+ 1: [self.__subdomain2_vertices[4],
+ self.__subdomain2_vertices[5]]
+ }
+
+ self.__subdomain3_outer_boundary_verts = {
+ 0: [self.__subdomain3_vertices[4],
+ self.__subdomain3_vertices[0]]
+ }
+
+ self.__subdomain4_outer_boundary_verts = None
+
+ self.__subdomain5_outer_boundary_verts = {
+ 0: [self.__subdomain5_vertices[2],
+ self.__subdomain5_vertices[3]]
+ }
+
+ self.__subdomain6_outer_boundary_verts = {
+ 0: [self.__subdomain6_vertices[6],
+ self.__subdomain6_vertices[0],
+ self.__subdomain6_vertices[1],
+ self.__subdomain6_vertices[2]]
+ }
+
+ # 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,
+ 3: self.__subdomain3_outer_boundary_verts,
+ 4: self.__subdomain4_outer_boundary_verts,
+ 5: self.__subdomain5_outer_boundary_verts,
+ 6: self.__subdomain6_outer_boundary_verts
+ }