Changeset 488 for inundation/ga/storm_surge/pyvolution/cornell_room.py

Timestamp:

Nov 5, 2004, 4:42:19 PM (20 years ago)

Author:

ole

Message:

Played with regridding of Cornell data

File:

• inundation/ga/storm_surge/pyvolution/cornell_room.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/cornell_room.py

@@ -12 +12 @@
      Transmissive_boundary, Time_boundary, Constant_height
 
-from mesh_factory import from_polyfile
+from mesh_factory import from_polyfile, rectangular
 from Numeric import array
+from math import sqrt
+from least_squares import Interpolation
     
 
 print 'Creating domain'
-#points, triangles, values = from_polyfile('cornell_room_medres')
-points, triangles, values = from_polyfile('hires2')
+#data_points, _, data_values = from_polyfile('cornell_room_medres')
+#points, triangles, values = from_polyfile('hires2')
+data_points, _, data_values = from_polyfile('hires2')
 
 
+#Regrid onto numerically stable mesh
+#
+#Compute regular mesh based on resolution and extent of data
+data_points = array(data_points)
+pmax = max(data_points)
+pmin = min(data_points)
+
+M = len(data_points)
+
+N = int(0.8*sqrt(M))
+
+#print N
+
+mesh_points, vertices, boundary = rectangular(N, N, 
+                                              len1=pmax[0]-pmin[0],
+                                              len2=pmax[1]-pmin[1],
+                                              origin = pmin)
+                                              
+
+#Compute smooth surface on new mesh based on values from old (regrid)
+print 'Interp'
+interp = Interpolation(mesh_points, vertices, data_points, alpha=0.1)
+mesh_values = interp.fit(data_values)
+print 'Len mesh values', len(mesh_values)
+print 'Len mesh points', len(mesh_points)
+
+    
 #Create shallow water domain
-domain = Domain(points, triangles)
-    
+print 'Creating domain'
+domain = Domain(mesh_points, vertices) #, boundary)
+
 domain.check_integrity()
 domain.default_order = 2
@@ -47 +78 @@
 
 print 'Field values'
-domain.set_quantity('elevation', values)
+domain.set_quantity('elevation', mesh_values)
 domain.set_quantity('friction', manning)