Changeset 2036 for inundation/analytical solutions/Non_symmetrical_dam_break.py

Timestamp:

Nov 17, 2005, 4:16:50 PM (19 years ago)

Author:

steve

Message:

Adding ways to create cut out regions

File:

• inundation/analytical solutions/Non_symmetrical_dam_break.py (modified) (8 diffs)

inundation/analytical solutions/Non_symmetrical_dam_break.py

@@ -4 +4 @@
    Christopher Zoppou, Stephen Roberts
    Australian National University
-   
+
 """
 
@@ -15 +15 @@
 from shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\
      Dirichlet_boundary
-from pmesh2domain import pmesh_to_domain_instance
+#from pmesh2domain import pmesh_to_domain_instance
 from util import Polygon_function
 from mesh_factory import rectangular_cross
+
+
+def cut_out_region(domain):
+    """
+    To do: make better comments!
+    Deal with passing the boundary info as well
+    """
+
+    points = domain.coordinates
+    elements = domain.triangles
+    boundary = domain.boundary
+    centroid_coordinates = domain.centroid_coordinates
+    N = domain.number_of_elements
+
+    elements_in  = []
+    elements_out = []
+    for i in range(N):
+        element = elements[i]
+        #print element
+        [x,y] = centroid_coordinates[i]
+        #print x,y
+        if x>10 and y>10:
+            #print i,'Out region'
+            elements_out.append(i)
+        else:
+            #print i,'In region'
+            elements_in.append(i)
+
+    #print elements_in
+    #print elements_out
+
+    points_in = {}
+    for i in elements_in:
+        #print i
+        [v0,v1,v2] = elements[i]
+        #print v0,v1,v2
+        points_in[v0] = v0
+        points_in[v1] = v1
+        points_in[v2] = v2
+
+
+    new_index = []
+    for i,value in enumerate(points_in):
+        #print i , value
+        points_in[value] = i
+
+
+    #print points_in
+    new_elements  = []
+    for i in elements_in:
+        #print i
+        [v0,v1,v2] = elements[i]
+        #print v0,v1,v2
+
+        nv0 = points_in[v0]
+        nv1 = points_in[v1]
+        nv2 = points_in[v2]
+
+        new_elements.append([nv0,nv1,nv2])
+
+    new_points = []
+    for key,value in points_in.iteritems():
+        [x,y] = points[key]
+        new_points.append([x,y])
+
+    #print new_points
+
+    return new_points, new_elements
+
+
+
+
+
 
 
@@ -28 +101 @@
 #---------------------------------------
 # Domain
-n = 19
-m = 40
-delta_x = 5.
-delta_y = 5.
+n = 60
+m = 60
+delta_x = .5
+delta_y = .5
 lenx = delta_x*n
 leny = delta_y*m
@@ -38 +111 @@
 points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
 domain = Domain(points, elements, boundary)
+
+new_points, new_elements = cut_out_region(domain)
+domain = Domain(new_points, new_elements)
+
+
 R = Reflective_boundary(domain)
 T = Transmissive_boundary(domain)
 D = Dirichlet_boundary([h1, 0.0, 0.0])
-domain.set_boundary({'left': R, 'right': T, 'top': R, 'bottom': R})
-
-n = 2
-m = 15
-delta_x = 5.
-delta_y = 5.
-lenx = delta_x*n
-leny = delta_y*m
-origin = (90.0, 90.0)
-
-points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
-domain = Domain(points, elements, boundary)
-R = Reflective_boundary(domain)
-T = Transmissive_boundary(domain)
-D = Dirichlet_boundary([h1, 0.0, 0.0])
-domain.set_boundary({'left': T, 'right': T, 'top': R, 'bottom': R})
+domain.set_boundary({'exterior': R})
 
 
-n = 19
-m = 40
-delta_x = 5.
-delta_y = 5.
-lenx = delta_x*n
-leny = delta_y*m
-origin = (100.0, 0.0)
 
-points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
-domain = Domain(points, elements, boundary)
-R = Reflective_boundary(domain)
-T = Transmissive_boundary(domain)
-D = Dirichlet_boundary([h1, 0.0, 0.0])
-domain.set_boundary({'left': T, 'right': R, 'top': R, 'bottom': R})
+
+
+
 
 print "Number of triangles = ", len(domain)
@@ -79 +132 @@
 #---------------------------------------
 #Initial condition
-domain.set_quantity('stage', h0)
-p0 = [[0.0, 0.0], [100.0, 0.0], [100.0, 200.0], [0.0, 200.0]]
-domain.set_quantity('stage',Polygon_function([(p0,h1)]))
+p0 = [[0.0, 0.0], [10.0, 0.0], [10.0, 20.0], [0.0, 20.0]]
+domain.set_quantity('stage',Polygon_function([(p0,h1)],default = h0))
 
 
@@ -90 +142 @@
 domain.filename = 'Non_symetrical_Dam_Break_second_order'
 
+
+# Visualization smoothing
+domain.smooth=False
+domain.visualise=True
 
 #---------------------------------------
@@ -108 +164 @@
 # Friction
 domain.set_quantity('friction', 0.0)
-       
-  
+
+
 #--------------------------------------
 # Evolution
 
 yieldstep = 0.1
-finaltime = 5.0
-    
+finaltime = 15.0
+
 domain.default_order = 2
 domain.smooth = True
@@ -123 +179 @@
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
-    
+
 print 'That took %.2f seconds' %(time.time()-t0)
 
-