source: anuga_work/development/classroom/spell-your-name-with-water_finished.py @ 7840

Last change on this file since 7840 was 7591, checked in by ole, 15 years ago

Added spelling your name exercise from Rowan
File size: 32.4 KB
RevLine 
[7591]1"""Simple water flow example using ANUGA
2
3This program can be used to write some text and water it become covered by water!
4"""
5######################THIS IS THE SECTION TO PUT YOUR TEXT INTO########################################################
6
7letter_number = 5 #the number of letters in the text you wish to spell with water (include spaces as letters when counting)
8
9your_text = ['A','N','U','G','A'] #put your text in here following the format shown use 'space' to put in a space
10
11#---------------------------------------------------------------------------------------------------------------------------------------------------------------
12##########Parameters - These values are intended to be experimented with to produce different results###########################################################
13#---------------------------------------------------------------------------------------------------------------------------------------------------------------
14
15depth = -50 #the elevation of the channel cut by the letters (negative number gives channels cutting into the surface,
16            #positive gives peaks rising up from the surface). Channels generally look better.
17
18letter_height = 80 #the height of each letter
19letter_width = 80 #the width of each letter
20letter_thickness = 7 #the thickness of each letter
21
22slope = -0.1        #the slope of the surface (negative number gives a downslope, positive gives an upslope)
23
24inflow_amount = 25 #controls the volume of water which flows into the model
25outflow_amount = -50 #controls the volume of water which flows out of the model
26
27amplitude = 50 #amplitude of wave (wave height m) if you want to use a wave change a boundary to Bt
28period = 10      #wave period (sec) if you want to use a wave change a boundary to Bt
29
30left_boundary = 'Bi'    #Sets up the boundary conditions for the model, Bt calls a wave input (use it to see a series of waves wash across the model)
31right_boundary = 'Bo'    # Bi calls a constant inflow amount (defined above) while Bo calls a constant outflow amount
32top_boundary = 'Br'      # Br causes any flow that contacts the boundary to be reflected back from it (the reflected flows momentum has equal magnitude but opposite direction to the incident flow)
33bottom_boundary = 'Br'
34
35simulation_speed = 1 #smaller number gives a slower, more precise simulation.
36                    #Very low numbers e.g. 0.1 will make the animation appear to run in slow motion
37simulation_length = 18*letter_number #the number of seconds that the simulation runs for
38
39#--------------------------------------------------------------------------------------------------------------------------------------------------------
40#########################################################################################################################################################
41#Please ensure that you understand what the program is doing and what you are changing before changing any of the values below this point################
42#########################################################################################################################################################
43#--------------------------------------------------------------------------------------------------------------------------------------------------------
44
45# Standard modules
46import os
47import time
48import sys
49#------------------------------------------------------------------------------
50# Import necessary modules
51#------------------------------------------------------------------------------
52from anuga.abstract_2d_finite_volumes.mesh_factory import *
53from anuga.abstract_2d_finite_volumes.domain import *
54from anuga.shallow_water import *
55from anuga.shallow_water.shallow_water_domain import *
56from anuga.shallow_water.data_manager import *
57from math import *
58from numpy import *
59from numpy import sin, cos, pi
60#--------------------------------------------------------------------------------------------------------------------------------------------------------
61# Setup computational domain
62#------------------------------------------------------------------------------
63length = (letter_width)*(letter_number) +0.2*letter_width   
64width = (letter_height)*1.3
65dx = dy = 5           # Resolution: Length of subdivisions on both axes
66
67points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
68                                               len1=length, len2=width)
69domain = Domain(points, vertices, boundary)
70domain.set_name('Spell_Your_Name_With_Water')# Output name
71domain.set_datadir('.')         #store sww output here
72
73#------------------------------------------------------------------------------
74# Setup initial conditions
75#------------------------------------------------------------------------------
76def topography(x,y):
77    """Complex topography defined by a function of vectors x and y."""
78
79    #general slope and buildings
80    z = slope*x + 0.1*(cos(x)+cos(y))+ 5
81    return z
82
83#BUILD TOPOGRAPHY from functions defining letters (above)############################
84domain.set_quantity('elevation', topography)           # elevation is a function
85
86N = len(your_text)                                      #the number of letters in the word to be spelt
87print "the length of the text is", N, "letters"
88
89print your_text
90
91for k in range(N):
92    current_letter = your_text[k]                       #the letter that is currently being plotted
93    print "the current letter is", current_letter
94   
95    #brute force approach
96    def A(x,y):
97        z = 0*x
98        N = len(x)
99        for i in range(N):
100            ymin = ((letter_height)/((letter_width)/2))*(x[i]- (k)*letter_width)
101            ymax = ((letter_height)/((letter_width)/2))*(x[i]- (k)*letter_width) + letter_thickness
102            xmin = 0
103            xmax = 0.5*letter_width
104            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
105                z[i] += depth
106        for i in range(N):
107            ymin = ((-letter_height)/((letter_width)/2))*(x[i]-(0.5*letter_width+ (k)*letter_width)) + letter_height
108            ymax = ((-letter_height)/((letter_width)/2))*(x[i]-(0.5*letter_width+ (k)*letter_width)) + letter_height + letter_thickness
109            xmin = 0.5*letter_width
110            xmax = letter_width
111            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
112                z[i] += depth
113        for i in range(N):
114            ymin = letter_height/2 - 0.5*letter_thickness
115            ymax = letter_height/2 + 0.5*letter_thickness
116            xmin = (letter_height/2)/((letter_height)/((letter_width)/2))
117            xmax = ((letter_height/2)- letter_height)/((-letter_height)/((letter_width)/2)) + 0.5*letter_width
118            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
119                z[i] += depth
120       
121        return z
122
123    def B(x,y):         
124        z = 0*x
125        N = len(x)
126        for i in range(N):
127            ymin = 0
128            ymax = letter_height
129            xmin = 0.5*letter_thickness
130            xmax = 2*letter_thickness
131            if ymin < y[i] < ymax and xmin + (k+0.05)*letter_width< x[i] < xmax+ (k+0.05)*letter_width:
132                z[i] += depth
133        for i in range(N):
134            if (x[i]-((letter_width/2)-10+ (k+0.05)*letter_width))**2+(y[i]-(0.25*letter_height))**2<((letter_height/4))**2 and letter_thickness + (k+0.05)*letter_width< x[i]:
135                z[i] += depth
136            if (x[i]-((letter_width/2)-10+ (k-0.05)*letter_width))**2+(y[i]-(0.25*letter_height))**2<(((letter_height/4)-0.6*letter_thickness))**2 and letter_thickness + (k-0.1)*letter_width< x[i]:
137                z[i] += -1*depth
138            if (x[i]-((letter_width/2)-10+ (k+0.05)*letter_width))**2+(y[i]-(0.75*letter_height))**2<((letter_height/4))**2 and letter_thickness + (k+0.05)*letter_width< x[i]:
139                z[i] += depth
140            if (x[i]-((letter_width/2)-10+ (k-0.05)*letter_width))**2+(y[i]-(0.75*letter_height))**2<(((letter_height/4)-0.6*letter_thickness))**2 and letter_thickness + (k-0.1)*letter_width< x[i]:
141                z[i] += -1*depth
142        for i in range(N):
143            ymin = (letter_height/2)-5
144            ymax = (letter_height/2)+5
145            xmin = 0.5*letter_thickness
146            xmax = 0.3*letter_width
147            if ymin < y[i] < ymax and xmin+ (k+0.05)*letter_width < x[i] < xmax+ (k+0.05)*letter_width:
148                z[i] += depth   
149        return z
150
151    def C(x,y):       
152        z = 0*x
153        N = len(x)
154        for i in range(N):
155            if (x[i]-(letter_width/2 + (k)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*(letter_height/2))**2 and 0 + (k)*letter_width < x[i] < 0.80*letter_width + (k)*letter_width:
156                z[i] += depth
157            if (x[i]-(letter_width/2 + (k)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*((letter_height/2)-letter_thickness))**2 and 0 + (k)*letter_width < x[i] < 0.80*letter_width + (k)*letter_width:
158                z[i] += -1*depth
159        return z
160
161    def D(x,y):       
162        z = 0*x
163        N = len(x)
164        for i in range(N):
165            ymin = 0
166            ymax = letter_height
167            xmin = 0.5*letter_thickness
168            xmax = 3*letter_thickness
169            if ymin < y[i] < ymax and xmin + (k)*letter_width < x[i] < xmax + (k)*letter_width:
170                z[i] += depth
171            if (x[i]-((letter_width/2)+ (k)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*(letter_height/2))**2 and 1.9*letter_thickness + (k)*letter_width < x[i] :
172                z[i] += depth
173            if (x[i]-((letter_width/2)+ (k-0.15)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*((letter_height/2)-1.3*letter_thickness))**2 and 1.9*letter_thickness + (k)*letter_width < x[i]:
174                z[i] += -1*depth
175        return z
176
177    def E(x,y):       
178        z = 0*x
179        N = len(x)
180        for i in range(N):
181            if 0 < y[i] < letter_height and 0.5*letter_thickness + (k)*letter_width < x[i] < 2*letter_thickness+ (k)*letter_width:
182                z[i] += depth
183            if (1-0.01*letter_thickness)*letter_height < y[i] < letter_height and 1.9*letter_thickness + (k)*letter_width < x[i] < 0.8*letter_width+ (k)*letter_width:
184                z[i] += depth
185            if (0.5-0.005*letter_thickness)*letter_height < y[i] < (0.5+0.005*letter_thickness)*letter_height and 1.9*letter_thickness + (k)*letter_width< x[i] < 0.8*letter_width+ (k)*letter_width:
186                z[i] += depth
187            if 0 < y[i] < 0.01*letter_thickness*letter_height and 1.9*letter_thickness + (k)*letter_width< x[i] < 0.8*letter_width+ (k)*letter_width:
188                z[i] += depth
189        return z
190
191    def F(x,y):       
192        z = 0*x
193        N = len(x)
194        for i in range(N):
195            if 0 < y[i] < letter_height and 0.5*letter_thickness + (k)*letter_width< x[i] < 2*letter_thickness+ (k)*letter_width:
196                z[i] += depth
197            if (1-0.01*letter_thickness)*letter_height < y[i] < letter_height and 1.9*letter_thickness + (k)*letter_width< x[i] < 0.8*letter_width+ (k)*letter_width:
198                z[i] += depth
199            if (0.5-0.005*letter_thickness)*letter_height < y[i] < (0.5+0.005*letter_thickness)*letter_height and 1.9*letter_thickness + (k)*letter_width< x[i] < 0.8*letter_width+ (k)*letter_width:
200                z[i] += depth
201        return z
202
203    def G(x,y):       
204        z = 0*x
205        N = len(x)
206        for i in range(N):
207            if (x[i]-(letter_width/2+ (k)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*(letter_height/2))**2 and 0 + (k)*letter_width< x[i] < 0.80*letter_width+ (k)*letter_width:
208                z[i] += depth
209            if (x[i]-(letter_width/2+ (k)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*((letter_height/2)-letter_thickness))**2 and 0 + (k)*letter_width< x[i] < 0.80*letter_width+ (k)*letter_width:
210                z[i] += -1*depth
211            if (0.4-0.005*letter_thickness)*letter_height < y[i] < (0.4+0.005*letter_thickness)*letter_height  and 0.5*letter_width + (k)*letter_width< x[i] < 0.95*letter_width+ (k)*letter_width:
212                z[i] += depth
213            if 0 < y[i] < 0.4*letter_height and (0.8-0.01*letter_thickness)*letter_width + (k)*letter_width< x[i] < 0.8*letter_width+ (k)*letter_width:
214                z[i] += depth
215        return z
216
217    def H(x,y):       
218        z = 0*x
219        N = len(x)
220        for i in range(N):
221            if 0 < y[i] < letter_height and 0.5*letter_thickness + (k)*letter_width< x[i] < 2*letter_thickness+ (k)*letter_width:
222                z[i] += depth
223            if 0 < y[i] < letter_height and (letter_width - 2*letter_thickness) + (k)*letter_width< x[i] < (letter_width -0.5*letter_thickness)+ (k)*letter_width:
224                z[i] += depth
225            if (0.5-0.005*letter_thickness)*letter_height < y[i] < (0.5+0.005*letter_thickness)*letter_height and 1.9*letter_thickness + (k)*letter_width< x[i] < (letter_width-1.9*letter_thickness)+ (k)*letter_width:
226                z[i] += depth
227        return z
228
229    def I(x,y):       
230        z = 0*x
231        N = len(x)
232        for i in range(N):
233            if 0 < y[i] < letter_height and (0.5-0.01*letter_thickness)*letter_width + (k)*letter_width< x[i] < (0.5+0.01*letter_thickness)*letter_width+ (k)*letter_width:
234                z[i] += depth
235            if (letter_height - 0.85*letter_thickness) < y[i] < letter_height and 0.3*letter_width + (k)*letter_width< x[i] < 0.7*letter_width+ (k)*letter_width:
236                z[i] += depth
237            if 0 < y[i] < 0.85*letter_thickness and 0.3*letter_width + (k)*letter_width< x[i] < 0.7*letter_width+ (k)*letter_width:
238                z[i] += depth
239        return z
240
241    def J(x,y):       
242        z = 0*x
243        N = len(x)
244        for i in range(N):
245            if y[i]< -(((0.9*(letter_height/2))**2 - (x[i]- (0.5*letter_width+ (k)*letter_width))**2)**0.5)+ 0.6*letter_height and 0<y[i]<0.3*letter_height:
246                z[i] += depth
247            if y[i]< -(((0.9*(letter_height/2))**2 - (x[i]- (0.5*letter_width+ (k)*letter_width))**2)**0.5)+ 0.5*letter_height and 0<y[i]<0.3*letter_height:
248                z[i] += -depth
249            if 0.2*letter_height < y[i] < letter_height and (0.83-0.01*letter_thickness)*letter_width + (k)*letter_width< x[i] < (0.83+0.01*letter_thickness)*letter_width+ (k)*letter_width:
250                z[i] += depth
251        return z
252
253    def K(x,y):       
254        z = 0*x
255        N = len(x)
256        for i in range(N):
257            if 0 < y[i] < letter_height and 0.5*letter_thickness + (k+0.2)*letter_width < x[i] < 2*letter_thickness + (k+0.2)*letter_width:
258                z[i] += depth
259        for i in range(N):
260            ymin = ((0.5*letter_height)/((letter_width)))*(x[i]-(1.9*letter_thickness+ (k+0.2)*letter_width)) + 0.5*letter_height
261            ymax = ((0.5*letter_height)/((letter_width)))*(x[i]-(1.9*letter_thickness+ (k+0.2)*letter_width)) + 0.5*letter_height + letter_thickness
262            xmin = 1.9*letter_thickness
263            xmax = letter_width
264            if ymin < y[i] < ymax and xmin + (k+0.2)*letter_width< x[i] < xmax+ (k+0.2)*letter_width:
265                z[i] += depth
266        for i in range(N):
267            ymin = (-(0.5*letter_height)/((letter_width)))*(x[i]-(1.9*letter_thickness+ (k+0.2)*letter_width)) + 0.5*letter_height - letter_thickness
268            ymax = (-(0.5*letter_height)/((letter_width)))*(x[i]-(1.9*letter_thickness+ (k+0.2)*letter_width)) + 0.5*letter_height
269            xmin = 1.9*letter_thickness
270            xmax = letter_width
271            if ymin < y[i] < ymax and xmin + (k+0.2)*letter_width< x[i] < xmax+ (k+0.2)*letter_width:
272                z[i] += depth   
273        return z
274
275    def L(x,y):       
276        z = 0*x
277        N = len(x)
278        for i in range(N):
279            if 0 < y[i] < letter_height and 0.5*letter_thickness + (k)*letter_width < x[i] < 2*letter_thickness+ (k)*letter_width:
280                z[i] += depth
281            if 0 < y[i] < 0.01*letter_thickness*letter_height and 1.9*letter_thickness + (k)*letter_width < x[i] < 0.95*letter_width+ (k)*letter_width:
282                z[i] += depth
283        return z
284
285    def M(x,y):
286        z = 0*x
287        N = len(x)
288        for i in range(N):
289            if 0 < y[i] < letter_height and 0.5*letter_thickness + (k)*letter_width < x[i] < 2*letter_thickness+ (k)*letter_width:
290                z[i] += depth
291            if 0 < y[i] < letter_height and (letter_width - 2*letter_thickness)+ (k)*letter_width < x[i] < (letter_width -0.5*letter_thickness) + (k)*letter_width:
292                z[i] += depth
293        for i in range(N):
294            ymin = (-(0.7*letter_height)/((0.5*letter_width)))*(x[i]-(1.9*letter_thickness+ (k)*letter_width)) + letter_height - 2*letter_thickness
295            ymax = (-(0.7*letter_height)/((0.5*letter_width)))*(x[i]-(1.9*letter_thickness+ (k)*letter_width)) + letter_height
296            xmin = 1.9*letter_thickness
297            xmax = 0.5*letter_width
298            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
299                z[i] += depth
300        for i in range(N):
301            ymin = (0.7*(letter_height)/(0.5*(letter_width)))*(x[i]-(0.5*letter_width+ (k)*letter_width)) + (-(0.7*letter_height)/((0.5*letter_width)))*((0.5*letter_width)-(1.9*letter_thickness)) + letter_height - 2*letter_thickness
302            ymax = (0.7*(letter_height)/(0.5*(letter_width)))*(x[i]-(0.5*letter_width+ (k)*letter_width)) + (-(0.7*letter_height)/((0.5*letter_width)))*((0.5*letter_width)-(1.9*letter_thickness)) + letter_height
303            xmin = 0.5*letter_width
304            xmax = letter_width - 1.9*letter_thickness
305            if ymin < y[i] < ymax and xmin + (k)*letter_width < x[i] < xmax+ (k)*letter_width:
306                z[i] += depth   
307        return z
308
309    def N(x,y):
310        z = 0*x
311        p = len(x)
312        for i in range(p):
313            if 0 < y[i] < letter_height and 0.5*letter_thickness + (k)*letter_width< x[i] < 2*letter_thickness+ (k)*letter_width:
314                z[i] += depth
315            if 0 < y[i] < letter_height and (letter_width - 2*letter_thickness)+ (k)*letter_width < x[i] < (letter_width - 0.5*letter_thickness)+ (k)*letter_width:
316                z[i] += depth
317        for i in range(p):
318            ymin = (-(letter_height)/((letter_width)-1.9*letter_thickness))*(x[i]-(1.9*letter_thickness+ (k)*letter_width)) + letter_height - 2*letter_thickness
319            ymax = (-(letter_height)/((letter_width)-1.9*letter_thickness))*(x[i]-(1.9*letter_thickness+ (k)*letter_width)) + letter_height
320            xmin = 1.9*letter_thickness
321            xmax = letter_width - 1.9*letter_thickness
322            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
323                z[i] += depth
324        return z
325
326    def O(x,y):       
327        z = 0*x
328        N = len(x)
329        for i in range(N):
330            if (x[i]-(letter_width/2+ (k)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*(letter_height/2))**2 and 0+ (k)*letter_width < x[i] < letter_width+ (k)*letter_width:
331                z[i] += depth
332            if (x[i]-(letter_width/2+ (k)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*((letter_height/2)-letter_thickness))**2 and 0+ (k)*letter_width < x[i] < letter_width+ (k)*letter_width:
333                z[i] += -1*depth
334        return z
335
336    def P(x,y):         
337        z = 0*x
338        N = len(x)
339        for i in range(N):
340            ymin = 0
341            ymax = letter_height
342            xmin = 0.5*letter_thickness
343            xmax = 2*letter_thickness
344            if ymin < y[i] < ymax and xmin + (k+0.2)*letter_width< x[i] < xmax+ (k+0.2)*letter_width:
345                z[i] += depth
346        for i in range(N):
347            if (x[i]-((letter_width/2)-10+ (k+0.2)*letter_width))**2+(y[i]-(0.75*letter_height))**2<((letter_height/4))**2 and letter_thickness + (k+0.2)*letter_width< x[i]:
348                z[i] += depth
349            if (x[i]-((letter_width/2)-10+ (k+0.2)*letter_width))**2+(y[i]-(0.75*letter_height))**2<(((letter_height/4)-letter_thickness))**2 and letter_thickness + (k+0.2)*letter_width< x[i]:
350                z[i] += -1*depth
351        for i in range(N):
352            ymin = (letter_height/2)-5
353            ymax = (letter_height/2)+5
354            xmin = 0
355            xmax = 0.3*letter_width
356            if ymin < y[i] < ymax and xmin + (k+0.2)*letter_width< x[i] < xmax+ (k+0.2)*letter_width:
357                z[i] += depth   
358        return z
359
360    def Q(x,y):       
361        z = 0*x
362        N = len(x)
363        for i in range(N):
364            if (x[i]-(letter_width/2+ (k)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*(letter_height/2))**2 and 0 + (k)*letter_width< x[i] < letter_width+ (k)*letter_width:
365                z[i] += depth
366            if (x[i]-(letter_width/2+ (k)*letter_width))**2+(y[i]-(letter_height/2))**2<(0.9*((letter_height/2)-letter_thickness))**2 and 0 + (k)*letter_width< x[i] < letter_width+ (k)*letter_width:
367                z[i] += -1*depth
368        for i in range(N):
369            ymin = (-(letter_height)/((letter_width)-1.9*letter_thickness))*(x[i]-(1.9*letter_thickness+ (k)*letter_width)) + letter_height - 2*letter_thickness
370            ymax = (-(letter_height)/((letter_width)-1.9*letter_thickness))*(x[i]-(1.9*letter_thickness+ (k)*letter_width)) + letter_height
371            xmin = 0.55*letter_width
372            xmax = 0.95*letter_width
373            if ymin < y[i] < ymax and xmin + (k)*letter_width < x[i] < xmax + (k)*letter_width:
374                z[i] += depth
375        return z
376
377    def R(x,y):         
378        z = 0*x
379        N = len(x)
380        for i in range(N):
381            ymin = 0
382            ymax = letter_height
383            xmin = 0.5*letter_thickness
384            xmax = 2*letter_thickness
385            if ymin < y[i] < ymax and xmin + (k)*letter_width < x[i] < xmax+ (k)*letter_width:
386                z[i] += depth
387        for i in range(N):
388            if (x[i]-((letter_width/2)-10+ (k)*letter_width))**2+(y[i]-(0.75*letter_height))**2<((letter_height/4))**2 and letter_thickness + (k)*letter_width< x[i]:
389                z[i] += depth
390            if (x[i]-((letter_width/2)-10+ (k)*letter_width))**2+(y[i]-(0.75*letter_height))**2<(((letter_height/4)-letter_thickness))**2 and letter_thickness + (k)*letter_width < x[i]:
391                z[i] += -1*depth
392        for i in range(N):
393            ymin = (letter_height/2)-5
394            ymax = (letter_height/2)+5
395            xmin = 0.5*letter_thickness
396            xmax = 0.3*letter_width
397            if ymin < y[i] < ymax and xmin + (k)*letter_width < x[i] < xmax+ (k)*letter_width:
398                z[i] += depth
399        for i in range(N):
400            ymin = (-(letter_height)/((letter_width)))*(x[i]-(1*letter_thickness+ (k)*letter_width)) + letter_height - 2*letter_thickness
401            ymax = (-(letter_height)/((letter_width)))*(x[i]-(1*letter_thickness+ (k)*letter_width)) + letter_height
402            xmin = 0.5*letter_width
403            xmax = 0.95*letter_width
404            if ymin < y[i] < ymax and xmin + (k)*letter_width < x[i] < xmax + (k)*letter_width:
405                z[i] += depth
406        return z
407
408    def S(x,y):       
409        z = 0*x
410        N = len(x)
411        for i in range(N):
412            if 0.3*letter_width*sin((2*pi*y[i]/(0.7*letter_height))-0.1*letter_height) +0.5*letter_width+ (k)*letter_width< x[i] < 0.3*letter_width*sin((2*pi*y[i]/(0.7*letter_height))-0.1*letter_height)+(0.5*letter_width+2*letter_thickness)+ (k)*letter_width:
413                z[i] += depth
414        return z
415
416    def T(x,y):       
417        z = 0*x
418        N = len(x)
419        for i in range(N):
420            if 0 < y[i] < letter_height and (0.5-0.005*letter_thickness)*letter_width + (k)*letter_width < x[i] < (0.5+0.005*letter_thickness)*letter_width + (k)*letter_width:
421                z[i] += depth
422            if letter_height - letter_thickness < y[i] < letter_height and 0.08*letter_width + (k)*letter_width < x[i] < 0.95*letter_width + (k)*letter_width:
423                z[i] += depth
424        return z
425
426    def U(x,y):       
427        z = 0*x
428        N = len(x)
429        for i in range(N):
430            if y[i]< -(((0.9*(letter_height/2))**2 - (x[i]- (0.5*letter_width+ (k)*letter_width))**2)**0.5)+ 0.6*letter_height and 0<y[i]<0.3*letter_height:
431                z[i] += depth
432            if y[i]< -(((0.9*(letter_height/2))**2 - (x[i]- (0.5*letter_width+ (k)*letter_width))**2)**0.5)+ 0.5*letter_height and 0<y[i]<0.3*letter_height:
433                z[i] += -depth
434            if 0.2*letter_height < y[i] < letter_height and (0.83-0.01*letter_thickness)*letter_width + (k)*letter_width< x[i] < (0.83+0.01*letter_thickness)*letter_width+ (k)*letter_width:
435                z[i] += depth
436            if 0.2*letter_height < y[i] < letter_height and (0.18-0.01*letter_thickness)*letter_width + (k)*letter_width< x[i] < (0.18+0.01*letter_thickness)*letter_width+ (k)*letter_width:
437                z[i] += depth
438        return z
439
440    def V(x,y):       
441        z = 0*x
442        N = len(x)
443        for i in range(N):
444            ymin = (-(letter_height)/((0.5*letter_width)))*(x[i]- (k)*letter_width) + (letter_height)
445            ymax = (-(letter_height)/((0.5*letter_width)))*(x[i]- (k)*letter_width) + (letter_height + 2*letter_thickness)
446            xmin = 0
447            xmax = 0.5*letter_width
448            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
449                z[i] += depth
450        for i in range(N):
451            ymin = ((letter_height)/((0.5*letter_width)))*(x[i]-(0.5*letter_width+ (k)*letter_width)) 
452            ymax = ((letter_height)/((0.5*letter_width)))*(x[i]-(0.5*letter_width+ (k)*letter_width))+ 2*letter_thickness
453            xmin = 0.5*letter_width
454            xmax = letter_width
455            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
456                z[i] += depth   
457        return z
458
459    def W(x,y):       
460        z = 0*x
461        N = len(x)
462        for i in range(N):
463            ymin = (-(letter_height)/((0.25*letter_width)))*(x[i]- (k)*letter_width) + (letter_height-letter_thickness)
464            ymax = (-(letter_height)/((0.25*letter_width)))*(x[i]- (k)*letter_width) + (letter_height + 2*letter_thickness)
465            xmin = 0
466            xmax = 0.25*letter_width
467            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
468                z[i] += depth
469        for i in range(N):
470            ymin = ((0.5*letter_height)/((0.25*letter_width)))*(x[i]-(0.25*letter_width+ (k)*letter_width)) 
471            ymax = ((0.5*letter_height)/((0.25*letter_width)))*(x[i]-(0.25*letter_width+ (k)*letter_width))+ 2*letter_thickness
472            xmin = 0.25*letter_width
473            xmax = 0.5*letter_width
474            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
475                z[i] += depth
476        for i in range(N):
477            ymin = (-(0.5*letter_height)/((0.25*letter_width)))*(x[i]-(0.5*letter_width+ (k)*letter_width)) + (0.5*letter_height)
478            ymax = (-(0.5*letter_height)/((0.25*letter_width)))*(x[i]-(0.5*letter_width+ (k)*letter_width)) + (0.5*letter_height + 2*letter_thickness)
479            xmin = 0.5*letter_width
480            xmax = 0.75*letter_width
481            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
482                z[i] += depth
483        for i in range(N):
484            ymin = ((letter_height)/((0.25*letter_width)))*(x[i]-(0.75*letter_width+ (k)*letter_width)) -letter_thickness
485            ymax = ((letter_height)/((0.25*letter_width)))*(x[i]-(0.75*letter_width+ (k)*letter_width))+ 2*letter_thickness
486            xmin = 0.75*letter_width
487            xmax = letter_width
488            if ymin < y[i] < ymax and xmin + (k)*letter_width< x[i] < xmax+ (k)*letter_width:
489                z[i] += depth   
490        return z
491
492    def X(x,y):       
493        z = 0*x
494        N = len(x)
495        for i in range(N):
496            ymin = (-(letter_height)/((letter_width)))*(x[i]- (k+0.2)*letter_width) + (letter_height-letter_thickness)
497            ymax = (-(letter_height)/((letter_width)))*(x[i]- (k+0.2)*letter_width) + (letter_height + letter_thickness)
498            xmin = 0
499            xmax = letter_width
500            if ymin < y[i] < ymax and xmin + (k+0.2)*letter_width< x[i] < xmax+ (k+0.2)*letter_width:
501                z[i] += depth
502        for i in range(N):
503            ymin = ((letter_height)/((letter_width)))*(x[i]- (k+0.2)*letter_width) -letter_thickness
504            ymax = ((letter_height)/((letter_width)))*(x[i]- (k+0.2)*letter_width)+ letter_thickness
505            xmin = 0
506            xmax = letter_width
507            if ymin < y[i] < ymax and xmin + (k+0.2)*letter_width< x[i] < xmax+ (k+0.2)*letter_width:
508                z[i] += depth   
509        return z
510
511    def Y(x,y):       
512        z = 0*x
513        N = len(x)
514        for i in range(N):
515            if 0 < y[i] < 0.5*letter_height and (0.5-0.005*letter_thickness)*letter_width + (k)*letter_width < x[i] < (0.5+0.005*letter_thickness)*letter_width + (k)*letter_width:
516                z[i] += depth
517        for i in range(N):
518            ymin = (-(letter_height)/((letter_width)))*(x[i]- (k)*letter_width) + (letter_height-0.6*letter_thickness)
519            ymax = (-(letter_height)/((letter_width)))*(x[i]- (k)*letter_width) + (letter_height + 0.6*letter_thickness)
520            xmin = 0
521            xmax = 0.5*letter_width
522            if ymin < y[i] < ymax and xmin + (k)*letter_width < x[i] < xmax + (k)*letter_width:
523                z[i] += depth
524        for i in range(N):
525            ymin = ((letter_height)/((letter_width)))*(x[i]- (k)*letter_width) -0.6*letter_thickness
526            ymax = ((letter_height)/((letter_width)))*(x[i]- (k)*letter_width)+ 0.6*letter_thickness
527            xmin = 0.5*letter_width
528            xmax = letter_width
529            if ymin < y[i] < ymax and xmin + (k)*letter_width < x[i] < xmax + (k)*letter_width:
530                z[i] += depth   
531        return z
532
533    def Z(x,y):       
534        z = 0*x
535        N = len(x)
536        for i in range(N):
537            ymin = ((letter_height)/((letter_width)))*(x[i]- (k)*letter_width) -0.6*letter_thickness
538            ymax = ((letter_height)/((letter_width)))*(x[i]- (k)*letter_width)+ 0.6*letter_thickness
539            xmin = 0
540            xmax = letter_width
541            if ymin < y[i] < ymax and xmin + (k)*letter_width < x[i] < xmax + (k)*letter_width:
542                z[i] += depth
543            if letter_height - letter_thickness < y[i] < letter_height and 0 + (k)*letter_width < x[i] < letter_width + (k)*letter_width:
544                z[i] += depth
545            if 0 < y[i] < letter_thickness and 0 + (k)*letter_width < x[i] < letter_width + (k)*letter_width:
546                z[i] += depth
547        return z
548
549    def space(x,y):
550        z=0*x
551        return z
552
553    def translation(x,y):                               #shifts each letter along one place so they are not plotted on top of each other
554        x = x + (i-1)*letter_width
555        return x
556        y = y
557        return y
558        z = 50*sin(x)
559       
560    current_letter_function = eval(current_letter)      #function of vectors x and y whose output forms the shape of a letter
561    domain.add_quantity('elevation', current_letter_function)               #which is input here to be used as elevation.
562
563###############################################################################################################################################
564
565domain.set_quantity('friction', 0.01)                  # Constant friction
566domain.set_quantity('stage', expression='elevation')   # Dry initial condition
567domain.add_quantity('stage', -1000)
568
569
570#thunderstorm = Rainfall(domain,rate=100,center=(50,50),radius=5)       #remove comments to add rainfall input to a specified position on the grid
571#domain.forcing_terms.append(thunderstorm)
572
573
574#------------------------------------------------------------------------------
575# Setup boundary conditions
576#------------------------------------------------------------------------------
577Bi = Dirichlet_boundary([inflow_amount, 0, 0])          # Inflow
578Br = Reflective_boundary(domain)              # Solid reflective wall
579Bo = Dirichlet_boundary([outflow_amount, 0, 0])           # Outflow
580
581def wave(t):
582
583    A = amplitude # Amplitude [m] (Wave height)
584    T = period  # Wave period [s]
585
586    if t < 3000:
587        return [A*sin(2*pi*t/T) + 1, 0, 0] 
588    else:
589        return [0.0, 0, 0]
590
591Bt = Time_boundary(domain, f=wave)
592
593left_b = eval(left_boundary)
594right_b = eval(right_boundary)
595top_b = eval(top_boundary)
596bottom_b = eval(bottom_boundary)
597
598domain.set_boundary({'left': left_b, 'right': right_b, 'top': top_b, 'bottom': bottom_b})
599
600#------------------------------------------------------------------------------
601# Evolve system through time
602#------------------------------------------------------------------------------
603t0 = time.time()
604for t in domain.evolve(yieldstep=simulation_speed, finaltime=simulation_length):
605    print domain.timestepping_statistics()
606print 'Computation took %.2f seconds' % (time.time()-t0)
607   
Note: See TracBrowser for help on using the repository browser.