Changeset 1155 for inundation/ga/storm_surge/pyvolution/data_manager.py

Timestamp:

Mar 29, 2005, 1:09:46 PM (20 years ago)

Author:

prow

Message:

sww2domian - weeding non-unique coordinates.

File:

• inundation/ga/storm_surge/pyvolution/data_manager.py (modified) (6 diffs)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -1835 +1835 @@
 
 
-def sww2domain(filename,t=None,fail_if_NaN=True,NaN_filler=0,verbose = True):
+def sww2domain(filename,boundary=None,t=None,\
+               fail_if_NaN=True,NaN_filler=0\
+               ,verbose = True,very_verbose = False):
     """
     Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
 
-    If the sww has stages, but not 
+    Boundary is not recommended if domian.smooth is not selected, as it       
+    uses unique coordinates, but not unique boundaries. This means that 
+    the boundary file will not be compatable with the coordiantes, and will
+    give a different final boundary, or crash.
     """
     NaN=9.969209968386869e+036
@@ -1845 +1850 @@
 
     from Scientific.IO.NetCDF import NetCDFFile
-    from domain import Domain
-    from Numeric import asarray, transpose
+    from shallow_water import Domain
+    from Numeric import asarray, transpose, resize
 
     if verbose: print 'Reading from ', filename
@@ -1891 +1896 @@
 
     if verbose: print '    building domain'
-    domain = Domain(coordinates, volumes,\
-                    conserved_quantities = conserved_quantities,\
-                    other_quantities = other_quantities,zone=zone,\
-                    xllcorner=xllcorner, yllcorner=yllcorner)
-    domain.starttime=starttime
-    domain.time=t
+#    From domain.Domain:
+#    domain = Domain(coordinates, volumes,\
+#                    conserved_quantities = conserved_quantities,\
+#                    other_quantities = other_quantities,zone=zone,\
+#                    xllcorner=xllcorner, yllcorner=yllcorner)
+
+#   From shallow_water.Domain:
+    coordinates=coordinates.tolist()
+    volumes=volumes.tolist()
+    #FIXME:should this be in mesh?(peter row)
+    unique,coordinates,volumes,boundary=weed(coordinates,volumes,boundary)
+
+    domain = Domain(coordinates, volumes, boundary)
+
+    if not boundary is None:
+        domain.boundary = boundary
+    domain.zone=zone
+    domain.xllcorner=float(xllcorner)
+    domain.yllcorner=float(yllcorner)
+
+    domain.starttime=float(starttime)+float(t)
+    domain.time=0.0
+
     for quantity in other_quantities:
         try:
@@ -1918 +1940 @@
                 data = (X<>NaN)
                 X = (X*data)+(data==0)*NaN_filler
+        if unique:
+            X = resize(X,(len(X)/3,3))
         domain.set_quantity(quantity,X)
 #
@@ -1941 +1965 @@
                 data = (X<>NaN)
                 X = (X*data)+(data==0)*NaN_filler
+        if unique:
+            X = resize(X,(X.shape[0]/3,3))
         domain.set_quantity(quantity,X)
     fid.close()
@@ -1985 +2011 @@
 
 
+def weed(coordinates,volumes,boundary = None):
+    if type(coordinates)=='array':
+        coordinates = coordinates.tolist()
+    if type(volumes)=='array':
+        volumes = volumes.tolist()
+
+    unique = False
+    point_dict = {}
+    same_point = {}
+    for i in range(len(coordinates)):
+        point = tuple(coordinates[i])
+        if point_dict.has_key(point):
+            unique = True
+            same_point[i]=point
+            #to change all point i references to point j
+        else: 
+            point_dict[point]=i
+            same_point[i]=point
+
+    coordinates = []
+    i = 0
+    for point in point_dict.keys():
+        point = tuple(point)
+        coordinates.append(list(point))
+        point_dict[point]=i
+        i+=1
+
+
+    for volume in volumes:
+        for i in range(len(volume)):
+            index = volume[i]
+            if index>-1:
+                volume[i]=point_dict[same_point[index]]
+
+    new_boundary = {}
+    if not boundary is None:
+        for segment in boundary.keys():
+            point0 = point_dict[same_point[segment[0]]]
+            point1 = point_dict[same_point[segment[1]]]
+            label = boundary[segment]
+            #FIXME should the bounday attributes be concaterated
+            #('exterior, pond') or replaced ('pond')(peter row)
+
+            if new_boundary.has_key((point0,point1)):
+                new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\
+                                              #+','+label
+
+            elif new_boundary.has_key((point1,point0)):
+                new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\
+                                              #+','+label
+            else: new_boundary[(point0,point1)]=label
+
+        boundary = new_boundary               
+
+    return unique,coordinates,volumes,boundary
+        
+        
 
 #OBSOLETE STUFF