source: anuga_validation/analytical solutions/Non_symmetrical_dam_break.py @ 3568

Last change on this file since 3568 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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.filename = '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.