diff --git a/domain_layered_soil.py b/domain_layered_soil.py
old mode 100644
new mode 100755
index 6ea0c8a788e26af8e0c79b2009cb09da03073c38..ba201a2ea45e359705b8ece149fb86b6fad931b8
--- a/domain_layered_soil.py
+++ b/domain_layered_soil.py
@@ -10,7 +10,7 @@ The resulting mesh is saved into files for later use.
 import dolfin as df
 import mshr
 
-set_log_level(1)
+df.set_log_level(1)
 
 domain = mshr.Rectangle(df.Point(0.0,0.0), df.Point(13.0,8.0))
 
@@ -28,13 +28,12 @@ domain.set_subdomain(1,sub_domain1)
 
 # subdomain 2
 sub_domain2_vertices = [df.Point(0.0, 5.0),
-                        df.Point(3.0, 5.0)
+                        df.Point(3.0, 5.0),
                         df.Point(6.5, 4.5),
                         df.Point(9.5, 5.0),
                         df.Point(11.5, 3.5),
                         df.Point(13.0, 3), # southern boundary, 23 interface
-                        df.Point(13.0, 6.5), # eastern boundary, outer boundary
-                        sub_domain1_vertices[4],
+                        sub_domain1_vertices[4], # eastern boundary, outer boundary
                         sub_domain1_vertices[3],
                         sub_domain1_vertices[2],
                         sub_domain1_vertices[1],
@@ -43,6 +42,39 @@ sub_domain2 = mshr.Polygon(sub_domain2_vertices)
 # set this as subdomain2
 domain.set_subdomain(2,sub_domain2)
 
+# subdomain 3
+sub_domain3_vertices = [df.Point(0.0, 2.0),
+                        df.Point(4.0, 2.0),
+                        df.Point(9.0, 2.5),
+                        df.Point(10.5, 2.0),
+                        df.Point(13.0, 1.5), # southern boundary, 34 interface
+                        sub_domain2_vertices[5], # eastern boundary, outer boundary
+                        sub_domain2_vertices[4],
+                        sub_domain2_vertices[3],
+                        sub_domain2_vertices[2],
+                        sub_domain2_vertices[1],
+                        sub_domain2_vertices[0] ] # northern boundary, 23 interface
+sub_domain3 = mshr.Polygon(sub_domain3_vertices)
+# set this as subdomain3
+domain.set_subdomain(3,sub_domain3)
+
+# subdomain 4
+sub_domain4_vertices = [df.Point(0.0, 0.0),
+                        df.Point(13.0, 0.0), # southern boundary, outer boundary
+                        sub_domain3_vertices[4],# eastern boundary, outer boundary
+                        sub_domain3_vertices[3],
+                        sub_domain3_vertices[2],
+                        sub_domain3_vertices[1],
+                        sub_domain3_vertices[0] ] # northern boundary, 34 interface
+sub_domain4 = mshr.Polygon(sub_domain4_vertices)
+# set this as subdomain4
+domain.set_subdomain(4,sub_domain4)
+
+#genererate mesh for the whole domains
+mesh = mshr.generate_mesh(domain,36)
+
+# Save mesh to file
+df.File('./domain_layered_soil.xml.gz') << mesh
 
 # Save sub domains to file
 #file = File("subdomains_layered_soil.xml")