source: anuga_validation/analytical_solutions/non_symmetrical_dam_break.py @ 7636

Last change on this file since 7636 was 7626, checked in by steve, 14 years ago

Cleaning up directory

File size: 4.0 KB
Line 
1"""Simulation of the nonsymmetrical dam dreak problem.
2
3   Copyright 2005
4   Christopher Zoppou, Stephen Roberts
5   Australian National University
6
7"""
8
9#---------------------------------------
10# Setup Path and import modules
11import sys
12from os import sep, path
13sys.path.append('..'+sep+'pyvolution')
14
15from shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\
16     Dirichlet_boundary
17#from pmesh2domain import pmesh_to_domain_instance
18from anuga.pyvolution.util import Polygon_function
19from mesh_factory import rectangular_cross
20
21
22def cut_out_region(domain):
23    """
24    To do: make better comments!
25    Deal with passing the boundary info as well
26    """
27
28    points = domain.coordinates
29    elements = domain.triangles
30    boundary = domain.boundary
31    centroid_coordinates = domain.centroid_coordinates
32    N = domain.number_of_elements
33
34    elements_in  = []
35    elements_out = []
36    for i in range(N):
37        element = elements[i]
38        #print element
39        [x,y] = centroid_coordinates[i]
40        #print x,y
41        if x>10 and y>10:
42            #print i,'Out region'
43            elements_out.append(i)
44        else:
45            #print i,'In region'
46            elements_in.append(i)
47
48    #print elements_in
49    #print elements_out
50
51    points_in = {}
52    for i in elements_in:
53        #print i
54        [v0,v1,v2] = elements[i]
55        #print v0,v1,v2
56        points_in[v0] = v0
57        points_in[v1] = v1
58        points_in[v2] = v2
59
60
61    new_index = []
62    for i,value in enumerate(points_in):
63        #print i , value
64        points_in[value] = i
65
66
67    #print points_in
68    new_elements  = []
69    for i in elements_in:
70        #print i
71        [v0,v1,v2] = elements[i]
72        #print v0,v1,v2
73
74        nv0 = points_in[v0]
75        nv1 = points_in[v1]
76        nv2 = points_in[v2]
77
78        new_elements.append([nv0,nv1,nv2])
79
80    new_points = []
81    for key,value in points_in.iteritems():
82        [x,y] = points[key]
83        new_points.append([x,y])
84
85    #print new_points
86
87    return new_points, new_elements
88
89
90
91
92
93
94
95#---------------------------------------
96# Boundary conditions
97h0 = 1.0
98h1 = 10.0
99print 'Boundary conditions'
100
101#---------------------------------------
102# Domain
103n = 100
104m = 100
105lenx = 20.0
106leny = 20.0
107delta_x = lenx/n
108delta_y = leny/m
109origin = (0.0, 0.0)
110
111points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
112domain = Domain(points, elements, boundary)
113
114new_points, new_elements = cut_out_region(domain)
115domain = Domain(new_points, new_elements)
116
117
118R = Reflective_boundary(domain)
119T = Transmissive_boundary(domain)
120D = Dirichlet_boundary([h1, 0.0, 0.0])
121domain.set_boundary({'exterior': R})
122
123
124
125
126
127
128
129print "Number of triangles = ", len(domain)
130
131
132#---------------------------------------
133#Initial condition
134p0 = [[0.0, 0.0], [10.0, 0.0], [10.0, 20.0], [0.0, 20.0]]
135domain.set_quantity('stage',Polygon_function([(p0,h1)],default = h0))
136
137
138#---------------------------------------
139# Provide file name for storing output
140domain.store = True
141domain.format = 'sww'
142domain.set_name('Non_symetrical_Dam_Break_second_order')
143
144
145# Visualization smoothing
146domain.visualise=True
147
148#---------------------------------------
149#Decide which quantities are to be stored at each timestep
150domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
151
152
153#---------------------------------------
154# Set bed elevation
155def x_slope(x,y):
156    n = x.shape[0]
157    z = 0*x
158    return z
159domain.set_quantity('elevation', x_slope)
160
161
162#----------------------------
163# Friction
164domain.set_quantity('friction', 0.0)
165
166
167#--------------------------------------
168# Evolution
169
170yieldstep = 0.1
171finaltime = 15.0
172
173domain.CFL = 0.75
174
175domain.default_order = 1
176domain.smooth = True
177
178import time
179t0 = time.time()
180for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
181    domain.write_time()
182
183print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.