define the whole domain and subdomains 1 and 2

domain_layered_soil.py

@@ -12,6 +12,38 @@ import mshr
 
 set_log_level(1)
 
+domain = mshr.Rectangle(df.Point(0.0,0.0), df.Point(13.0,8.0))
+
+# define vertices of the top most subdomain.
+sub_domain1_vertices = [df.Point(0.0, 7.0),
+                        df.Point(9.0, 7.0),
+                        df.Point(10.5, 7.5),
+                        df.Point(12.0, 7.0),
+                        df.Point(13.0, 6.5), # southern boundary, 12 interface
+                        df.Point(13.0, 8.0), # eastern boundary, outer boundary
+                        df.Point(0.0, 8.0) ] # northern boundary, outer on_boundary
+sub_domain1 = mshr.Polygon(sub_domain1_vertices)
+# set this as subdomain1
+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(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[3],
+                        sub_domain1_vertices[2],
+                        sub_domain1_vertices[1],
+                        sub_domain1_vertices[0] ] # northern boundary, 12 interface
+sub_domain2 = mshr.Polygon(sub_domain2_vertices)
+# set this as subdomain2
+domain.set_subdomain(2,sub_domain2)
+
+
 # Save sub domains to file
 #file = File("subdomains_layered_soil.xml")
 #file << sub_domains