Changeset 8881

Timestamp:

Jun 5, 2013, 4:08:23 PM (11 years ago)

Author:

steve

Message:

Made the computation of topography faster using numpy procedures
in numerical_varying_width.py

File:

• trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_varying_width/numerical_varying_width.py (modified) (2 diffs)

trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_varying_width/numerical_varying_width.py

@@ -12 +12 @@
 from anuga import Domain as Domain
 from math import cos
-from numpy import zeros, ones, array, interp, polyval
+from numpy import zeros, ones, array, interp, polyval, ones_like, zeros_like
+from numpy import where, logical_and
 from time import localtime, strftime, gmtime
 from scipy.interpolate import interp1d
@@ -71 +72 @@
             45.,40.,40.,30.,40.,40.,5.,40.,35.,25.,40.,40.])/2.
 
-print "Start here"
+
+depth = interp1d(XX, ZZ)
+width = interp1d(XX, WW)
+
 def bed_elevation(x,y):
-    z = 25.*ones(len(x), float)
-    for i in range(len(XX)-1):
-        for j in range(len(x)):
-            if XX[i] <= x[j] <= XX[i+1]:
-                v1 = [XX[i+1]-XX[i], WW[i+1]-WW[i]]
-                v2 = [XX[i+1]-x[j], WW[i+1]-abs(y[j])]
-                xp = v1[0]*v2[1] - v1[1]*v2[0]
-                #if 0.0 <= abs(y[j]) <= interp1d([XX[i],XX[i+1]], [WW[i],WW[i+1]])(y[j]):
-                if xp >= 0:
-                    z[j] = interp1d([XX[i],XX[i+1]], [ZZ[i],ZZ[i+1]])(x[j])
+
+    z = 25.0*ones_like(x)
+
+    wid = width(x)
+    dep = depth(x)
+
+    z = where( logical_and(y < wid, y>-wid), dep, z)
+
     return z
+
 domain.set_quantity('elevation', bed_elevation)
-print "Start now"
+
 
 #-----------------------------------------------------------------------------