Changeset 1155

Timestamp:

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

Author:

prow

Message:

sww2domian - weeding non-unique coordinates.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• data_manager.py (modified) (6 diffs)
• test_data_manager.py (modified) (10 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

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -1341 +1341 @@
 
         #Create basic mesh
-        points, vertices, boundary = rectangular(2,2)
-
+
+        yiel=0.01
+        points, vertices, boundary = rectangular(10,10)
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
@@ -1358 +1359 @@
         Bt = Transmissive_boundary(domain)
         Bd = Dirichlet_boundary([0.2,0.,0.])
-        Bw = Time_boundary(domain=domain,
+        Bw = Time_boundary(domain=domain,
                            f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
 
         domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
+
         domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
         #Initial condition
@@ -1367 +1369 @@
         elevation = domain.quantities['elevation'].vertex_values
         domain.set_quantity('stage', elevation + h)
-        #elevation = domain.get_quantity('elevation',location='unique vertices')
-        #domain.set_quantity('stage', elevation + h,location='unique vertices')
+        #elevation = domain.get_quantity('elevation')
+        #domain.set_quantity('stage', elevation + h)
 
         domain.check_integrity()
-        dir(domain)
         #Evolution
-        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
+        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
         #    domain.write_time()
             pass
@@ -1383 +1384 @@
         from data_manager import sww2domain
         from Numeric import allclose
-        import os
-
-        filename = domain.datadir+os.sep+domain.filename+'.sww'
-
-        domain2 = sww2domain(filename,fail_if_NaN=False,verbose = False)
+
+        filename = domain.datadir+'\\'+domain.filename+'.sww'
+        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
+        #points, vertices, boundary = rectangular(15,15)
+        #domain2.boundary = boundary
+        ###################
+        ##NOW TEST IT!!!
+        ###################
+
+        bits = ['xllcorner','yllcorner','vertex_coordinates']
+
+        for quantity in ['elevation']+domain.quantities_to_be_stored:
+            bits.append('quantities["%s"].get_integral()'%quantity)
+            bits.append('get_quantity("%s")'%quantity)
+
+        for bit in bits:
+        #    print 'testing that domain.'+bit+' has been restored'
+        #    print bit
+        #print 'done'
+            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
+
+        ######################################
+        #Now evolve them both, just to be sure
+        visualise = False
+        #visualise = True
+        domain.visualise = visualise
+        domain.time = 0.
+        from time import sleep
+
+        final = .1
+        domain.set_quantity('friction', 0.1)
+        domain.store = False
+        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
+
+
+        for t in domain.evolve(yieldstep = yiel, finaltime = final):
+            if visualise: sleep(1.)
+            #domain.write_time()
+            pass
+
+        domain2.smooth = False
+        domain2.visualise = visualise
+        domain2.store = False
+        domain2.default_order=2
+        domain2.set_quantity('friction', 0.1)
+        #Bed-slope and friction
+        # Boundary conditions
+        Bd2=Dirichlet_boundary([0.2,0.,0.])
+        domain2.boundary = domain.boundary
+        #print 'domain2.boundary'
+        #print domain2.boundary
+        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
+        #domain2.set_boundary({'exterior': Bd})
+
+        domain2.check_integrity()
+
+        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
+            if visualise: sleep(1.)
+            #domain2.write_time()
+            pass
+
 
         ###################
@@ -1393 +1450 @@
         ##################
 
-        bits = ['xllcorner','yllcorner','vertex_coordinates','time','starttime']
-
-        for quantity in ['elevation']+domain.quantities_to_be_stored:
+        bits = ['xllcorner','yllcorner','vertex_coordinates']
+
+        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
+            bits.append('quantities["%s"].get_integral()'%quantity)
             bits.append('get_quantity("%s")'%quantity)
 
         for bit in bits:
-        #    print 'testing that domain.'+bit+' has been restored'
+            #print bit
             assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
 
-        #print 'passed'
-
+        
 
     def test_sww2domain2(self):
@@ -1443 +1500 @@
         elevation = domain.quantities['elevation'].vertex_values
         domain.set_quantity('stage', elevation + h)
-        #elevation = domain.get_quantity('elevation',location='unique vertices')
-        #domain.set_quantity('stage', elevation + h,location='unique vertices')
 
         domain.check_integrity()
@@ -1451 +1506 @@
             pass
             #domain.write_time()
+
+        
 
         ##################################
@@ -1459 +1516 @@
         import os
 
-
         filename = domain.datadir+os.sep+domain.filename+'.sww'
 
         #Fail because NaNs are present
         try:
-            domain2 = sww2domain(filename,fail_if_NaN=True,verbose=False)
+            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
             assert True == False
         except:
             #Now import it, filling NaNs to be 0
             filler = 0
-            domain2 = sww2domain(filename,fail_if_NaN=False,NaN_filler = filler,verbose=False)
-
-        bits = ['xllcorner','yllcorner','vertex_coordinates','time','starttime']
+            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
+
+        bits = ['xllcorner','yllcorner','vertex_coordinates']
 
         for quantity in ['elevation']+domain.quantities_to_be_stored:
+            bits.append('quantities["%s"].get_integral()'%quantity)
             bits.append('get_quantity("%s")'%quantity)
 
@@ -1479 +1536 @@
         #    print 'testing that domain.'+bit+' has been restored'
             assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
-
-        #print max(max(domain2.get_quantity('xmomentum')))
-        #print min(min(domain2.get_quantity('xmomentum')))
-        #print max(max(domain2.get_quantity('ymomentum')))
-        #print min(min(domain2.get_quantity('ymomentum')))
 
         assert max(max(domain2.get_quantity('xmomentum')))==filler
@@ -1496 +1548 @@
         #os.remove(domain.datadir+'/'+domain.filename+'.sww')
 
+
+    #def test_weed(self):
+        #from data_manager import weed
+
+        #coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
+        #volumes1 = [[0,1,2],[3,4,5]]
+        #boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
+        #unique2,coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
+        #print 'coordinates2'
+        #print coordinates2
+        #print 'volumes2'
+        #print volumes2
+        #print 'boundary2'
+        #print boundary2
+
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Data_Manager,'test')    
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2asc_mis')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2asc_stage')
     runner = unittest.TextTestRunner()
     runner.run(suite)