Changeset 612 for inundation/ga/storm_surge/pyvolution/wind_example_variable.py

Timestamp:

Nov 22, 2004, 5:29:10 PM (20 years ago)

Author:

ole

Message:

Work towards wind_field - it now accepts vect values fuctions or scalars

File:

• inundation/ga/storm_surge/pyvolution/wind_example_variable.py (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/wind_example_variable.py

@@ -13 +13 @@
 #Create basic mesh
 N = 80
-len = 100
-points, vertices, boundary = rectangular(N, N, len, len)
+length = 100
+points, vertices, boundary = rectangular(N, N, length, length)
 
 #Create shallow water domain
@@ -35 +35 @@
     Low speeds at center and edges
     """
+
+    from math import sqrt, exp, cos, pi
+
+    N = len(x)
+    s = 0*x  #New array
+
+    c = (length/2, length/2)
     
-    from math import sqrt, exp, cos, pi
-    
-    c = (len1/2, len2/2)
-    r = sqrt((x - c[0])**2 + (y - c[1])**2)
+    for k in range(N):
+        r = sqrt((x[k] - c[0])**2 + (y[k] - c[1])**2)
 
-    r /= max(len1,len2)
+        r /= max(length,length)
 
-    factor = exp( -(r-0.15)**2 )
+        factor = exp( -(r-0.15)**2 )
         
-    #print x,y,factor
+        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
 
-    return 4000 * factor * (cos(t*2*pi/150) + 2)
+    return s
 
 
@@ -55 +60 @@
     from math import sqrt, atan, pi
 
-    c = (len1/2, len2/2)
-    r = sqrt((x - c[0])**2 + (y - c[1])**2)    
+    N = len(x)
+    a = 0*x  #New array
+
+    c = (length/2, length/2)
+    for k in range(N):
+        r = sqrt((x[k] - c[0])**2 + (y[k] - c[1])**2)
+
+        xx = (x[k] - c[0])/length
+        yy = (y[k] - c[1])/length
+
+        angle = atan(yy/xx)
+
+        if xx < 0:
+            angle+=pi #atan in ]-pi/2; pi/2[
+
+        #Take normal direction
+        angle -= pi/2
     
-    xx = (x - c[0])/len1
-    yy = (y - c[1])/len2
+        #Ensure positive radians    
+        if angle < 0:
+            angle += 2*pi
 
-    angle = atan(yy/xx)
-
-    if xx < 0:
-        angle+=pi #atan in ]-pi/2; pi/2[
-
-    #Take normal direction
-    angle -= pi/2
-    
-    #Ensure positive radians    
-    if angle < 0:
-        angle += 2*pi
-
-    return angle/pi*180
+        a[k] = angle/pi*180
+        
+    return a
 
     
@@ -82 +93 @@
 def gust(x,y,t): 
     from math import sin, pi
+    from Numeric import zeros, ones, Float
+
+    N = len(x)
 
     tt = sin(2*pi*t/200)
 
     if tt > 0.9:
-        return 6000*tt
+        return 6000*tt*ones(N, Float)
     else:
-        return 0
+        return zeros(N, Float)
     
 domain.forcing_terms.append(Wind_stress(gust, 25))