Changeset 1018
- Timestamp:
- Mar 7, 2005, 9:32:07 AM (20 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 17 edited
-
data_manager.py (modified) (98 diffs)
-
test_advection.py (modified) (10 diffs)
-
test_all.py (modified) (4 diffs)
-
test_cg_solve.py (modified) (11 diffs)
-
test_data_manager.py (modified) (50 diffs)
-
test_domain.py (modified) (26 diffs)
-
test_general_mesh.py (modified) (4 diffs)
-
test_generic_boundary_conditions.py (modified) (14 diffs)
-
test_interpolate_sww.py (modified) (14 diffs)
-
test_least_squares.py (modified) (81 diffs)
-
test_mesh.py (modified) (51 diffs)
-
test_pmesh2domain.py (modified) (8 diffs)
-
test_quad.py (modified) (9 diffs)
-
test_quantity.py (modified) (63 diffs)
-
test_region.py (modified) (13 diffs)
-
test_sparse.py (modified) (14 diffs)
-
test_util.py (modified) (42 diffs)
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/data_manager.py
r995 r1018 5 5 2: Variable data: Conserved quantities. Stored once per timestep. 6 6 7 All data is assumed to reside at vertex locations. 7 All data is assumed to reside at vertex locations. 8 8 9 9 … … 38 38 PTS + TSH -> TSH with elevation: Least squares fit 39 39 40 TSH -> SWW: Conversion of TSH to sww viewable using Swollen 40 TSH -> SWW: Conversion of TSH to sww viewable using Swollen 41 41 42 42 TSH + Boundary SWW -> SWW: SImluation using pyvolution … … 45 45 """ 46 46 47 from Numeric import concatenate 47 from Numeric import concatenate 48 48 49 49 … … 56 56 s = s.replace('(', '') 57 57 s = s.replace(')', '') 58 s = s.replace('__', '_') 59 58 s = s.replace('__', '_') 59 60 60 return s 61 61 … … 63 63 def check_dir(path, verbose=None): 64 64 """Check that specified path exists. 65 If path does not exist it will be created if possible 65 If path does not exist it will be created if possible 66 66 67 67 USAGE: … … 74 74 RETURN VALUE: 75 75 Verified path including trailing separator 76 76 77 77 """ 78 78 … … 86 86 87 87 88 if path[-1] != os.sep: 88 if path[-1] != os.sep: 89 89 path = path + os.sep # Add separator for directories 90 90 … … 100 100 else: 101 101 pass # FIXME: What about acces rights under Windows? 102 102 103 103 if verbose: print 'MESSAGE: Directory', path, 'created.' 104 104 105 105 except: 106 106 print 'WARNING: Directory', path, 'could not be created.' … … 109 109 else: 110 110 path = 'C:' 111 111 112 112 print 'Using directory %s instead' %path 113 113 114 114 return(path) 115 115 116 116 117 117 118 118 def del_dir(path): … … 128 128 129 129 if os.path.isdir(X) and not os.path.islink(X): 130 del_dir(X) 130 del_dir(X) 131 131 else: 132 132 try: … … 138 138 139 139 140 140 141 141 def create_filename(datadir, filename, format, size=None, time=None): 142 142 143 143 import os 144 144 #from config import data_dir 145 145 146 146 FN = check_dir(datadir) + filename 147 148 if size is not None: 147 148 if size is not None: 149 149 FN += '_size%d' %size 150 150 151 151 if time is not None: 152 152 FN += '_time%.2f' %time 153 153 154 154 FN += '.' + format 155 155 return FN 156 156 157 157 158 158 def get_files(datadir, filename, format, size): … … 164 164 import os 165 165 #from config import data_dir 166 166 167 167 dir = check_dir(datadir) 168 168 … … 177 177 """ 178 178 179 179 180 180 def __init__(self, domain, extension, mode = 'w'): 181 181 assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\ 182 182 ''' 'w' (write)'''+\ 183 183 ''' 'r' (read)''' +\ 184 ''' 'a' (append)''' 184 ''' 'a' (append)''' 185 185 186 186 #Create filename 187 #self.filename = create_filename(domain.get_datadir(), 188 #domain.get_name(), extension, len(domain))189 190 191 self.filename = create_filename(domain.get_datadir(), 192 domain.get_name(), extension) 187 #self.filename = create_filename(domain.get_datadir(), 188 #domain.get_name(), extension, len(domain)) 189 190 191 self.filename = create_filename(domain.get_datadir(), 192 domain.get_name(), extension) 193 193 194 194 #print 'F', self.filename … … 196 196 self.number_of_volumes = len(domain) 197 197 self.domain = domain 198 199 200 #FIXME: Should we have a general set_precision function?201 198 199 200 #FIXME: Should we have a general set_precision function? 201 202 202 203 203 … … 207 207 """ 208 208 209 209 210 210 def __init__(self, domain, mode = 'w'): 211 211 from Scientific.IO.NetCDF import NetCDFFile 212 from Numeric import Int, Float, Float32 212 from Numeric import Int, Float, Float32 213 213 214 214 self.precision = Float32 #Use single precision 215 215 216 216 Data_format.__init__(self, domain, 'sww', mode) 217 217 … … 219 219 # NetCDF file definition 220 220 fid = NetCDFFile(self.filename, mode) 221 221 222 222 if mode == 'w': 223 223 224 224 #Create new file 225 225 fid.institution = 'Geoscience Australia' 226 226 fid.description = 'Output from pyvolution suitable for plotting' 227 227 228 228 if domain.smooth: 229 229 fid.smoothing = 'Yes' 230 else: 230 else: 231 231 fid.smoothing = 'No' 232 232 233 fid.order = domain.default_order 234 235 236 237 233 fid.order = domain.default_order 234 235 #Reference point 236 #Start time in seconds since the epoch (midnight 1/1/1970) 237 fid.starttime = domain.starttime 238 238 fid.xllcorner = domain.xllcorner 239 239 fid.yllcorner = domain.yllcorner … … 241 241 242 242 # dimension definitions 243 fid.createDimension('number_of_volumes', self.number_of_volumes) 243 fid.createDimension('number_of_volumes', self.number_of_volumes) 244 244 fid.createDimension('number_of_vertices', 3) 245 245 246 246 if domain.smooth is True: 247 247 fid.createDimension('number_of_points', len(domain.vertexlist)) 248 else: 248 else: 249 249 fid.createDimension('number_of_points', 3*self.number_of_volumes) 250 250 251 251 fid.createDimension('number_of_timesteps', None) #extensible 252 252 … … 255 255 fid.createVariable('y', self.precision, ('number_of_points',)) 256 256 fid.createVariable('elevation', self.precision, ('number_of_points',)) 257 257 258 258 #FIXME: Backwards compatibility 259 259 fid.createVariable('z', self.precision, ('number_of_points',)) 260 260 ################################# 261 261 262 262 fid.createVariable('volumes', Int, ('number_of_volumes', 263 263 'number_of_vertices')) … … 265 265 fid.createVariable('time', self.precision, 266 266 ('number_of_timesteps',)) 267 267 268 268 fid.createVariable('stage', self.precision, 269 269 ('number_of_timesteps', … … 279 279 280 280 #Close 281 fid.close() 281 fid.close() 282 282 283 283 … … 297 297 #Get NetCDF 298 298 fid = NetCDFFile(self.filename, 'a') #Open existing file for append 299 299 300 300 # Get the variables 301 301 x = fid.variables['x'] 302 302 y = fid.variables['y'] 303 303 z = fid.variables['elevation'] 304 304 305 305 volumes = fid.variables['volumes'] 306 306 … … 309 309 X,Y,Z,V = Q.get_vertex_values(xy=True, 310 310 precision = self.precision) 311 311 312 312 313 313 … … 315 315 y[:] = Y.astype(self.precision) 316 316 z[:] = Z.astype(self.precision) 317 317 318 318 #FIXME: Backwards compatibility 319 319 z = fid.variables['z'] 320 z[:] = Z.astype(self.precision) 320 z[:] = Z.astype(self.precision) 321 321 ################################ 322 322 323 323 volumes[:] = V.astype(volumes.typecode()) 324 324 325 325 #Close 326 326 fid.close() … … 334 334 from time import sleep 335 335 336 336 337 337 #Get NetCDF 338 338 retries = 0 339 339 file_open = False 340 340 while not file_open and retries < 10: 341 try: 341 try: 342 342 fid = NetCDFFile(self.filename, 'a') #Open existing file 343 343 except IOError: … … 352 352 else: 353 353 file_open = True 354 354 355 355 if not file_open: 356 356 msg = 'File %s could not be opened for append' %self.filename 357 357 raise msg 358 359 358 359 360 360 domain = self.domain 361 361 362 362 # Get the variables 363 363 time = fid.variables['time'] 364 364 stage = fid.variables['stage'] 365 365 xmomentum = fid.variables['xmomentum'] 366 ymomentum = fid.variables['ymomentum'] 366 ymomentum = fid.variables['ymomentum'] 367 367 i = len(time) 368 368 … … 374 374 names = [names] 375 375 376 for name in names: 376 for name in names: 377 377 # Get quantity 378 378 Q = domain.quantities[name] … … 386 386 xmomentum[i,:] = A.astype(self.precision) 387 387 elif name == 'ymomentum': 388 ymomentum[i,:] = A.astype(self.precision) 389 388 ymomentum[i,:] = A.astype(self.precision) 389 390 390 #As in.... 391 391 #eval( name + '[i,:] = A.astype(self.precision)' ) 392 392 #FIXME: But we need a UNIT test for that before refactoring 393 394 393 394 395 395 396 396 #Flush and close … … 404 404 """ 405 405 406 406 407 407 def __init__(self, domain, mode = 'w'): 408 408 from Scientific.IO.NetCDF import NetCDFFile 409 from Numeric import Int, Float, Float 409 from Numeric import Int, Float, Float 410 410 411 411 self.precision = Float #Use full precision 412 412 413 413 Data_format.__init__(self, domain, 'sww', mode) 414 414 … … 416 416 # NetCDF file definition 417 417 fid = NetCDFFile(self.filename, mode) 418 418 419 419 if mode == 'w': 420 420 #Create new file 421 421 fid.institution = 'Geoscience Australia' 422 422 fid.description = 'Checkpoint data' 423 #fid.smooth = domain.smooth 424 fid.order = domain.default_order 423 #fid.smooth = domain.smooth 424 fid.order = domain.default_order 425 425 426 426 # dimension definitions 427 fid.createDimension('number_of_volumes', self.number_of_volumes) 427 fid.createDimension('number_of_volumes', self.number_of_volumes) 428 428 fid.createDimension('number_of_vertices', 3) 429 429 … … 438 438 fid.createVariable('y', self.precision, ('number_of_points',)) 439 439 440 440 441 441 fid.createVariable('volumes', Int, ('number_of_volumes', 442 442 'number_of_vertices')) … … 452 452 453 453 #Close 454 fid.close() 454 fid.close() 455 455 456 456 … … 471 471 #Get NetCDF 472 472 fid = NetCDFFile(self.filename, 'a') #Open existing file for append 473 473 474 474 # Get the variables 475 475 x = fid.variables['x'] 476 476 y = fid.variables['y'] 477 477 478 478 volumes = fid.variables['volumes'] 479 479 … … 482 482 X,Y,Z,V = Q.get_vertex_values(xy=True, 483 483 precision = self.precision) 484 484 485 485 486 486 … … 490 490 491 491 volumes[:] = V 492 492 493 493 #Close 494 494 fid.close() … … 500 500 from Scientific.IO.NetCDF import NetCDFFile 501 501 from time import sleep 502 502 503 503 #Get NetCDF 504 504 retries = 0 505 505 file_open = False 506 506 while not file_open and retries < 10: 507 try: 507 try: 508 508 fid = NetCDFFile(self.filename, 'a') #Open existing file 509 509 except IOError: … … 518 518 else: 519 519 file_open = True 520 520 521 521 if not file_open: 522 522 msg = 'File %s could not be opened for append' %self.filename 523 523 raise msg 524 525 524 525 526 526 domain = self.domain 527 527 528 528 # Get the variables 529 529 time = fid.variables['time'] 530 stage = fid.variables['stage'] 530 stage = fid.variables['stage'] 531 531 i = len(time) 532 532 … … 538 538 A,V = Q.get_vertex_values(xy=False, 539 539 precision = self.precision) 540 540 541 541 stage[i,:] = A.astype(self.precision) 542 542 … … 560 560 def __init__(self, domain, mode = 'w'): 561 561 from Scientific.IO.NetCDF import NetCDFFile 562 from Numeric import Int, Float, Float32 562 from Numeric import Int, Float, Float32 563 563 564 564 self.precision = Float32 #Use single precision 565 565 566 566 Data_format.__init__(self, domain, 'xya', mode) 567 568 569 570 #FIXME -This is the old xya format 567 568 569 570 #FIXME -This is the old xya format 571 571 def store_all(self): 572 572 """Specialisation of store all for xya format … … 579 579 580 580 domain = self.domain 581 581 582 582 fd = open(self.filename, 'w') 583 583 … … 585 585 if domain.smooth is True: 586 586 number_of_points = len(domain.vertexlist) 587 else: 587 else: 588 588 number_of_points = 3*self.number_of_volumes 589 589 … … 613 613 s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict) 614 614 fd.write(s) 615 615 616 616 #close 617 617 fd.close() … … 622 622 """ 623 623 pass 624 624 625 625 626 626 #Auxiliary … … 634 634 635 635 e.g. the x coordinate of vertex 1 of element i is in x[i,1] 636 636 637 637 """ 638 638 #print 'Writing obj to %s' % filename … … 645 645 else: 646 646 FN = filename + '.obj' 647 647 648 648 649 649 outfile = open(FN, 'wb') … … 669 669 """Convert netcdf based data output to obj 670 670 """ 671 from Scientific.IO.NetCDF import NetCDFFile 671 from Scientific.IO.NetCDF import NetCDFFile 672 672 673 673 from Numeric import Float, zeros … … 678 678 fid = NetCDFFile(FN, 'r') #Open existing file for read 679 679 680 680 681 681 # Get the variables 682 682 x = fid.variables['x'] 683 683 y = fid.variables['y'] 684 z = fid.variables['elevation'] 684 z = fid.variables['elevation'] 685 685 time = fid.variables['time'] 686 686 stage = fid.variables['stage'] 687 687 688 688 M = size #Number of lines 689 689 xx = zeros((M,3), Float) … … 696 696 yy[i,j] = y[i+j*M] 697 697 zz[i,j] = z[i+j*M] 698 698 699 699 #Write obj for bathymetry 700 FN = create_filename('.', basefilename, 'obj', size) 700 FN = create_filename('.', basefilename, 'obj', size) 701 701 write_obj(FN,xx,yy,zz) 702 702 … … 716 716 #Write obj for variable data 717 717 #FN = create_filename(basefilename, 'obj', size, time=t) 718 FN = create_filename('.', basefilename[:5], 'obj', size, time=t) 718 FN = create_filename('.', basefilename[:5], 'obj', size, time=t) 719 719 write_obj(FN,xx,yy,zz) 720 720 … … 727 727 import glob, os 728 728 from config import data_dir 729 730 729 730 731 731 #Get bathymetry and x,y's 732 732 lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines() 733 733 734 734 from Numeric import zeros, Float 735 735 … … 742 742 for i, line in enumerate(lines): 743 743 tokens = line.split() 744 values = map(float,tokens) 744 values = map(float,tokens) 745 745 746 746 for j in range(3): 747 747 x[i,j] = values[j*3] 748 748 y[i,j] = values[j*3+1] 749 z[i,j] = values[j*3+2] 750 749 z[i,j] = values[j*3+2] 750 751 751 ##i += 1 752 752 … … 774 774 #Skip bathymetry file 775 775 continue 776 777 i0 += 6 #Position where time starts 778 i1 = filename.find('.dat') 776 777 i0 += 6 #Position where time starts 778 i1 = filename.find('.dat') 779 779 780 780 if i1 > i0: … … 782 782 else: 783 783 raise 'Hmmmm' 784 785 786 784 785 786 787 787 ##i = 0 788 788 for i, line in enumerate(lines): 789 789 tokens = line.split() 790 values = map(float,tokens) 790 values = map(float,tokens) 791 791 792 792 for j in range(3): 793 z[i,j] = values[j] 794 793 z[i,j] = values[j] 794 795 795 ##i += 1 796 796 … … 803 803 nettcdf file filename2 804 804 """ 805 from Scientific.IO.NetCDF import NetCDFFile 806 805 from Scientific.IO.NetCDF import NetCDFFile 806 807 807 #Get NetCDF 808 808 infile = NetCDFFile(filename1, 'r') #Open existing file for read 809 809 outfile = NetCDFFile(filename2, 'w') #Open new file 810 810 811 811 812 812 #Copy dimensions … … 817 817 var = infile.variables[name] 818 818 outfile.createVariable(name, var.typecode(), var.dimensions) 819 820 819 820 821 821 #Copy the static variables 822 822 for name in infile.variables: … … 825 825 else: 826 826 #Copy 827 outfile.variables[name][:] = infile.variables[name][:] 828 829 #Copy selected timesteps 827 outfile.variables[name][:] = infile.variables[name][:] 828 829 #Copy selected timesteps 830 830 time = infile.variables['time'] 831 831 stage = infile.variables['stage'] 832 832 833 833 newtime = outfile.variables['time'] 834 newstage = outfile.variables['stage'] 834 newstage = outfile.variables['stage'] 835 835 836 836 if last is None: … … 841 841 print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j]) 842 842 newtime[i] = time[j] 843 newstage[i,:] = stage[j,:] 844 843 newstage[i,:] = stage[j,:] 844 845 845 #Close 846 846 infile.close() … … 873 873 points: (Nx2) Float array 874 874 elevation: N Float array 875 """ 875 """ 876 876 877 877 import os 878 878 from Scientific.IO.NetCDF import NetCDFFile 879 from Numeric import Float, arrayrange, concatenate 879 from Numeric import Float, arrayrange, concatenate 880 880 881 881 root = basename_in … … 885 885 886 886 if verbose: print 'Reading DEM from %s' %(root + '.dem') 887 887 888 888 ncols = infile.ncols[0] 889 889 nrows = infile.nrows[0] … … 902 902 datum = infile.datum 903 903 units = infile.units 904 904 905 905 906 906 #Get output file … … 909 909 else: 910 910 ptsname = basename_out + '.pts' 911 911 912 912 if verbose: print 'Store to NetCDF file %s' %ptsname 913 913 # NetCDF file definition 914 914 outfile = NetCDFFile(ptsname, 'w') 915 915 916 916 #Create new file 917 917 outfile.institution = 'Geoscience Australia' … … 925 925 outfile.false_easting = false_easting 926 926 outfile.false_northing =false_northing 927 927 928 928 outfile.projection = projection 929 929 outfile.datum = datum … … 935 935 outfile.nrows = nrows 936 936 937 937 938 938 # dimension definitions 939 outfile.createDimension('number_of_points', nrows*ncols) 939 outfile.createDimension('number_of_points', nrows*ncols) 940 940 outfile.createDimension('number_of_dimensions', 2) #This is 2d data 941 941 942 942 # variable definitions 943 943 outfile.createVariable('points', Float, ('number_of_points', … … 953 953 for i in range(nrows): 954 954 if verbose: print 'Processing row %d of %d' %(i, nrows) 955 956 y = (nrows-i)*cellsize 955 956 y = (nrows-i)*cellsize 957 957 for j in range(ncols): 958 958 index = i*ncols + j 959 959 960 960 x = j*cellsize 961 961 points[index, :] = [x,y] 962 elevation[index] = dem_elevation[i, j] 963 964 infile.close() 962 elevation[index] = dem_elevation[i, j] 963 964 infile.close() 965 965 outfile.close() 966 966 967 967 968 968 def convert_dem_from_ascii2netcdf(basename_in, basename_out = None, 969 969 verbose=False): … … 988 988 989 989 The prj format is assumed to be as 990 990 991 991 Projection UTM 992 992 Zone 56 … … 998 998 Yshift 10000000.0000000000 999 999 Parameters 1000 """ 1000 """ 1001 1001 1002 1002 import os … … 1010 1010 # Read Meta data 1011 1011 if verbose: print 'Reading METADATA from %s' %root + '.prj' 1012 metadatafile = open(root + '.prj') 1012 metadatafile = open(root + '.prj') 1013 1013 metalines = metadatafile.readlines() 1014 1014 metadatafile.close() … … 1020 1020 L = metalines[1].strip().split() 1021 1021 assert L[0].strip().lower() == 'zone' 1022 zone = int(L[1].strip()) 1022 zone = int(L[1].strip()) 1023 1023 1024 1024 L = metalines[2].strip().split() 1025 1025 assert L[0].strip().lower() == 'datum' 1026 datum = L[1].strip() #TEXT 1027 1028 L = metalines[3].strip().split() 1026 datum = L[1].strip() #TEXT 1027 1028 L = metalines[3].strip().split() 1029 1029 assert L[0].strip().lower() == 'zunits' #IGNORE 1030 1030 zunits = L[1].strip() #TEXT … … 1039 1039 1040 1040 L = metalines[6].strip().split() 1041 assert L[0].strip().lower() == 'xshift' 1041 assert L[0].strip().lower() == 'xshift' 1042 1042 false_easting = float(L[1].strip()) 1043 1043 1044 1044 L = metalines[7].strip().split() 1045 1045 assert L[0].strip().lower() == 'yshift' 1046 false_northing = float(L[1].strip()) 1046 false_northing = float(L[1].strip()) 1047 1047 1048 1048 #print false_easting, false_northing, zone, datum … … 1051 1051 ########################################### 1052 1052 #Read DEM data 1053 1053 1054 1054 datafile = open(basename_in + '.asc') 1055 1055 … … 1059 1059 1060 1060 if verbose: print 'Got', len(lines), ' lines' 1061 1061 1062 1062 ncols = int(lines[0].split()[1].strip()) 1063 1063 nrows = int(lines[1].split()[1].strip()) … … 1076 1076 netcdfname = root + '.dem' 1077 1077 else: 1078 netcdfname = basename_out + '.dem' 1079 1078 netcdfname = basename_out + '.dem' 1079 1080 1080 if verbose: print 'Store to NetCDF file %s' %netcdfname 1081 1081 # NetCDF file definition 1082 1082 fid = NetCDFFile(netcdfname, 'w') 1083 1083 1084 1084 #Create new file 1085 1085 fid.institution = 'Geoscience Australia' … … 1096 1096 fid.zone = zone 1097 1097 fid.false_easting = false_easting 1098 fid.false_northing = false_northing 1098 fid.false_northing = false_northing 1099 1099 fid.projection = projection 1100 1100 fid.datum = datum 1101 1101 fid.units = units 1102 1102 1103 1103 1104 1104 # dimension definitions 1105 1105 fid.createDimension('number_of_rows', nrows) 1106 fid.createDimension('number_of_columns', ncols) 1106 fid.createDimension('number_of_columns', ncols) 1107 1107 1108 1108 # variable definitions … … 1122 1122 fid.close() 1123 1123 1124 1124 1125 1125 1126 1126 def ferret2sww(basename_in, basename_out = None, … … 1153 1153 ferret2sww uses grid points as vertices in a triangular grid 1154 1154 counting vertices from lower left corner upwards, then right 1155 """ 1155 """ 1156 1156 1157 1157 import os 1158 1158 from Scientific.IO.NetCDF import NetCDFFile 1159 1159 from Numeric import Float, Int, Int32, searchsorted, zeros, array 1160 precision = Float 1160 precision = Float 1161 1161 1162 1162 … … 1165 1165 file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm) 1166 1166 file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s) 1167 file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s) 1167 file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s) 1168 1168 1169 1169 if basename_out is None: 1170 1170 swwname = basename_in + '.sww' 1171 1171 else: 1172 swwname = basename_out + '.sww' 1172 swwname = basename_out + '.sww' 1173 1173 1174 1174 times = file_h.variables['TIME'] … … 1178 1178 if mint == None: 1179 1179 jmin = 0 1180 else: 1180 else: 1181 1181 jmin = searchsorted(times, mint) 1182 1182 1183 1183 if maxt == None: 1184 1184 jmax=len(times) 1185 else: 1185 else: 1186 1186 jmax = searchsorted(times, maxt) 1187 1187 1188 1188 if minlat == None: 1189 1189 kmin=0 … … 1198 1198 if minlon == None: 1199 1199 lmin=0 1200 else: 1200 else: 1201 1201 lmin = searchsorted(longitudes, minlon) 1202 1202 1203 1203 if maxlon == None: 1204 1204 lmax = len(longitudes) 1205 1205 else: 1206 lmax = searchsorted(longitudes, maxlon) 1207 1208 times = times[jmin:jmax] 1206 lmax = searchsorted(longitudes, maxlon) 1207 1208 times = times[jmin:jmax] 1209 1209 latitudes = latitudes[kmin:kmax] 1210 1210 longitudes = longitudes[lmin:lmax] … … 1214 1214 amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax] 1215 1215 xspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] 1216 yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] 1216 yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] 1217 1217 1218 1218 number_of_times = times.shape[0] … … 1222 1222 assert amplitudes.shape[0] == number_of_times 1223 1223 assert amplitudes.shape[1] == number_of_latitudes 1224 assert amplitudes.shape[2] == number_of_longitudes 1224 assert amplitudes.shape[2] == number_of_longitudes 1225 1225 1226 1226 #print times … … 1229 1229 1230 1230 #print 'MIN', min(min(min(amplitudes))) 1231 #print 'MAX', max(max(max(amplitudes))) 1231 #print 'MAX', max(max(max(amplitudes))) 1232 1232 1233 1233 #print number_of_latitudes, number_of_longitudes 1234 1234 number_of_points = number_of_latitudes*number_of_longitudes 1235 1235 number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2 1236 1236 1237 1237 #print file_h.dimensions.keys() 1238 #print file_h.variables.keys() 1238 #print file_h.variables.keys() 1239 1239 1240 1240 file_h.close() 1241 1241 file_u.close() 1242 file_v.close() 1242 file_v.close() 1243 1243 1244 1244 … … 1246 1246 # NetCDF file definition 1247 1247 outfile = NetCDFFile(swwname, 'w') 1248 1248 1249 1249 #Create new file 1250 1250 outfile.institution = 'Geoscience Australia' … … 1256 1256 1257 1257 #For sww compatibility 1258 outfile.smoothing = 'Yes' 1258 outfile.smoothing = 'Yes' 1259 1259 outfile.order = 1 1260 1260 1261 1261 #Start time in seconds since the epoch (midnight 1/1/1970) 1262 1262 outfile.starttime = times[0] 1263 1263 1264 1264 # dimension definitions 1265 1265 outfile.createDimension('number_of_volumes', number_of_volumes) … … 1267 1267 outfile.createDimension('number_of_vertices', 3) 1268 1268 outfile.createDimension('number_of_points', number_of_points) 1269 1270 1269 1270 1271 1271 #outfile.createDimension('number_of_timesteps', len(times)) 1272 1272 outfile.createDimension('number_of_timesteps', len(times)) … … 1276 1276 outfile.createVariable('y', precision, ('number_of_points',)) 1277 1277 outfile.createVariable('elevation', precision, ('number_of_points',)) 1278 1278 1279 1279 #FIXME: Backwards compatibility 1280 1280 outfile.createVariable('z', precision, ('number_of_points',)) 1281 1281 ################################# 1282 1282 1283 1283 outfile.createVariable('volumes', Int, ('number_of_volumes', 1284 1284 'number_of_vertices')) 1285 1285 1286 1286 outfile.createVariable('time', precision, 1287 1287 ('number_of_timesteps',)) 1288 1288 1289 1289 outfile.createVariable('stage', precision, 1290 1290 ('number_of_timesteps', … … 1294 1294 ('number_of_timesteps', 1295 1295 'number_of_points')) 1296 1296 1297 1297 outfile.createVariable('ymomentum', precision, 1298 1298 ('number_of_timesteps', … … 1308 1308 1309 1309 #Check zone boundaries 1310 refzone, _, _ = redfearn(latitudes[0],longitudes[0]) 1310 refzone, _, _ = redfearn(latitudes[0],longitudes[0]) 1311 1311 1312 1312 vertices = {} 1313 for k, lat in enumerate(latitudes): 1314 for l, lon in enumerate(longitudes): 1313 for k, lat in enumerate(latitudes): 1314 for l, lon in enumerate(longitudes): 1315 1315 1316 1316 vertices[l,k] = i 1317 1317 1318 zone, easting, northing = redfearn(lat,lon) 1318 zone, easting, northing = redfearn(lat,lon) 1319 1319 1320 1320 msg = 'Zone boundary crossed at longitude =', lon … … 1331 1331 for k in range(number_of_latitudes-1): 1332 1332 v1 = vertices[l,k+1] 1333 v2 = vertices[l,k] 1334 v3 = vertices[l+1,k+1] 1335 v4 = vertices[l+1,k] 1333 v2 = vertices[l,k] 1334 v3 = vertices[l+1,k+1] 1335 v4 = vertices[l+1,k] 1336 1336 1337 1337 volumes.append([v1,v2,v3]) #Upper element 1338 volumes.append([v4,v3,v2]) #Lower element 1339 1340 volumes = array(volumes) 1338 volumes.append([v4,v3,v2]) #Lower element 1339 1340 volumes = array(volumes) 1341 1341 1342 1342 if origin == None: 1343 zone = refzone 1343 zone = refzone 1344 1344 xllcorner = min(x) 1345 1345 yllcorner = min(y) … … 1347 1347 zone = origin[0] 1348 1348 xllcorner = origin[1] 1349 yllcorner = origin[2] 1350 1351 1349 yllcorner = origin[2] 1350 1351 1352 1352 outfile.xllcorner = xllcorner 1353 1353 outfile.yllcorner = yllcorner … … 1358 1358 outfile.variables['z'][:] = 0.0 1359 1359 outfile.variables['elevation'][:] = 0.0 1360 outfile.variables['time'][:] = times 1360 outfile.variables['time'][:] = times 1361 1361 outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 1362 1362 1363 1363 1364 1364 … … 1366 1366 stage = outfile.variables['stage'] 1367 1367 xmomentum = outfile.variables['xmomentum'] 1368 ymomentum = outfile.variables['ymomentum'] 1368 ymomentum = outfile.variables['ymomentum'] 1369 1369 1370 1370 for j in range(len(times)): 1371 1371 i = 0 1372 1372 for k in range(number_of_latitudes): 1373 for l in range(number_of_longitudes): 1373 for l in range(number_of_longitudes): 1374 1374 h = zscale*amplitudes[j,k,l]/100 + mean_stage 1375 1375 stage[j,i] = h 1376 1376 xmomentum[j,i] = xspeed[j,k,l]/100*h 1377 ymomentum[j,i] = yspeed[j,k,l]/100*h 1377 ymomentum[j,i] = yspeed[j,k,l]/100*h 1378 1378 i += 1 1379 1380 outfile.close() 1381 1382 1379 1380 outfile.close() 1381 1382 1383 1383 def extent_sww(file_name): 1384 1384 """ 1385 1385 Read in an sww file. 1386 1386 1387 1387 Input; 1388 1388 file_name - the sww file 1389 1389 1390 1390 Output; 1391 1391 z - Vector of bed elevation … … 1395 1395 stage - array with respect to time and vertices (x,y) 1396 1396 """ 1397 1398 1399 from Scientific.IO.NetCDF import NetCDFFile 1400 1397 1398 1399 from Scientific.IO.NetCDF import NetCDFFile 1400 1401 1401 #Check contents 1402 1402 #Get NetCDF 1403 fid = NetCDFFile(file_name, 'r') 1404 1403 fid = NetCDFFile(file_name, 'r') 1404 1405 1405 # Get the variables 1406 1406 x = fid.variables['x'][:] … … 1413 1413 1414 1414 fid.close() 1415 1415 1416 1416 return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)] 1417 1417 … … 1426 1426 1427 1427 M, N = v0.shape 1428 1428 1429 1429 FN = create_filename(filename, 'cpt', M, t) 1430 1430 #print 'Writing to %s' %FN 1431 1431 1432 1432 fid = open(FN, 'w') 1433 1433 for i in range(M): 1434 for j in range(N):1435 fid.write('%.16e ' %v0[i,j])1436 1434 for j in range(N): 1437 fid.write('%.16e ' %v 1[i,j])1435 fid.write('%.16e ' %v0[i,j]) 1438 1436 for j in range(N): 1439 fid.write('%.16e ' %v2[i,j]) 1440 1437 fid.write('%.16e ' %v1[i,j]) 1438 for j in range(N): 1439 fid.write('%.16e ' %v2[i,j]) 1440 1441 1441 fid.write('\n') 1442 1442 fid.close() … … 1448 1448 1449 1449 M, N = v0.shape 1450 1450 1451 1451 FN = create_filename(filename, 'cpt', M, t) 1452 1452 #print 'Reading from %s' %FN 1453 1453 1454 1454 fid = open(FN) 1455 1455 … … 1462 1462 v1[i,j] = float(values[3+j]) 1463 1463 v2[i,j] = float(values[6+j]) 1464 1464 1465 1465 fid.close() 1466 1466 … … 1474 1474 1475 1475 print X0 1476 import sys; sys.exit() 1476 import sys; sys.exit() 1477 1477 FN = create_filename(filename, 'cpt', M) 1478 1478 print 'Writing to %s' %FN 1479 1479 1480 1480 fid = open(FN, 'w') 1481 1481 for i in range(M): 1482 for j in range(2): 1482 for j in range(2): 1483 1483 fid.write('%.16e ' %X0[i,j]) #x, y 1484 for j in range(N): 1484 for j in range(N): 1485 1485 fid.write('%.16e ' %v0[i,j]) #z,z,z, 1486 1487 for j in range(2): 1486 1487 for j in range(2): 1488 1488 fid.write('%.16e ' %X1[i,j]) #x, y 1489 for j in range(N): 1490 fid.write('%.16e ' %v1[i,j]) 1491 1492 for j in range(2): 1489 for j in range(N): 1490 fid.write('%.16e ' %v1[i,j]) 1491 1492 for j in range(2): 1493 1493 fid.write('%.16e ' %X2[i,j]) #x, y 1494 for j in range(N): 1495 fid.write('%.16e ' %v2[i,j]) 1496 1494 for j in range(N): 1495 fid.write('%.16e ' %v2[i,j]) 1496 1497 1497 fid.write('\n') 1498 1498 fid.close() … … 1508 1508 1509 1509 M, N = v0.shape 1510 1510 1511 1511 FN = create_filename(filename, 'dat', M) 1512 1512 #print 'Writing to %s' %FN 1513 1513 1514 1514 fid = open(FN, 'w') 1515 1515 for i in range(M): 1516 for j in range(2): 1516 for j in range(2): 1517 1517 fid.write('%f ' %X0[i,j]) #x, y 1518 1518 fid.write('%f ' %v0[i,0]) #z 1519 1520 for j in range(2): 1519 1520 for j in range(2): 1521 1521 fid.write('%f ' %X1[i,j]) #x, y 1522 1522 fid.write('%f ' %v1[i,0]) #z 1523 1524 for j in range(2): 1523 1524 for j in range(2): 1525 1525 fid.write('%f ' %X2[i,j]) #x, y 1526 1526 fid.write('%f ' %v2[i,0]) #z 1527 1527 1528 1528 fid.write('\n') 1529 1529 fid.close() … … 1536 1536 1537 1537 M, N = v0.shape 1538 1538 1539 1539 FN = create_filename(filename, 'dat', M, t) 1540 1540 #print 'Writing to %s' %FN 1541 1541 1542 1542 fid = open(FN, 'w') 1543 1543 for i in range(M): 1544 1544 fid.write('%.4f ' %v0[i,0]) 1545 1545 fid.write('%.4f ' %v1[i,0]) 1546 fid.write('%.4f ' %v2[i,0]) 1546 fid.write('%.4f ' %v2[i,0]) 1547 1547 1548 1548 fid.write('\n') 1549 1549 fid.close() 1550 -
inundation/ga/storm_surge/pyvolution/test_advection.py
r773 r1018 8 8 from advection import * 9 9 10 11 class Test Case(unittest.TestCase):10 11 class Test_Advection(unittest.TestCase): 12 12 def setUp(self): 13 13 pass 14 14 15 15 def tearDown(self): 16 16 pass … … 25 25 26 26 points = [a, b, c, d, e, f] 27 #bac, bce, ecf, dbe, daf, dae 27 #bac, bce, ecf, dbe, daf, dae 28 28 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]] 29 30 domain = Domain(points, vertices) 29 30 domain = Domain(points, vertices) 31 31 domain.check_integrity() 32 32 33 33 assert domain.quantities.has_key('stage') 34 34 … … 38 38 def test_flux_1_triangle0(self): 39 39 a = [0.0, 0.5] 40 b = [0.0, 0.0] 40 b = [0.0, 0.0] 41 41 c = [0.5, 0.5] 42 42 43 43 points = [a, b, c] 44 44 vertices = [ [0,1,2] ] 45 domain = Domain(points, vertices) 45 domain = Domain(points, vertices) 46 46 domain.check_integrity() 47 47 48 48 49 49 #Populate boundary array with dirichlet conditions. 50 domain.neighbours = array([[-1,-2,-3]]) 50 domain.neighbours = array([[-1,-2,-3]]) 51 51 domain.quantities['stage'].boundary_values[:] = 1.0 52 52 … … 54 54 55 55 domain.distribute_to_vertices_and_edges() #Use first order default 56 56 57 57 domain.check_integrity() 58 58 … … 61 61 R = -0.5/domain.areas[0] 62 62 63 assert U==R, '%s %s' %(U, R) 63 assert U==R, '%s %s' %(U, R) 64 64 65 65 66 66 def test_flux_1_triangle1(self): 67 67 68 68 a = [0.0, 0.5] 69 b = [0.0, 0.0] 69 b = [0.0, 0.0] 70 70 c = [0.5, 0.5] 71 71 72 72 points = [a, b, c] 73 73 vertices = [ [0,1,2] ] 74 domain = Domain(points, vertices) 74 domain = Domain(points, vertices) 75 75 domain.check_integrity() 76 76 … … 79 79 domain.distribute_to_vertices_and_edges() 80 80 domain.check_integrity() 81 82 81 82 83 83 domain.compute_fluxes() 84 84 U = -domain.quantities['stage'].explicit_update 85 85 R = 0.5/domain.areas[0] 86 86 87 assert U==R, '%s %s' %(U, R) 87 assert U==R, '%s %s' %(U, R) 88 88 89 89 90 90 91 91 def test_flux_1_triangle2(self): 92 92 93 93 a = [0.0, 0.5] 94 b = [0.0, 0.0] 94 b = [0.0, 0.0] 95 95 c = [0.5, 0.5] 96 96 97 97 points = [a, b, c] 98 98 vertices = [ [0,1,2] ] 99 domain = Domain(points, vertices) 99 domain = Domain(points, vertices) 100 100 domain.check_integrity() 101 101 102 102 103 103 #Populate boundary array with dirichlet conditions. 104 domain.neighbours = array([[-1,-2,-3]]) 104 domain.neighbours = array([[-1,-2,-3]]) 105 105 domain.quantities['stage'].boundary_values[0] = 1.0 106 106 107 107 domain.distribute_to_vertices_and_edges() #Use first order default 108 108 109 109 domain.check_integrity() 110 110 … … 113 113 assert allclose(U, 0) 114 114 115 115 116 116 117 117 … … 122 122 123 123 a = [0.0, 0.5] 124 b = [0.0, 0.0] 124 b = [0.0, 0.0] 125 125 c = [0.5, 0.5] 126 d = [0.5, 0.0] 126 d = [0.5, 0.0] 127 127 128 128 points = [a, b, c, d] 129 129 vertices = [ [0,1,2], [3,2,1] ] 130 domain = Domain(points, vertices) 130 domain = Domain(points, vertices) 131 131 domain.check_integrity() 132 132 … … 137 137 domain.distribute_to_vertices_and_edges() 138 138 139 domain.compute_fluxes() 139 domain.compute_fluxes() 140 140 141 141 X = domain.quantities['stage'].explicit_update 142 142 assert X[0] == -X[1] 143 143 144 144 145 145 def test_advection_example(self): 146 146 #Test that system can evolve 147 147 148 148 from mesh_factory import rectangular 149 149 150 150 points, vertices, boundary = rectangular(6, 6) 151 151 152 #Create advection domain with direction (1,-1) 152 #Create advection domain with direction (1,-1) 153 153 domain = Domain(points, vertices, boundary, velocity=[1.0, -1.0]) 154 154 155 155 # Initial condition is zero by default 156 156 157 157 #Boundaries 158 158 T = Transmissive_boundary(domain) 159 159 D = Dirichlet_boundary(array([3.1415])) 160 160 161 161 domain.set_boundary( {'left': D, 'right': T, 'bottom': T, 'top': T} ) 162 162 domain.check_integrity() … … 166 166 if allclose(domain.quantities['stage'].centroid_values, 3.1415): 167 167 break 168 168 169 169 assert allclose(domain.quantities['stage'].centroid_values, 3.1415) 170 170 171 171 172 172 #------------------------------------------------------------- 173 173 if __name__ == "__main__": 174 suite = unittest.makeSuite(Test Case,'test')174 suite = unittest.makeSuite(Test_Advection,'test') 175 175 runner = unittest.TextTestRunner() 176 176 runner.run(suite) -
inundation/ga/storm_surge/pyvolution/test_all.py
r892 r1018 23 23 24 24 import sys 25 25 26 26 files = os.listdir(path) 27 27 … … 38 38 else: 39 39 pass 40 return test_files 40 return test_files 41 41 42 42 43 43 44 44 def regressionTest(): 45 import sys, os, re, unittest 45 import sys, os, re, unittest 46 46 path = os.path.split(sys.argv[0])[0] or os.getcwd() 47 47 … … 53 53 #test = re.compile('^test_[\w]*.py$', re.IGNORECASE) 54 54 #files = filter(test.search, files) 55 55 56 56 57 57 try: 58 58 files.remove(__file__) #Remove self from list (Ver 2.3. or later) 59 59 except: 60 files.remove('test_all.py') 60 files.remove('test_all.py') 61 61 62 62 63 if globals().has_key('exclude'): 63 if globals().has_key('exclude'): 64 64 for file in exclude: 65 65 files.remove(file) 66 66 print 'WARNING: File '+ file + ' excluded from testing' 67 67 68 68 69 69 filenameToModuleName = lambda f: os.path.splitext(f)[0] … … 76 76 77 77 os.system('python compile.py') #Attempt to compile all extensions 78 78 79 #print regressionTest() 79 80 unittest.main(defaultTest='regressionTest') 80 -
inundation/ga/storm_surge/pyvolution/test_cg_solve.py
r599 r1018 9 9 10 10 11 class Test Case(unittest.TestCase):11 class Test_CG_Solve(unittest.TestCase): 12 12 13 13 def test_sparse_solve(self): 14 14 """Small Sparse Matrix""" 15 15 16 16 A = [[2.0, -1.0, 0.0, 0.0 ], 17 17 [-1.0, 2.0, -1.0, 0.0], 18 18 [0.0, -1.0, 2.0, -1.0], 19 19 [0.0,0.0, -1.0, 2.0]] 20 20 21 21 A = Sparse(A) 22 22 … … 35 35 n = 50 36 36 A = Sparse(n,n) 37 37 38 38 for i in arange(0,n): 39 39 A[i,i] = 1.0 … … 42 42 if i < n-1 : 43 43 A[i,i+1] = -0.5 44 44 45 45 xe = ones( (n,), Float) 46 46 … … 52 52 def test_solve_large_2d(self): 53 53 """Standard 2d laplacian""" 54 54 55 55 n = 20 56 56 m = 10 … … 70 70 if j < m-1 : 71 71 A[I,I+1] = -1.0 72 72 73 73 xe = ones( (n*m,), Float) 74 74 … … 81 81 """Standard 2d laplacian with csr format 82 82 """ 83 83 84 84 n = 100 85 85 m = 100 … … 99 99 if j < m-1 : 100 100 A[I,I+1] = -1.0 101 101 102 102 xe = ones( (n*m,), Float) 103 103 … … 114 114 def test_solve_large_2d_with_default_guess(self): 115 115 """Standard 2d laplacian using default first guess""" 116 116 117 117 n = 20 118 118 m = 10 … … 132 132 if j < m-1 : 133 133 A[I,I+1] = -1.0 134 134 135 135 xe = ones( (n*m,), Float) 136 136 … … 143 143 def test_vector_shape_error(self): 144 144 """Raise VectorShapeError""" 145 145 146 146 A = [[2.0, -1.0, 0.0, 0.0 ], 147 147 [-1.0, 2.0, -1.0, 0.0], 148 148 [0.0, -1.0, 2.0, -1.0], 149 149 [0.0,0.0, -1.0, 2.0]] 150 150 151 151 A = Sparse(A) 152 152 … … 161 161 raise msg 162 162 163 163 164 164 #------------------------------------------------------------- 165 165 if __name__ == "__main__": 166 suite = unittest.makeSuite(Test Case,'test')166 suite = unittest.makeSuite(Test_CG_Solve,'test') 167 167 runner = unittest.TextTestRunner() #(verbosity=2) 168 168 runner.run(suite) 169 169 170 171 172 170 173 171 -
inundation/ga/storm_surge/pyvolution/test_data_manager.py
r995 r1018 11 11 from config import epsilon 12 12 13 class dataTestCase(unittest.TestCase):13 class Test_Data_Manager(unittest.TestCase): 14 14 def setUp(self): 15 15 import time … … 24 24 domain.default_order=2 25 25 26 26 27 27 #Set some field values 28 domain.set_quantity('elevation', lambda x,y: -x) 28 domain.set_quantity('elevation', lambda x,y: -x) 29 29 domain.set_quantity('friction', 0.03) 30 30 31 31 32 ###################### 32 ###################### 33 33 # Boundary conditions 34 34 B = Transmissive_boundary(domain) 35 35 domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B}) 36 37 38 ###################### 36 37 38 ###################### 39 39 #Initial condition - with jumps 40 40 … … 45 45 h = 0.3 46 46 for i in range(stage.shape[0]): 47 if i % 2 == 0: 47 if i % 2 == 0: 48 48 stage[i,:] = bed[i,:] + h 49 49 else: 50 50 stage[i,:] = bed[i,:] 51 51 52 52 domain.set_quantity('stage', stage) 53 53 self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values) 54 54 55 domain.distribute_to_vertices_and_edges() 55 domain.distribute_to_vertices_and_edges() 56 56 57 57 … … 61 61 self.X = C[:,0:6:2].copy() 62 62 self.Y = C[:,1:6:2].copy() 63 63 64 64 self.F = bed 65 65 66 66 67 67 def tearDown(self): 68 68 pass 69 69 70 70 71 71 72 72 73 73 # def test_xya(self): … … 79 79 80 80 # domain = self.domain 81 81 82 82 # domain.filename = 'datatest' + str(time.time()) 83 83 # domain.format = 'xya' 84 84 # domain.smooth = True 85 85 86 86 # xya = get_dataobject(self.domain) 87 87 # xya.store_all() … … 99 99 # if domain.smooth: 100 100 # self.failUnless(lFile[0] == '9 3 # <vertex #> <x> <y> [attributes]') 101 # else: 102 # self.failUnless(lFile[0] == '24 3 # <vertex #> <x> <y> [attributes]') 101 # else: 102 # self.failUnless(lFile[0] == '24 3 # <vertex #> <x> <y> [attributes]') 103 103 104 104 # #Get smoothed field values with X and Y … … 108 108 109 109 # Q,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities', 110 # indices = (0,), precision = Float) 111 112 113 110 # indices = (0,), precision = Float) 111 112 113 114 114 # for i, line in enumerate(lFile[1:]): 115 115 # fields = line.split() … … 121 121 # assert allclose(float(fields[2]), F[i,0]) 122 122 # assert allclose(float(fields[3]), Q[i,0]) 123 # assert allclose(float(fields[4]), F[i,1]) 123 # assert allclose(float(fields[4]), F[i,1]) 124 124 125 125 … … 133 133 import time, os 134 134 from Numeric import array, zeros, allclose, Float, concatenate 135 from Scientific.IO.NetCDF import NetCDFFile 135 from Scientific.IO.NetCDF import NetCDFFile 136 136 137 137 self.domain.filename = 'datatest' + str(id(self)) 138 138 self.domain.format = 'sww' 139 139 self.domain.smooth = False 140 140 141 141 sww = get_dataobject(self.domain) 142 142 sww.store_connectivity() … … 145 145 #Get NetCDF 146 146 fid = NetCDFFile(sww.filename, 'r') #Open existing file for append 147 147 148 148 # Get the variables 149 149 x = fid.variables['x'] … … 151 151 z = fid.variables['elevation'] 152 152 153 volumes = fid.variables['volumes'] 153 volumes = fid.variables['volumes'] 154 154 155 155 … … 164 164 assert V[k, 0] == 3*k 165 165 assert V[k, 1] == 3*k+1 166 assert V[k, 2] == 3*k+2 167 168 169 fid.close() 170 166 assert V[k, 2] == 3*k+2 167 168 169 fid.close() 170 171 171 #Cleanup 172 172 os.remove(sww.filename) … … 180 180 import time, os 181 181 from Numeric import array, zeros, allclose, Float, concatenate 182 from Scientific.IO.NetCDF import NetCDFFile 183 184 self.domain.filename = 'datatest' + str(id(self)) 182 from Scientific.IO.NetCDF import NetCDFFile 183 184 self.domain.filename = 'datatest' + str(id(self)) 185 185 self.domain.format = 'sww' 186 186 self.domain.smooth = True 187 187 188 188 sww = get_dataobject(self.domain) 189 sww.store_connectivity() 189 sww.store_connectivity() 190 190 191 191 #Check contents 192 192 #Get NetCDF 193 193 fid = NetCDFFile(sww.filename, 'r') #Open existing file for append 194 194 195 195 # Get the variables 196 196 x = fid.variables['x'] … … 198 198 z = fid.variables['elevation'] 199 199 200 volumes = fid.variables['volumes'] 200 volumes = fid.variables['volumes'] 201 201 202 202 X = x[:] 203 203 Y = y[:] 204 204 205 205 assert allclose([X[0], Y[0]], array([0.0, 0.0])) 206 206 assert allclose([X[1], Y[1]], array([0.0, 0.5])) … … 209 209 assert allclose([X[4], Y[4]], array([0.5, 0.5])) 210 210 211 assert allclose([X[7], Y[7]], array([1.0, 0.5])) 211 assert allclose([X[7], Y[7]], array([1.0, 0.5])) 212 212 213 213 Z = z[:] … … 221 221 assert V[4,0] == 6 222 222 assert V[4,1] == 7 223 assert V[4,2] == 3 224 225 226 fid.close() 227 223 assert V[4,2] == 3 224 225 226 fid.close() 227 228 228 #Cleanup 229 229 os.remove(sww.filename) … … 237 237 import time, os 238 238 from Numeric import array, zeros, allclose, Float, concatenate 239 from Scientific.IO.NetCDF import NetCDFFile 239 from Scientific.IO.NetCDF import NetCDFFile 240 240 241 241 self.domain.filename = 'datatest' + str(id(self)) 242 242 self.domain.format = 'sww' 243 243 self.domain.smooth = True 244 self.domain.reduction = mean 245 244 self.domain.reduction = mean 245 246 246 sww = get_dataobject(self.domain) 247 247 sww.store_connectivity() … … 251 251 #Get NetCDF 252 252 fid = NetCDFFile(sww.filename, 'r') #Open existing file for append 253 254 255 # Get the variables 256 x = fid.variables['x'] 257 y = fid.variables['y'] 258 z = fid.variables['elevation'] 259 time = fid.variables['time'] 260 stage = fid.variables['stage'] 261 262 263 Q = self.domain.quantities['stage'] 264 Q0 = Q.vertex_values[:,0] 265 Q1 = Q.vertex_values[:,1] 266 Q2 = Q.vertex_values[:,2] 267 268 A = stage[0,:] 269 #print A[0], (Q2[0,0] + Q1[1,0])/2 270 assert allclose(A[0], (Q2[0] + Q1[1])/2) 271 assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3) 272 assert allclose(A[2], Q0[3]) 273 assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3) 274 275 #Center point 276 assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\ 277 Q0[5] + Q2[6] + Q1[7])/6) 278 279 280 281 fid.close() 282 283 #Cleanup 284 os.remove(sww.filename) 285 286 287 def test_sww_variable2(self): 288 """Test that sww information can be written correctly 289 multiple timesteps. Use average as reduction operator 290 """ 291 292 import time, os 293 from Numeric import array, zeros, allclose, Float, concatenate 294 from Scientific.IO.NetCDF import NetCDFFile 295 296 self.domain.filename = 'datatest' + str(id(self)) 297 self.domain.format = 'sww' 298 self.domain.smooth = True 299 300 self.domain.reduction = mean 301 302 sww = get_dataobject(self.domain) 303 sww.store_connectivity() 304 sww.store_timestep('stage') 305 self.domain.evolve_to_end(finaltime = 0.01) 306 sww.store_timestep('stage') 307 308 309 #Check contents 310 #Get NetCDF 311 fid = NetCDFFile(sww.filename, 'r') #Open existing file for append 253 312 254 313 255 # Get the variables … … 316 258 z = fid.variables['elevation'] 317 259 time = fid.variables['time'] 318 stage = fid.variables['stage'] 319 320 #Check values 321 Q = self.domain.quantities['stage'] 260 stage = fid.variables['stage'] 261 262 263 Q = self.domain.quantities['stage'] 322 264 Q0 = Q.vertex_values[:,0] 323 265 Q1 = Q.vertex_values[:,1] 324 266 Q2 = Q.vertex_values[:,2] 325 267 326 A = stage[1,:] 327 assert allclose(A[0], (Q2[0] + Q1[1])/2) 268 A = stage[0,:] 269 #print A[0], (Q2[0,0] + Q1[1,0])/2 270 assert allclose(A[0], (Q2[0] + Q1[1])/2) 328 271 assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3) 329 272 assert allclose(A[2], Q0[3]) … … 333 276 assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\ 334 277 Q0[5] + Q2[6] + Q1[7])/6) 335 336 337 338 339 340 fid.close() 341 278 279 280 281 fid.close() 282 342 283 #Cleanup 343 284 os.remove(sww.filename) 344 345 def test_sww_variable3(self): 285 286 287 def test_sww_variable2(self): 346 288 """Test that sww information can be written correctly 347 multiple timesteps using a different reduction operator (min)289 multiple timesteps. Use average as reduction operator 348 290 """ 349 291 … … 352 294 from Scientific.IO.NetCDF import NetCDFFile 353 295 354 self.domain.filename = 'datatest' + str(id(self)) 296 self.domain.filename = 'datatest' + str(id(self)) 355 297 self.domain.format = 'sww' 356 298 self.domain.smooth = True 357 self.domain.reduction = min 358 299 300 self.domain.reduction = mean 301 359 302 sww = get_dataobject(self.domain) 360 303 sww.store_connectivity() 361 304 sww.store_timestep('stage') 362 363 305 self.domain.evolve_to_end(finaltime = 0.01) 364 306 sww.store_timestep('stage') 365 307 366 308 367 368 309 #Check contents 369 310 #Get NetCDF 370 fid = NetCDFFile(sww.filename, 'r') 371 311 fid = NetCDFFile(sww.filename, 'r') #Open existing file for append 372 312 373 313 # Get the variables … … 376 316 z = fid.variables['elevation'] 377 317 time = fid.variables['time'] 378 stage = fid.variables['stage'] 318 stage = fid.variables['stage'] 379 319 380 320 #Check values 381 #Check values 382 Q = self.domain.quantities['stage'] 321 Q = self.domain.quantities['stage'] 383 322 Q0 = Q.vertex_values[:,0] 384 323 Q1 = Q.vertex_values[:,1] … … 386 325 387 326 A = stage[1,:] 388 assert allclose(A[0], min(Q2[0], Q1[1])) 327 assert allclose(A[0], (Q2[0] + Q1[1])/2) 328 assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3) 329 assert allclose(A[2], Q0[3]) 330 assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3) 331 332 #Center point 333 assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\ 334 Q0[5] + Q2[6] + Q1[7])/6) 335 336 337 338 339 340 fid.close() 341 342 #Cleanup 343 os.remove(sww.filename) 344 345 def test_sww_variable3(self): 346 """Test that sww information can be written correctly 347 multiple timesteps using a different reduction operator (min) 348 """ 349 350 import time, os 351 from Numeric import array, zeros, allclose, Float, concatenate 352 from Scientific.IO.NetCDF import NetCDFFile 353 354 self.domain.filename = 'datatest' + str(id(self)) 355 self.domain.format = 'sww' 356 self.domain.smooth = True 357 self.domain.reduction = min 358 359 sww = get_dataobject(self.domain) 360 sww.store_connectivity() 361 sww.store_timestep('stage') 362 363 self.domain.evolve_to_end(finaltime = 0.01) 364 sww.store_timestep('stage') 365 366 367 368 #Check contents 369 #Get NetCDF 370 fid = NetCDFFile(sww.filename, 'r') 371 372 373 # Get the variables 374 x = fid.variables['x'] 375 y = fid.variables['y'] 376 z = fid.variables['elevation'] 377 time = fid.variables['time'] 378 stage = fid.variables['stage'] 379 380 #Check values 381 #Check values 382 Q = self.domain.quantities['stage'] 383 Q0 = Q.vertex_values[:,0] 384 Q1 = Q.vertex_values[:,1] 385 Q2 = Q.vertex_values[:,2] 386 387 A = stage[1,:] 388 assert allclose(A[0], min(Q2[0], Q1[1])) 389 389 assert allclose(A[1], min(Q0[1], Q1[3], Q2[2])) 390 390 assert allclose(A[2], Q0[3]) … … 394 394 assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\ 395 395 Q0[5], Q2[6], Q1[7])) 396 397 398 fid.close() 399 396 397 398 fid.close() 399 400 400 #Cleanup 401 401 os.remove(sww.filename) … … 408 408 import time, os, config 409 409 from Numeric import array, zeros, allclose, Float, concatenate 410 from Scientific.IO.NetCDF import NetCDFFile 410 from Scientific.IO.NetCDF import NetCDFFile 411 411 412 412 self.domain.filename = 'synctest' … … 426 426 if t == 0.0: 427 427 assert allclose(stage, self.initial_stage) 428 assert allclose(stage_file[:], stage.flat) 428 assert allclose(stage_file[:], stage.flat) 429 429 else: 430 430 assert not allclose(stage, self.initial_stage) 431 assert not allclose(stage_file[:], stage.flat) 431 assert not allclose(stage_file[:], stage.flat) 432 432 433 433 fid.close() … … 442 442 import time, os 443 443 from Numeric import array, zeros, allclose, Float, concatenate 444 from Scientific.IO.NetCDF import NetCDFFile 445 446 self.domain.filename = 'datatest' + str(id(self)) 444 from Scientific.IO.NetCDF import NetCDFFile 445 446 self.domain.filename = 'datatest' + str(id(self)) 447 447 self.domain.format = 'sww' 448 448 self.domain.smooth = True 449 self.domain.reduction = mean 450 449 self.domain.reduction = mean 450 451 451 sww = get_dataobject(self.domain) 452 452 sww.store_connectivity() … … 455 455 #Check contents 456 456 #Get NetCDF 457 fid = NetCDFFile(sww.filename, 'r') 458 457 fid = NetCDFFile(sww.filename, 'r') 458 459 459 # Get the variables 460 460 x = fid.variables['x'] … … 462 462 z = fid.variables['elevation'] 463 463 464 volumes = fid.variables['volumes'] 464 volumes = fid.variables['volumes'] 465 465 time = fid.variables['time'] 466 466 467 # 2D 468 stage = fid.variables['stage'] 467 # 2D 468 stage = fid.variables['stage'] 469 469 470 470 X = x[:] … … 474 474 T = time[:] 475 475 S = stage[:,:] 476 476 477 477 # print "****************************" 478 478 # print "X ",X … … 489 489 # print "****************************" 490 490 491 492 fid.close() 493 491 492 fid.close() 493 494 494 #Cleanup 495 495 os.remove(sww.filename) … … 522 522 ref_elevation = [] 523 523 for i in range(6): 524 y = (6-i)*25.0 524 y = (6-i)*25.0 525 525 for j in range(5): 526 526 x = j*25.0 … … 530 530 ref_elevation.append(z) 531 531 fid.write('%f ' %z) 532 fid.write('\n') 532 fid.write('\n') 533 533 534 534 fid.close() … … 550 550 """) 551 551 fid.close() 552 552 553 553 #Convert to NetCDF xya 554 554 convert_dem_from_ascii2netcdf(root) … … 558 558 #Get NetCDF 559 559 fid = NetCDFFile(root+'.pts', 'r') 560 560 561 561 # Get the variables 562 562 #print fid.variables.keys() … … 572 572 #print attributes[:] 573 573 #print ref_elevation 574 assert allclose(elevation, ref_elevation) 575 576 #Cleanup 577 fid.close() 578 574 assert allclose(elevation, ref_elevation) 575 576 #Cleanup 577 fid.close() 578 579 579 580 580 os.remove(root + '.pts') 581 os.remove(root + '.dem') 582 os.remove(root + '.asc') 583 os.remove(root + '.prj') 581 os.remove(root + '.dem') 582 os.remove(root + '.asc') 583 os.remove(root + '.prj') 584 584 585 585 … … 590 590 from Scientific.IO.NetCDF import NetCDFFile 591 591 592 #The test file has 592 #The test file has 593 593 # LON = 150.66667, 150.83334, 151, 151.16667 594 594 # LAT = -34.5, -34.33333, -34.16667, -34 ; … … 597 597 # First value (index=0) in small_ha.nc is 0.3400644 cm, 598 598 # Fourth value (index==3) is -6.50198 cm 599 599 600 600 601 601 from coordinate_transforms.redfearn import redfearn … … 605 605 first_value = fid.variables['HA'][:][0,0,0] 606 606 fourth_value = fid.variables['HA'][:][0,0,3] 607 607 608 608 #Call conversion (with zero origin) 609 609 ferret2sww('small', verbose=False, 610 610 origin = (56, 0, 0)) 611 611 612 612 613 613 #Work out the UTM coordinates for first point 614 614 zone, e, n = redfearn(-34.5, 150.66667) 615 615 #print zone, e, n 616 616 617 617 #Read output file 'small.sww' 618 618 fid = NetCDFFile('small.sww') 619 619 620 620 x = fid.variables['x'][:] 621 y = fid.variables['y'][:] 621 y = fid.variables['y'][:] 622 622 623 623 #Check that first coordinate is correctly represented … … 631 631 632 632 #Check fourth value 633 assert allclose(stage[0,3], fourth_value/100) #Meters 633 assert allclose(stage[0,3], fourth_value/100) #Meters 634 634 635 635 fid.close() … … 639 639 os.remove('small.sww') 640 640 641 641 642 642 643 643 def test_ferret2sww_2(self): … … 647 647 from Scientific.IO.NetCDF import NetCDFFile 648 648 649 #The test file has 649 #The test file has 650 650 # LON = 150.66667, 150.83334, 151, 151.16667 651 651 # LAT = -34.5, -34.33333, -34.16667, -34 ; … … 654 654 # First value (index=0) in small_ha.nc is 0.3400644 cm, 655 655 # Fourth value (index==3) is -6.50198 cm 656 656 657 657 658 658 from coordinate_transforms.redfearn import redfearn … … 665 665 lat_index = 0 666 666 lon_index = 2 667 667 668 668 test_value = fid.variables['HA'][:][time_index, lat_index, lon_index] 669 669 test_time = fid.variables['TIME'][:][time_index] 670 670 test_lat = fid.variables['LAT'][:][lat_index] 671 test_lon = fid.variables['LON'][:][lon_index] 672 673 linear_point_index = lat_index*4 + lon_index 674 fid.close() 675 671 test_lon = fid.variables['LON'][:][lon_index] 672 673 linear_point_index = lat_index*4 + lon_index 674 fid.close() 675 676 676 #Call conversion (with zero origin) 677 677 ferret2sww('small', verbose=False, 678 678 origin = (56, 0, 0)) 679 679 680 680 681 681 #Work out the UTM coordinates for test point 682 682 zone, e, n = redfearn(test_lat, test_lon) 683 683 684 684 #Read output file 'small.sww' 685 685 fid = NetCDFFile('small.sww') 686 686 687 687 x = fid.variables['x'][:] 688 y = fid.variables['y'][:] 688 y = fid.variables['y'][:] 689 689 690 690 #Check that test coordinate is correctly represented … … 705 705 706 706 707 707 708 708 def test_sww_extent(self): 709 709 """Not a test, rather a look at the sww format … … 712 712 import time, os 713 713 from Numeric import array, zeros, allclose, Float, concatenate 714 from Scientific.IO.NetCDF import NetCDFFile 715 716 self.domain.filename = 'datatest' + str(id(self)) 714 from Scientific.IO.NetCDF import NetCDFFile 715 716 self.domain.filename = 'datatest' + str(id(self)) 717 717 self.domain.format = 'sww' 718 718 self.domain.smooth = True 719 self.domain.reduction = mean 719 self.domain.reduction = mean 720 720 self.domain.set_datadir('.') 721 721 722 722 723 723 sww = get_dataobject(self.domain) 724 724 sww.store_connectivity() … … 736 736 [xmin, xmax, ymin, ymax, stagemin, stagemax] = \ 737 737 extent_sww(file_and_extension_name ) 738 739 assert allclose(xmin, 0.0) 740 assert allclose(xmax, 1.0) 741 assert allclose(ymin, 0.0) 738 739 assert allclose(xmin, 0.0) 740 assert allclose(xmax, 1.0) 741 assert allclose(ymin, 0.0) 742 742 assert allclose(ymax, 1.0) 743 assert allclose(stagemin, -0.85) 743 assert allclose(stagemin, -0.85) 744 744 assert allclose(stagemax, 0.15) 745 745 … … 747 747 #Cleanup 748 748 os.remove(sww.filename) 749 750 749 750 751 751 def test_ferret2sww_nz_origin(self): 752 752 from Scientific.IO.NetCDF import NetCDFFile 753 from coordinate_transforms.redfearn import redfearn 753 from coordinate_transforms.redfearn import redfearn 754 754 755 755 #Call conversion (with nonzero origin) … … 765 765 766 766 x = fid.variables['x'][:] 767 y = fid.variables['y'][:] 767 y = fid.variables['y'][:] 768 768 769 769 #Check that first coordinate is correctly represented … … 777 777 os.remove('small.sww') 778 778 779 779 780 780 #------------------------------------------------------------- 781 781 if __name__ == "__main__": 782 suite = unittest.makeSuite( dataTestCase,'test')782 suite = unittest.makeSuite(Test_Data_Manager,'test') 783 783 runner = unittest.TextTestRunner() 784 784 runner.run(suite) -
inundation/ga/storm_surge/pyvolution/test_domain.py
r773 r1018 18 18 #print 'bottom - indexes',elements 19 19 domain.set_quantity('friction', 0.09, indexes = elements) 20 20 21 21 def set_top_friction(tag, elements, domain): 22 22 if tag == "top": 23 23 #print 'top - indexes',elements 24 24 domain.set_quantity('friction', 1., indexes = elements) 25 26 25 26 27 27 def set_all_friction(tag, elements, domain): 28 28 if tag == "all": 29 29 new_values = domain.get_quantity('friction', indexes = elements) + 10.0 30 30 31 31 domain.set_quantity('friction', new_values, indexes = elements) 32 33 34 class Test Case(unittest.TestCase):32 33 34 class Test_Domain(unittest.TestCase): 35 35 def setUp(self): 36 36 pass 37 37 38 38 39 39 def tearDown(self): 40 40 pass 41 41 42 42 43 43 def test_simple(self): 44 44 a = [0.0, 0.0] … … 50 50 51 51 points = [a, b, c, d, e, f] 52 #bac, bce, ecf, dbe, daf, dae 52 #bac, bce, ecf, dbe, daf, dae 53 53 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]] 54 54 55 55 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 56 56 other_quantities = ['elevation', 'friction'] 57 57 58 58 domain = Domain(points, vertices, None, 59 conserved_quantities, other_quantities) 59 conserved_quantities, other_quantities) 60 60 domain.check_integrity() 61 61 … … 64 64 65 65 66 assert domain.get_conserved_quantities(0, edge=1) == 0. 66 assert domain.get_conserved_quantities(0, edge=1) == 0. 67 67 68 68 69 69 def test_conserved_quantities(self): 70 71 a = [0.0, 0.0] 72 b = [0.0, 2.0] 73 c = [2.0,0.0] 74 d = [0.0, 4.0] 75 e = [2.0, 2.0] 76 f = [4.0,0.0] 77 78 points = [a, b, c, d, e, f] 79 #bac, bce, ecf, dbe, daf, dae 70 71 a = [0.0, 0.0] 72 b = [0.0, 2.0] 73 c = [2.0,0.0] 74 d = [0.0, 4.0] 75 e = [2.0, 2.0] 76 f = [4.0,0.0] 77 78 points = [a, b, c, d, e, f] 79 #bac, bce, ecf, dbe, daf, dae 80 80 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]] 81 81 … … 87 87 domain.set_quantity('stage', [[1,2,3], [5,5,5], 88 88 [0,0,9], [-6, 3, 3], 89 [0,0,0], [0,0,0]]) 89 [0,0,0], [0,0,0]]) 90 90 91 91 domain.set_quantity('xmomentum', [[1,2,3], [5,5,5], 92 92 [0,0,9], [-6, 3, 3], 93 [0,0,0], [0,0,0]]) 93 [0,0,0], [0,0,0]]) 94 94 95 95 domain.check_integrity() … … 100 100 101 101 q = domain.get_conserved_quantities(1) 102 assert allclose(q, [5., 5., 0.]) 102 assert allclose(q, [5., 5., 0.]) 103 103 104 104 q = domain.get_conserved_quantities(2) … … 106 106 107 107 q = domain.get_conserved_quantities(3) 108 assert allclose(q, [0., 0., 0.]) 109 108 assert allclose(q, [0., 0., 0.]) 109 110 110 111 111 #Edges … … 135 135 assert allclose(q, [-1.5, -1.5, 0.]) 136 136 q = domain.get_conserved_quantities(3, edge=2) 137 assert allclose(q, [-1.5, -1.5, 0.]) 137 assert allclose(q, [-1.5, -1.5, 0.]) 138 138 139 139 … … 142 142 from config import default_boundary_tag 143 143 144 144 145 145 a = [0.0, 0.5] 146 b = [0.0, 0.0] 146 b = [0.0, 0.0] 147 147 c = [0.5, 0.5] 148 148 … … 153 153 domain.set_boundary( {default_boundary_tag: Dirichlet_boundary([5,2,1])} ) 154 154 155 156 domain.check_integrity() 157 158 assert allclose(domain.neighbours, [[-1,-2,-3]]) 159 155 156 domain.check_integrity() 157 158 assert allclose(domain.neighbours, [[-1,-2,-3]]) 159 160 160 161 161 162 162 def test_boundary_conditions(self): 163 163 164 164 a = [0.0, 0.0] 165 165 b = [0.0, 2.0] … … 177 177 (2, 1): 'Second', 178 178 (3, 1): 'Second', 179 (3, 2): 'Second'} 180 179 (3, 2): 'Second'} 180 181 181 182 182 domain = Domain(points, vertices, boundary, 183 183 conserved_quantities =\ 184 ['stage', 'xmomentum', 'ymomentum']) 184 ['stage', 'xmomentum', 'ymomentum']) 185 185 domain.check_integrity() 186 186 … … 205 205 domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5) 206 206 assert domain.quantities['stage'].boundary_values[5] ==\ 207 domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5) 207 domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5) 208 208 209 209 #Check enumeration … … 217 217 """Domain implements a default first order gradient limiter 218 218 """ 219 219 220 220 a = [0.0, 0.0] 221 221 b = [0.0, 2.0] … … 233 233 (2, 1): 'Second', 234 234 (3, 1): 'Second', 235 (3, 2): 'Third'} 236 235 (3, 2): 'Third'} 236 237 237 238 238 domain = Domain(points, vertices, boundary, 239 239 conserved_quantities =\ 240 ['stage', 'xmomentum', 'ymomentum']) 240 ['stage', 'xmomentum', 'ymomentum']) 241 241 domain.check_integrity() 242 242 … … 247 247 assert allclose( domain.quantities['stage'].centroid_values, 248 248 [2,5,3,0] ) 249 249 250 250 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], 251 251 [3,3,3], [4, 4, 4]]) 252 252 253 253 domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], 254 [30,30,30], [40, 40, 40]]) 254 [30,30,30], [40, 40, 40]]) 255 255 256 256 … … 283 283 (2, 1): 'Second', 284 284 (3, 1): 'Second', 285 (3, 2): 'Third'} 286 285 (3, 2): 'Third'} 286 287 287 288 288 domain = Domain(points, vertices, boundary, 289 289 conserved_quantities =\ 290 ['stage', 'xmomentum', 'ymomentum']) 290 ['stage', 'xmomentum', 'ymomentum']) 291 291 domain.check_integrity() 292 292 … … 295 295 domain.set_quantity('xmomentum', [1,2,3,4], 'centroids') 296 296 domain.set_quantity('ymomentum', [1,2,3,4], 'centroids') 297 297 298 298 299 299 #Assign some values to update vectors … … 303 303 domain.quantities[name].explicit_update = array([4.,3.,2.,1.]) 304 304 domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.]) 305 305 306 306 307 307 #Update with given timestep (assuming no other forcing terms) … … 316 316 317 317 for name in domain.conserved_quantities: 318 assert allclose(domain.quantities[name].centroid_values, x) 318 assert allclose(domain.quantities[name].centroid_values, x) 319 319 320 320 … … 322 322 """Set quantities for sub region 323 323 """ 324 324 325 325 a = [0.0, 0.0] 326 326 b = [0.0, 2.0] … … 338 338 (2, 1): 'Second', 339 339 (3, 1): 'Second', 340 (3, 2): 'Third'} 340 (3, 2): 'Third'} 341 341 342 342 domain = Domain(points, vertices, boundary, 343 343 conserved_quantities =\ 344 ['stage', 'xmomentum', 'ymomentum']) 344 ['stage', 'xmomentum', 'ymomentum']) 345 345 domain.check_integrity() 346 346 … … 350 350 assert allclose( domain.quantities['stage'].centroid_values, 351 351 [2,5,3,0] ) 352 352 353 353 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], 354 354 [3,3,3], [4, 4, 4]]) 355 355 356 356 domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], 357 [30,30,30], [40, 40, 40]]) 357 [30,30,30], [40, 40, 40]]) 358 358 359 359 … … 369 369 domain.build_tagged_elements_dictionary({'mound':[0,1]}) 370 370 domain.set_region([add_to_verts]) 371 371 372 372 self.failUnless(domain.test == "Mound", 373 373 'set region failed') 374 374 375 375 376 376 … … 382 382 from shallow_water import Domain 383 383 from Numeric import zeros, Float 384 384 385 385 #Create basic mesh 386 386 points, vertices, boundary = rectangular(1, 3) … … 392 392 'all':[0,1,2,3,4,5]}) 393 393 394 394 395 395 #Set friction 396 396 manning = 0.07 … … 406 406 [ 1.0, 1.0, 1.0], 407 407 [ 1.0, 1.0, 1.0]]) 408 408 409 409 domain.set_region([set_all_friction]) 410 410 #print domain.quantities['friction'].get_values() … … 416 416 [ 11.0, 11.0, 11.0], 417 417 [ 11.0, 11.0, 11.0]]) 418 419 418 420 #------------------------------------------------------------- 419 421 if __name__ == "__main__": 420 suite = unittest.makeSuite(Test Case,'test')422 suite = unittest.makeSuite(Test_Domain,'test') 421 423 runner = unittest.TextTestRunner() 422 424 runner.run(suite) -
inundation/ga/storm_surge/pyvolution/test_general_mesh.py
r1011 r1018 9 9 from Numeric import allclose, array, ones, Float 10 10 11 12 class Test Case(unittest.TestCase):11 12 class Test_General_Mesh(unittest.TestCase): 13 13 def setUp(self): 14 14 pass 15 15 16 16 def tearDown(self): 17 17 pass … … 24 24 from shallow_water import Domain 25 25 from Numeric import zeros, Float 26 26 27 27 #Create basic mesh 28 28 points, vertices, boundary = rectangular(1, 3) 29 29 domain = Domain(points, vertices, boundary) 30 30 31 31 value = [7] 32 32 indexes = [1] 33 assert domain.get_vertices() == domain.triangles 33 assert domain.get_vertices() == domain.triangles 34 34 assert domain.get_vertices([0,4]) == [domain.triangles[0], 35 35 domain.triangles[4]] … … 38 38 from shallow_water import Domain 39 39 from Numeric import zeros, Float 40 40 41 41 #Create basic mesh 42 42 points, vertices, boundary = rectangular(1, 3) 43 43 domain = Domain(points, vertices, boundary) 44 44 45 assert domain.get_area() == 1.0 45 assert domain.get_area() == 1.0 46 46 47 47 … … 53 53 from shallow_water import Domain 54 54 from Numeric import zeros, Float 55 55 56 56 #Create basic mesh 57 57 points, vertices, boundary = rectangular(1, 3) 58 58 domain = Domain(points, vertices, boundary) 59 59 60 60 assert domain.get_unique_vertices() == [0,1,2,3,4,5,6,7] 61 61 unique_vertices = domain.get_unique_vertices([0,1,4]) 62 62 unique_vertices.sort() 63 63 assert unique_vertices == [0,1,2,4,5,6,7] 64 64 65 65 unique_vertices = domain.get_unique_vertices([0,4]) 66 66 unique_vertices.sort() 67 67 assert unique_vertices == [0,2,4,5,6,7] 68 68 69 69 #------------------------------------------------------------- 70 70 if __name__ == "__main__": 71 suite = unittest.makeSuite(Test Case,'test')71 suite = unittest.makeSuite(Test_General_Mesh,'test') 72 72 runner = unittest.TextTestRunner() 73 73 runner.run(suite) 74 74 75 -
inundation/ga/storm_surge/pyvolution/test_generic_boundary_conditions.py
r850 r1018 8 8 from Numeric import allclose, array 9 9 10 11 class Test Case(unittest.TestCase):10 11 class Test_Generic_Boundary_Conditions(unittest.TestCase): 12 12 def setUp(self): 13 13 pass 14 14 #print " Setting up" 15 15 16 16 def tearDown(self): 17 17 pass … … 28 28 else: 29 29 raise 'Should have raised exception' 30 30 31 31 32 32 def test_dirichlet_empty(self): … … 50 50 from domain import Domain 51 51 from quantity import Quantity 52 52 53 53 a = [0.0, 0.0] 54 54 b = [0.0, 2.0] … … 62 62 #bac, bce, ecf, dbe 63 63 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 64 64 65 65 domain = Domain(points, elements) 66 66 domain.check_integrity() 67 67 68 domain.conserved_quantities = ['stage', 'ymomentum'] 68 domain.conserved_quantities = ['stage', 'ymomentum'] 69 69 domain.quantities['stage'] =\ 70 70 Quantity(domain, [[1,2,3], [5,5,5], … … 110 110 from math import sin, pi 111 111 from config import time_format 112 112 113 113 a = [0.0, 0.0] 114 114 b = [0.0, 2.0] … … 122 122 #bac, bce, ecf, dbe 123 123 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 124 124 125 125 domain = Domain(points, elements) 126 126 domain.conserved_quantities = ['stage', 'ymomentum'] … … 134 134 135 135 domain.check_integrity() 136 136 137 137 138 138 #Write file … … 145 145 t_string = time.strftime(time_format, time.gmtime(t)) 146 146 147 fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10))) 147 fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10))) 148 148 fid.close() 149 150 149 150 151 151 F = File_boundary(filename, domain) 152 152 os.remove(filename) 153 153 154 154 155 155 #Check that midpoint coordinates at boundary are correctly computed 156 assert allclose( F.midpoint_coordinates, 157 [[0.0, 3.0], [1.0, 3.0], [0.0, 1.0], 156 assert allclose( F.midpoint_coordinates, 157 [[0.0, 3.0], [1.0, 3.0], [0.0, 1.0], 158 158 [1.0, 0.0], [3.0, 0.0], [3.0, 1.0]]) 159 159 160 160 #assert allclose(F.midpoint_coordinates[(3,2)], [0.0, 3.0]) 161 161 #assert allclose(F.midpoint_coordinates[(3,1)], [1.0, 3.0]) … … 163 163 #assert allclose(F.midpoint_coordinates[(0,0)], [1.0, 0.0]) 164 164 #assert allclose(F.midpoint_coordinates[(2,0)], [3.0, 0.0]) 165 #assert allclose(F.midpoint_coordinates[(2,1)], [3.0, 1.0]) 165 #assert allclose(F.midpoint_coordinates[(2,1)], [3.0, 1.0]) 166 166 167 167 … … 177 177 domain.time = 2.5*5*60 #Half way between steps 2 and 3 178 178 q = F.evaluate() 179 assert allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2]) 179 assert allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2]) 180 180 181 181 … … 191 191 from math import sin, pi 192 192 from config import time_format 193 193 194 194 a = [0.0, 0.0] 195 195 b = [0.0, 2.0] … … 203 203 #bac, bce, ecf, dbe 204 204 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 205 205 206 206 domain = Domain(points, elements) 207 207 domain.conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] … … 228 228 t_string = time.strftime(time_format, time.gmtime(t)) 229 229 230 fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10))) 230 fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10))) 231 231 fid.close() 232 232 233 233 234 234 try: … … 242 242 243 243 244 244 245 245 #------------------------------------------------------------- 246 246 if __name__ == "__main__": 247 suite = unittest.makeSuite(Test Case,'test')247 suite = unittest.makeSuite(Test_Generic_Boundary_Conditions,'test') 248 248 runner = unittest.TextTestRunner() 249 249 runner.run(suite) 250 250 251 251 252 252 -
inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py
r907 r1018 16 16 17 17 18 class dataTestCase(unittest.TestCase):18 class Test_Interpolate_sww(unittest.TestCase): 19 19 def setUp(self): 20 20 import time … … 30 30 domain.beta_h = 0 31 31 32 32 33 33 #Set some field values 34 domain.set_quantity('elevation', lambda x,y: -x) 34 domain.set_quantity('elevation', lambda x,y: -x) 35 35 domain.set_quantity('friction', 0.03) 36 36 37 37 38 ###################### 38 ###################### 39 39 # Boundary conditions 40 40 B = Transmissive_boundary(domain) 41 41 domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B}) 42 43 44 ###################### 42 43 44 ###################### 45 45 #Initial condition - with jumps 46 46 … … 50 50 h = 0.3 51 51 for i in range(stage.shape[0]): 52 if i % 2 == 0: 52 if i % 2 == 0: 53 53 stage[i,:] = bed[i,:] + h 54 54 else: 55 55 stage[i,:] = bed[i,:] 56 56 57 57 domain.set_quantity('stage', stage) 58 58 59 domain.distribute_to_vertices_and_edges() 59 domain.distribute_to_vertices_and_edges() 60 60 61 61 … … 65 65 self.X = C[:,0:6:2].copy() 66 66 self.Y = C[:,1:6:2].copy() 67 67 68 68 self.F = bed 69 69 70 70 71 71 def tearDown(self): 72 72 pass … … 78 78 import time, os 79 79 from Numeric import array, zeros, allclose, Float, concatenate 80 from Scientific.IO.NetCDF import NetCDFFile 80 from Scientific.IO.NetCDF import NetCDFFile 81 81 82 82 self.domain.filename = 'datatest' + str(time.time()) 83 83 self.domain.format = 'sww' 84 84 self.domain.smooth = True 85 self.domain.reduction = mean 86 85 self.domain.reduction = mean 86 87 87 sww = get_dataobject(self.domain) 88 88 sww.store_connectivity() 89 89 sww.store_timestep('stage') 90 self.domain.time = 2. 90 self.domain.time = 2. 91 91 sww.store_timestep('stage') 92 92 93 93 #Check contents 94 94 #Get NetCDF 95 fid = NetCDFFile(sww.filename, 'r') 96 95 fid = NetCDFFile(sww.filename, 'r') 96 97 97 # Get the variables 98 98 x = fid.variables['x'] … … 100 100 z = fid.variables['elevation'] 101 101 102 volumes = fid.variables['volumes'] 102 volumes = fid.variables['volumes'] 103 103 time = fid.variables['time'] 104 104 105 # 2D 106 stage = fid.variables['stage'] 105 # 2D 106 stage = fid.variables['stage'] 107 107 108 108 X = x[:] … … 127 127 print "Stage ",S 128 128 print "****************************" 129 130 131 132 129 130 131 132 133 133 fid.close() 134 134 135 135 #Cleanup 136 136 os.remove(sww.filename) 137 137 138 138 def test_interpolate_sww(self): 139 """Not reaa unit test, rather a system test for 139 """Not reaa unit test, rather a system test for 140 140 """ 141 141 … … 143 143 from Numeric import array, zeros, allclose, Float, concatenate, \ 144 144 transpose 145 from Scientific.IO.NetCDF import NetCDFFile 145 from Scientific.IO.NetCDF import NetCDFFile 146 146 import tempfile 147 147 from load_mesh.loadASCII import load_points_file, \ … … 151 151 self.domain.format = 'sww' 152 152 self.domain.smooth = True 153 self.domain.reduction = mean 154 153 self.domain.reduction = mean 154 155 155 sww = get_dataobject(self.domain) 156 156 sww.store_connectivity() 157 157 sww.store_timestep('stage') 158 self.domain.time = 2. 158 self.domain.time = 2. 159 159 sww.store_timestep('stage') 160 160 161 161 #print "self.domain.filename",self.domain.filename 162 162 interp = Interpolate_sww(sww.filename, 'depth') 163 163 164 164 assert allclose(interp.time,[0.0,2.0]) 165 165 … … 186 186 point_file_out = tempfile.mktemp(".xya") 187 187 interp.write_depth_xya(point_file_out) 188 188 189 189 #check the output file 190 190 xya_dict = load_points_file(point_file_out) 191 191 #print "xya_dict", xya_dict 192 title_list,point_attributes = concatinate_attributelist(xya_dict['attributelist']) 192 title_list,point_attributes = concatinate_attributelist(xya_dict['attributelist']) 193 193 assert allclose(interp.point_coordinates, xya_dict['pointlist']) 194 194 assert allclose(interp.interpolated_quantity, … … 201 201 #time_list = [] 202 202 203 # Try another quantity 203 # Try another quantity 204 204 interp = Interpolate_sww(sww.filename, 'stage') 205 205 interp.interpolate_xya(point_file) … … 209 209 assert allclose(interp.interpolated_quantity,answer) 210 210 211 211 212 212 # look at error catching 213 213 try: … … 218 218 self.failUnless(0==1, 219 219 'bad key did not raise an error!') 220 220 221 221 # look at error catching 222 222 try: … … 227 227 self.failUnless(0==1, 228 228 'bad key did not raise an error!') 229 229 230 230 #Cleanup 231 231 os.remove(sww.filename) 232 232 os.remove(point_file_out) 233 233 os.remove(point_file) 234 234 235 235 #------------------------------------------------------------- 236 236 if __name__ == "__main__": 237 suite = unittest.makeSuite( dataTestCase,'test')237 suite = unittest.makeSuite(Test_Interpolate_sww,'test') 238 238 runner = unittest.TextTestRunner() 239 239 runner.run(suite) -
inundation/ga/storm_surge/pyvolution/test_least_squares.py
r1003 r1018 11 11 12 12 def distance(x, y): 13 return sqrt( sum( (array(x)-array(y))**2 )) 13 return sqrt( sum( (array(x)-array(y))**2 )) 14 14 15 15 def linear_function(point): 16 16 point = array(point) 17 17 return point[:,0]+point[:,1] 18 19 20 class Test Case(unittest.TestCase):18 19 20 class Test_Least_Squares(unittest.TestCase): 21 21 22 22 def setUp(self): 23 23 pass 24 24 25 25 def tearDown(self): 26 26 pass … … 33 33 vertices = [ [1,0,2] ] #bac 34 34 35 data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point 36 35 data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point 36 37 37 interp = Interpolation(points, vertices, data) 38 assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]]) 39 38 assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]]) 39 40 40 41 41 def test_quad_tree(self): … … 50 50 p8 = [-66.0, 20.0] 51 51 p9 = [10.0, -66.0] 52 52 53 53 points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9] 54 54 triangles = [ [0, 1, 2], … … 60 60 [7, 4, 5], 61 61 [8, 0, 2], 62 [0, 9, 1]] 63 64 data = [ [4,4] ] 62 [0, 9, 1]] 63 64 data = [ [4,4] ] 65 65 interp = Interpolation(points, triangles, data, alpha = 0.0, 66 66 max_points_per_cell = 4) … … 75 75 assert allclose(interp.get_A(), answer) 76 76 77 77 78 78 #point outside of quad tree root cell 79 79 interp.set_point_coordinates([[-70, -70]]) … … 82 82 0., 0. , 0., 0., 0., 0.]] 83 83 assert allclose(interp.get_A(), answer) 84 84 85 85 def test_expand_search(self): 86 86 p0 = [-10.0, -10.0] … … 94 94 p8 = [-66.0, 20.0] 95 95 p9 = [10.0, -66.0] 96 96 97 97 points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9] 98 98 triangles = [ [0, 1, 2], … … 104 104 [7, 4, 5], 105 105 [8, 0, 2], 106 [0, 9, 1]] 106 [0, 9, 1]] 107 107 108 108 data = [ [4,4], … … 134 134 [50,50], 135 135 [30,0], 136 [-20,-20]] 136 [-20,-20]] 137 137 point_attributes = [ -400000, 138 138 10, … … 163 163 10, 164 164 10, 165 99] 166 165 99] 166 167 167 interp = Interpolation(points, triangles, data, 168 168 alpha=0.0, expand_search=False, #verbose = True, #False, … … 170 170 calc = interp.fit_points(point_attributes, ) 171 171 #print "calc",calc 172 172 173 173 # the point at 4,4 is ignored. An expanded search has to be done 174 174 # to fine which triangel it's in. … … 194 194 p4 = [10.0, 60.0] 195 195 p5 = [60.0, 60.0] 196 196 197 197 points = [p0, p1, p2, p3, p4, p5] 198 198 triangles = [ [0, 1, 2], … … 201 201 [0, 2, 4], 202 202 [4, 2, 5], 203 [5, 2, 3]] 203 [5, 2, 3]] 204 204 205 205 data = [ [-26.0,-26.0] ] … … 216 216 assert allclose(interp.get_A(), answer) 217 217 218 218 219 219 #point outside of quad tree root cell 220 220 interp.set_point_coordinates([[-70, -70]]) … … 223 223 0., 0. ]] 224 224 assert allclose(interp.get_A(), answer) 225 225 226 226 227 227 def test_datapoints_at_vertices(self): 228 228 """Test that data points coinciding with vertices yield a diagonal matrix 229 229 """ 230 230 231 231 a = [0.0, 0.0] 232 232 b = [0.0, 2.0] … … 235 235 vertices = [ [1,0,2] ] #bac 236 236 237 data = points #Use data at vertices 238 237 data = points #Use data at vertices 238 239 239 interp = Interpolation(points, vertices, data) 240 240 assert allclose(interp.get_A(), [[1., 0., 0.], 241 241 [0., 1., 0.], 242 [0., 0., 1.]]) 243 242 [0., 0., 1.]]) 243 244 244 245 245 … … 248 248 each point should affect two matrix entries equally 249 249 """ 250 250 251 251 a = [0.0, 0.0] 252 252 b = [0.0, 2.0] … … 256 256 257 257 data = [ [0., 1.], [1., 0.], [1., 1.] ] 258 258 259 259 interp = Interpolation(points, vertices, data) 260 261 assert allclose(interp.get_A(), [[0.5, 0.5, 0.0], #Affects vertex 1 and 0 260 261 assert allclose(interp.get_A(), [[0.5, 0.5, 0.0], #Affects vertex 1 and 0 262 262 [0.5, 0.0, 0.5], #Affects vertex 0 and 2 263 263 [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2 … … 268 268 each point should affect two matrix entries in proportion 269 269 """ 270 270 271 271 a = [0.0, 0.0] 272 272 b = [0.0, 2.0] … … 276 276 277 277 data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ] 278 278 279 279 interp = Interpolation(points, vertices, data) 280 280 281 281 assert allclose(interp.get_A(), [[0.25, 0.75, 0.0], #Affects vertex 1 and 0 282 282 [0.25, 0.0, 0.75], #Affects vertex 0 and 2 … … 288 288 289 289 from Numeric import sum 290 290 291 291 a = [0.0, 0.0] 292 292 b = [0.0, 2.0] … … 296 296 297 297 data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ] 298 298 299 299 interp = Interpolation(points, vertices, data) 300 300 #print "interp.get_A()", interp.get_A() 301 301 assert allclose(sum(interp.get_A(), axis=1), 1.0) 302 302 303 303 def test_arbitrary_datapoints_some_outside(self): 304 304 """Try arbitrary datapoints one outside the triangle. … … 307 307 308 308 from Numeric import sum 309 309 310 310 a = [0.0, 0.0] 311 311 b = [0.0, 2.0] … … 322 322 interp = Interpolation(points, vertices, data, precrop = False) 323 323 assert allclose(sum(interp.get_A(), axis=1), [1,1,1,0]) 324 325 326 327 # this causes a memory error in scipy.sparse 324 325 326 327 # this causes a memory error in scipy.sparse 328 328 def test_more_triangles(self): 329 329 330 330 a = [-1.0, 0.0] 331 331 b = [3.0, 4.0] … … 340 340 #Data points 341 341 data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ] 342 interp = Interpolation(points, triangles, data_points) 343 344 answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d 342 interp = Interpolation(points, triangles, data_points) 343 344 answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d 345 345 [0.5, 0.0, 0.0, 0.5, 0.0, 0.0], #Affects points a and d 346 346 [0.75, 0.25, 0.0, 0.0, 0.0, 0.0], #Affects points a and b … … 354 354 for j in range(A.shape[1]): 355 355 if not allclose(A[i,j], answer[i][j]): 356 print i,j,':',A[i,j], answer[i][j] 356 print i,j,':',A[i,j], answer[i][j] 357 357 358 358 … … 368 368 points = [a, b, c] 369 369 triangles = [ [1,0,2] ] #bac 370 370 371 371 d1 = [1.0, 1.0] 372 372 d2 = [1.0, 3.0] … … 375 375 z2 = 4 376 376 z3 = 4 377 data_coords = [d1, d2, d3] 378 377 data_coords = [d1, d2, d3] 378 379 379 interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20) 380 380 z = [z1, z2, z3] … … 384 384 #print "f\n",f 385 385 #print "answer\n",answer 386 386 387 387 assert allclose(f, answer, atol=1e-7) 388 389 388 389 390 390 def test_smooth_att_to_meshII(self): 391 391 392 392 a = [0.0, 0.0] 393 393 b = [0.0, 5.0] … … 395 395 points = [a, b, c] 396 396 triangles = [ [1,0,2] ] #bac 397 397 398 398 d1 = [1.0, 1.0] 399 399 d2 = [1.0, 2.0] 400 400 d3 = [3.0,1.0] 401 data_coords = [d1, d2, d3] 401 data_coords = [d1, d2, d3] 402 402 z = linear_function(data_coords) 403 403 interp = Interpolation(points, triangles, data_coords, alpha=0.0) 404 404 f = interp.fit(z) 405 405 answer = linear_function(points) 406 406 407 407 assert allclose(f, answer) 408 408 409 409 def test_smooth_attributes_to_meshIII(self): 410 410 411 411 a = [-1.0, 0.0] 412 412 b = [3.0, 4.0] … … 441 441 #print "answer\n",answer 442 442 assert allclose(f, answer) 443 443 444 444 445 445 def test_smooth_attributes_to_meshIV(self): 446 446 """ Testing 2 attributes smoothed to the mesh 447 447 """ 448 448 449 449 a = [0.0, 0.0] 450 450 b = [0.0, 5.0] … … 452 452 points = [a, b, c] 453 453 triangles = [ [1,0,2] ] #bac 454 454 455 455 d1 = [1.0, 1.0] 456 456 d2 = [1.0, 3.0] … … 459 459 z2 = [4, 8] 460 460 z3 = [4, 8] 461 data_coords = [d1, d2, d3] 462 461 data_coords = [d1, d2, d3] 462 463 463 interp = Interpolation(points, triangles, data_coords, alpha=0.0) 464 464 z = [z1, z2, z3] … … 466 466 answer = [[0,0], [5., 10.], [5., 10.]] 467 467 assert allclose(f, answer) 468 468 469 469 def test_interpolate_attributes_to_points(self): 470 470 v0 = [0.0, 0.0] 471 471 v1 = [0.0, 5.0] 472 472 v2 = [5.0, 0.0] 473 473 474 474 vertices = [v0, v1, v2] 475 475 triangles = [ [1,0,2] ] #bac 476 476 477 477 d0 = [1.0, 1.0] 478 478 d1 = [1.0, 2.0] 479 479 d2 = [3.0, 1.0] 480 480 point_coords = [ d0, d1, d2] 481 481 482 482 interp = Interpolation(vertices, triangles, point_coords) 483 483 f = linear_function(vertices) 484 484 z = interp.interpolate(f) 485 485 answer = linear_function(point_coords) 486 486 487 487 488 488 assert allclose(z, answer) 489 489 490 490 491 491 def test_interpolate_attributes_to_pointsII(self): 492 492 a = [-1.0, 0.0] … … 500 500 triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc 501 501 502 502 503 503 point_coords = [[-2.0, 2.0], 504 504 [-1.0, 1.0], … … 513 513 [0.5, -1.9], 514 514 [3.0, 1.0]] 515 515 516 516 interp = Interpolation(vertices, triangles, point_coords) 517 517 f = linear_function(vertices) … … 521 521 #print "answer",answer 522 522 assert allclose(z, answer) 523 523 524 524 def test_interpolate_attributes_to_pointsIII(self): 525 525 """Test linear interpolation of known values at vertices to … … 530 530 c = [5.0, 0.0] 531 531 d = [5.0, 5.0] 532 532 533 533 vertices = [a, b, c, d] 534 534 triangles = [ [1,0,2], [2,3,0] ] #bac, cdb … … 545 545 d4 = [2.5, 2.5] 546 546 d5 = [4.0, 1.0] 547 547 548 548 #Point on common vertex 549 549 d6 = [0., 5.] 550 550 551 551 552 552 point_coords = [d0, d1, d2, d3, d4, d5, d6] 553 553 554 554 interp = Interpolation(vertices, triangles, point_coords) 555 555 556 556 #Known values at vertices 557 #Functions are x+y, x+2y, 2x+y, x-y-5 558 f = [ [0., 0., 0., -5.], # (0,0) 557 #Functions are x+y, x+2y, 2x+y, x-y-5 558 f = [ [0., 0., 0., -5.], # (0,0) 559 559 [5., 10., 5., -10.], # (0,5) 560 560 [5., 5., 10.0, 0.], # (5,0) 561 561 [10., 15., 15., -5.]] # (5,5) 562 562 563 563 z = interp.interpolate(f) 564 564 answer = [ [2., 3., 3., -5.], # (1,1) … … 576 576 577 577 assert allclose(z, answer) 578 578 579 579 def test_interpolate_attributes_to_pointsIV(self): 580 580 a = [-1.0, 0.0] … … 588 588 triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc 589 589 590 590 591 591 point_coords = [[-2.0, 2.0], 592 592 [-1.0, 1.0], … … 601 601 [0.5, -1.9], 602 602 [3.0, 1.0]] 603 603 604 604 interp = Interpolation(vertices, triangles, point_coords) 605 605 f = array([linear_function(vertices),2*linear_function(vertices) ]) … … 607 607 #print "f",f 608 608 z = interp.interpolate(f) 609 answer = [linear_function(point_coords), 609 answer = [linear_function(point_coords), 610 610 2*linear_function(point_coords) ] 611 611 answer = transpose(answer) … … 613 613 #print "answer",answer 614 614 assert allclose(z, answer) 615 615 616 616 def test_smooth_attributes_to_mesh_function(self): 617 617 """ Testing 2 attributes smoothed to the mesh 618 618 """ 619 619 620 620 a = [0.0, 0.0] 621 621 b = [0.0, 5.0] … … 623 623 points = [a, b, c] 624 624 triangles = [ [1,0,2] ] #bac 625 625 626 626 d1 = [1.0, 1.0] 627 627 d2 = [1.0, 3.0] … … 632 632 data_coords = [d1, d2, d3] 633 633 z = [z1, z2, z3] 634 634 635 635 f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0) 636 636 answer = [[0, 0], [5., 10.], [5., 10.]] 637 637 638 638 assert allclose(f, answer) 639 639 640 640 641 641 642 642 def test_pts2rectangular(self): 643 643 … … 645 645 FN = 'xyatest' + str(time.time()) + '.xya' 646 646 fid = open(FN, 'w') 647 fid.write(' %s \n' %('elevation')) 647 fid.write(' %s \n' %('elevation')) 648 648 fid.write('%f %f %f\n' %(1,1,2) ) 649 649 fid.write('%f %f %f\n' %(1,3,4) ) 650 fid.write('%f %f %f\n' %(3,1,4) ) 650 fid.write('%f %f %f\n' %(3,1,4) ) 651 651 fid.close() 652 652 653 653 points, triangles, boundary, attributes =\ 654 654 pts2rectangular(FN, 4, 4, format = 'asc') 655 655 656 656 657 657 data_coords = [ [1,1], [1,3], [3,1] ] 658 658 z = [2, 4, 4] 659 659 660 ref = fit_to_mesh(points, triangles, data_coords, z) 660 ref = fit_to_mesh(points, triangles, data_coords, z) 661 661 662 662 #print attributes … … 665 665 666 666 os.remove(FN) 667 668 669 #Tests of smoothing matrix 667 668 669 #Tests of smoothing matrix 670 670 def test_smoothing_matrix_one_triangle(self): 671 671 from Numeric import dot … … 674 674 c = [2.0,0.0] 675 675 points = [a, b, c] 676 676 677 677 vertices = [ [1,0,2] ] #bac 678 678 … … 689 689 # int 1 dx dy = area = 2 690 690 assert dot(dot(f, interp.get_D()), f) == 2 691 691 692 692 #Define f(x,y) = y 693 693 f = array([0,2,0]) #Value at global vertex 1 … … 702 702 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 703 703 # int 2 dx dy = 2*area = 4 704 assert dot(dot(f, interp.get_D()), f) == 4 705 704 assert dot(dot(f, interp.get_D()), f) == 4 705 706 706 707 707 708 708 def test_smoothing_matrix_more_triangles(self): 709 709 from Numeric import dot 710 710 711 711 a = [0.0, 0.0] 712 712 b = [0.0, 2.0] … … 717 717 718 718 points = [a, b, c, d, e, f] 719 #bac, bce, ecf, dbe, daf, dae 719 #bac, bce, ecf, dbe, daf, dae 720 720 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 721 721 … … 728 728 729 729 #Define f(x,y) = x 730 f = array([0,0,2,0,2,4]) #f evaluated at points a-f 730 f = array([0,0,2,0,2,4]) #f evaluated at points a-f 731 731 732 732 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 733 733 # int 1 dx dy = total area = 8 734 734 assert dot(dot(f, interp.get_D()), f) == 8 735 735 736 736 #Define f(x,y) = y 737 f = array([0,2,0,4,2,0]) #f evaluated at points a-f 737 f = array([0,2,0,4,2,0]) #f evaluated at points a-f 738 738 739 739 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = … … 742 742 743 743 #Define f(x,y) = x+y 744 f = array([0,2,2,4,4,4]) #f evaluated at points a-f 744 f = array([0,2,2,4,4,4]) #f evaluated at points a-f 745 745 746 746 #Check that int (df/dx)**2 + (df/dy)**2 dx dy = 747 747 # int 2 dx dy = 2*area = 16 748 assert dot(dot(f, interp.get_D()), f) == 16 748 assert dot(dot(f, interp.get_D()), f) == 16 749 749 750 750 751 751 def test_fit_and_interpolation(self): 752 752 from mesh import Mesh 753 753 754 754 a = [0.0, 0.0] 755 755 b = [0.0, 2.0] … … 760 760 761 761 points = [a, b, c, d, e, f] 762 #bac, bce, ecf, dbe, daf, dae 762 #bac, bce, ecf, dbe, daf, dae 763 763 triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 764 764 … … 773 773 [ 1.0, 1.0], 774 774 [ 1.0, 2.0], 775 [ 1.0, 3.0], 775 [ 1.0, 3.0], 776 776 [ 2.0, 1.0], 777 777 [ 3.0, 0.0], 778 [ 3.0, 1.0]] 779 778 [ 3.0, 1.0]] 779 780 780 interp = Interpolation(points, triangles, data_points, alpha=0.0) 781 781 … … 797 797 798 798 def test_smoothing_and_interpolation(self): 799 799 800 800 a = [0.0, 0.0] 801 801 b = [0.0, 2.0] … … 806 806 807 807 points = [a, b, c, d, e, f] 808 #bac, bce, ecf, dbe, daf, dae 808 #bac, bce, ecf, dbe, daf, dae 809 809 triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 810 810 … … 816 816 817 817 z = linear_function(data_points) 818 answer = linear_function(points) 818 answer = linear_function(points) 819 819 820 820 #Make interpolator with too few data points and no smoothing 821 821 interp = Interpolation(points, triangles, data_points, alpha=0.0) 822 #Must raise an exception 822 #Must raise an exception 823 823 try: 824 824 f = interp.fit(z) … … 828 828 #Now try with smoothing parameter 829 829 interp = Interpolation(points, triangles, data_points, alpha=1.0e-13) 830 830 831 831 f = interp.fit(z) 832 832 #f will be different from answerr due to smoothing … … 855 855 856 856 points = [a, b, c, d, e, f] 857 #bac, bce, ecf, dbe, daf, dae 857 #bac, bce, ecf, dbe, daf, dae 858 858 triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 859 859 … … 869 869 [ 15, -17], #Outside mesh 870 870 [ 1.0, 2.0], 871 [ 1.0, 3.0], 871 [ 1.0, 3.0], 872 872 [ 2.0, 1.0], 873 873 [ 3.0, 0.0], … … 889 889 [ 1.0, 0.5], 890 890 [ 2.0, 0.4], 891 [ 2.8, 1.2]] 892 891 [ 2.8, 1.2]] 892 893 893 894 894 … … 907 907 import Numeric 908 908 z1 = Numeric.take(z1, [0,1,2,4,5,6]) 909 answer = Numeric.take(answer, [0,1,2,4,5,6]) 909 answer = Numeric.take(answer, [0,1,2,4,5,6]) 910 910 assert allclose(z1, answer) 911 911 … … 942 942 943 943 points = [a, b, c, d, e, f] 944 #bac, bce, ecf, dbe, daf, dae 944 #bac, bce, ecf, dbe, daf, dae 945 945 triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 946 946 … … 955 955 [ 1.0, 1.0], 956 956 [ 1.0, 2.0], 957 [ 1.0, 3.0], 957 [ 1.0, 3.0], 958 958 [ 2.0, 1.0], 959 959 [ 3.0, 0.0], … … 963 963 #First check that things are OK when using same origin 964 964 mesh_origin = (56, 290000, 618000) #zone, easting, northing 965 data_origin = (56, 290000, 618000) #zone, easting, northing 965 data_origin = (56, 290000, 618000) #zone, easting, northing 966 966 967 967 … … 971 971 data_origin = data_origin, 972 972 mesh_origin = mesh_origin) 973 973 974 974 z = linear_function(data_points1) #Example z-values 975 975 f = interp.fit(z) #Fitted values at vertices … … 982 982 [ 1.0, 0.5], 983 983 [ 2.0, 0.4], 984 [ 2.8, 1.2]] 985 984 [ 2.8, 1.2]] 985 986 986 987 987 #Build new A matrix based on new points … … 1000 1000 #Then check situtaion where points are relative to a different 1001 1001 #origin (same zone, though, until we figure that out (FIXME)) 1002 1002 1003 1003 mesh_origin = (56, 290000, 618000) #zone, easting, northing 1004 data_origin = (56, 10000, 10000) #zone, easting, northing 1004 data_origin = (56, 10000, 10000) #zone, easting, northing 1005 1005 1006 1006 #Shift datapoints according to new origin … … 1012 1012 for k in range(len(data_points2)): 1013 1013 data_points2[k][0] += mesh_origin[1] - data_origin[1] 1014 data_points2[k][1] += mesh_origin[2] - data_origin[2] 1015 1016 1017 1014 data_points2[k][1] += mesh_origin[2] - data_origin[2] 1015 1016 1017 1018 1018 #Fit surface to mesh 1019 1019 interp = Interpolation(points, triangles, data_points1, … … 1021 1021 data_origin = data_origin, 1022 1022 mesh_origin = mesh_origin) 1023 1023 1024 1024 f1 = interp.fit(z) #Fitted values at vertices (using same z as before) 1025 1025 … … 1055 1055 def test_fit_to_mesh_file(self): 1056 1056 from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \ 1057 export_mesh_file 1057 export_mesh_file 1058 1058 import tempfile 1059 1059 import os 1060 1060 1061 1061 # create a .tsh file, no user outline 1062 1062 mesh_dic = {} … … 1072 1072 mesh_dic['segment_tags'] = ['external', 1073 1073 'external', 1074 'external'] 1074 'external'] 1075 1075 mesh_file = tempfile.mktemp(".tsh") 1076 1076 export_mesh_file(mesh_file,mesh_dic) 1077 1077 1078 1078 # create an .xya file 1079 1079 point_file = tempfile.mktemp(".xya") … … 1081 1081 fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n") 1082 1082 fd.close() 1083 1083 1084 1084 mesh_output_file = "new_trianlge.tsh" 1085 1085 fit_to_mesh_file(mesh_file, 1086 1086 point_file, 1087 1087 mesh_output_file, 1088 alpha = 0.0) 1088 alpha = 0.0) 1089 1089 # load in the .tsh file we just wrote 1090 1090 mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file) 1091 #print "mesh_dic",mesh_dic 1091 #print "mesh_dic",mesh_dic 1092 1092 ans =[[0.0, 0.0], 1093 1093 [5.0, 10.0], 1094 1094 [5.0,10.0]] 1095 1095 assert allclose(mesh_dic['vertex_attributes'],ans) 1096 1096 1097 1097 self.failUnless(mesh_dic['vertex_attribute_titles'] == 1098 1098 ['elevation','stage'], 1099 1099 'test_fit_to_mesh_file failed') 1100 1100 1101 1101 #clean up 1102 1102 os.remove(mesh_file) 1103 1103 os.remove(point_file) 1104 os.remove(mesh_output_file) 1105 1104 os.remove(mesh_output_file) 1105 1106 1106 def test_fit_to_mesh_fileII(self): 1107 1107 from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \ 1108 export_mesh_file 1108 export_mesh_file 1109 1109 import tempfile 1110 1110 import os 1111 1111 1112 1112 # create a .tsh file, no user outline 1113 1113 mesh_dic = {} … … 1123 1123 mesh_dic['segment_tags'] = ['external', 1124 1124 'external', 1125 'external'] 1125 'external'] 1126 1126 mesh_file = tempfile.mktemp(".tsh") 1127 1127 export_mesh_file(mesh_file,mesh_dic) 1128 1128 1129 1129 # create an .xya file 1130 1130 point_file = tempfile.mktemp(".xya") … … 1132 1132 fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n") 1133 1133 fd.close() 1134 1134 1135 1135 mesh_output_file = "new_triangle.tsh" 1136 1136 fit_to_mesh_file(mesh_file, 1137 1137 point_file, 1138 1138 mesh_output_file, 1139 alpha = 0.0) 1139 alpha = 0.0) 1140 1140 # load in the .tsh file we just wrote 1141 1141 mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file) … … 1145 1145 [1.0, 2.0,5.0, 10.0], 1146 1146 [1.0, 2.0,5.0,10.0]]) 1147 1147 1148 1148 self.failUnless(mesh_dic['vertex_attribute_titles'] == 1149 1149 ['density', 'temp','elevation','stage'], 1150 1150 'test_fit_to_mesh_file failed') 1151 1151 1152 1152 #clean up 1153 1153 os.remove(mesh_file) 1154 os.remove(mesh_output_file) 1154 os.remove(mesh_output_file) 1155 1155 os.remove(point_file) 1156 1156 1157 1157 def test_fit_to_msh_netcdf_fileII(self): 1158 from load_mesh.loadASCII import mesh_file_to_mesh_dictionary,export_mesh_file 1158 from load_mesh.loadASCII import mesh_file_to_mesh_dictionary,export_mesh_file 1159 1159 import tempfile 1160 1160 import os 1161 1161 1162 1162 # create a .tsh file, no user outline 1163 1163 mesh_dic = {} … … 1173 1173 mesh_dic['segment_tags'] = ['external', 1174 1174 'external', 1175 'external'] 1175 'external'] 1176 1176 mesh_file = tempfile.mktemp(".msh") 1177 1177 export_mesh_file(mesh_file,mesh_dic) 1178 1178 1179 1179 # create an .xya file 1180 1180 point_file = tempfile.mktemp(".xya") … … 1182 1182 fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n") 1183 1183 fd.close() 1184 1184 1185 1185 mesh_output_file = "new_triangle.msh" 1186 1186 fit_to_mesh_file(mesh_file, 1187 1187 point_file, 1188 1188 mesh_output_file, 1189 alpha = 0.0) 1189 alpha = 0.0) 1190 1190 # load in the .tsh file we just wrote 1191 1191 mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file) … … 1195 1195 [1.0, 2.0,5.0, 10.0], 1196 1196 [1.0, 2.0,5.0,10.0]]) 1197 1197 1198 1198 self.failUnless(mesh_dic['vertex_attribute_titles'] == 1199 1199 ['density', 'temp','elevation','stage'], 1200 1200 'test_fit_to_mesh_file failed') 1201 1201 1202 1202 #clean up 1203 1203 os.remove(mesh_file) 1204 os.remove(mesh_output_file) 1204 os.remove(mesh_output_file) 1205 1205 os.remove(point_file) 1206 1206 1207 1207 #------------------------------------------------------------- 1208 1208 if __name__ == "__main__": 1209 suite = unittest.makeSuite(Test Case,'test')1210 1211 #suite = unittest.makeSuite(Test Case,'test_fit_to_msh_netcdf_fileII')1212 #suite = unittest.makeSuite(Test Case,'test_fit_to_mesh_fileII')1209 suite = unittest.makeSuite(Test_Least_Squares,'test') 1210 1211 #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII') 1212 #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII') 1213 1213 runner = unittest.TextTestRunner(verbosity=1) 1214 1214 runner.run(suite) 1215 1215 1216 1217 1218 1219 1220 1216 1217 1218 1219 -
inundation/ga/storm_surge/pyvolution/test_mesh.py
r984 r1018 12 12 13 13 def distance(x, y): 14 return sqrt( sum( (array(x)-array(y))**2 )) 15 16 class Test Case(unittest.TestCase):14 return sqrt( sum( (array(x)-array(y))**2 )) 15 16 class Test_Mesh(unittest.TestCase): 17 17 def setUp(self): 18 18 pass 19 19 20 20 def tearDown(self): 21 21 pass … … 32 32 msg = 'Should have raised exception' 33 33 raise msg 34 35 34 35 36 36 def test_basic_triangle(self): 37 37 … … 39 39 b = [4.0, 0.0] 40 40 c = [0.0, 3.0] 41 41 42 42 points = [a, b, c] 43 43 vertices = [[0,1,2]] 44 mesh = Mesh(points, vertices) 44 mesh = Mesh(points, vertices) 45 45 46 46 #Centroid … … 57 57 assert allclose(normals[0, 0:2], [3.0/5, 4.0/5]) 58 58 assert allclose(normals[0, 2:4], [-1.0, 0.0]) 59 assert allclose(normals[0, 4:6], [0.0, -1.0]) 59 assert allclose(normals[0, 4:6], [0.0, -1.0]) 60 60 61 61 assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5]) 62 62 assert allclose(mesh.get_normal(0,1), [-1.0, 0.0]) 63 assert allclose(mesh.get_normal(0,2), [0.0, -1.0]) 64 63 assert allclose(mesh.get_normal(0,2), [0.0, -1.0]) 64 65 65 #Edge lengths 66 66 assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0]) … … 78 78 79 79 V2 = mesh.get_vertex_coordinate(0, 2) 80 assert allclose(V2, [0.0, 3.0]) 80 assert allclose(V2, [0.0, 3.0]) 81 81 82 82 … … 86 86 mesh.check_integrity() 87 87 88 88 89 89 #Test that the centroid is located 2/3 of the way 90 90 #from each vertex to the midpoint of the opposite side 91 91 92 92 V = mesh.get_vertex_coordinates() 93 93 94 94 x0 = V[0,0] 95 95 y0 = V[0,1] … … 98 98 x2 = V[0,4] 99 99 y2 = V[0,5] 100 100 101 101 m0 = [(x1 + x2)/2, (y1 + y2)/2] 102 102 m1 = [(x0 + x2)/2, (y0 + y2)/2] … … 108 108 # 109 109 d0 = distance(centroid, [x1, y1]) 110 d1 = distance(m1, [x1, y1]) 110 d1 = distance(m1, [x1, y1]) 111 111 assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3) 112 112 113 113 d0 = distance(centroid, [x2, y2]) 114 114 d1 = distance(m2, [x2, y2]) 115 assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3) 116 115 assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3) 116 117 117 #Radius 118 118 d0 = distance(centroid, m0) 119 119 assert d0 == 5.0/6 120 120 121 121 d1 = distance(centroid, m1) 122 assert d1 == sqrt(73.0/36) 123 124 d2 = distance(centroid, m2) 125 assert d2 == sqrt(13.0/9) 126 122 assert d1 == sqrt(73.0/36) 123 124 d2 = distance(centroid, m2) 125 assert d2 == sqrt(13.0/9) 126 127 127 assert mesh.radii[0] == min(d0, d1, d2) 128 128 assert mesh.radii[0] == 5.0/6 … … 134 134 vertices = [[0,3,2], [2,3,1], [1,3,0]] 135 135 new_mesh = Mesh(points, vertices) 136 136 137 137 assert new_mesh.areas[0] == new_mesh.areas[1] 138 assert new_mesh.areas[1] == new_mesh.areas[2] 139 assert new_mesh.areas[1] == new_mesh.areas[2] 140 141 assert new_mesh.areas[1] == mesh.areas[0]/3 138 assert new_mesh.areas[1] == new_mesh.areas[2] 139 assert new_mesh.areas[1] == new_mesh.areas[2] 140 141 assert new_mesh.areas[1] == mesh.areas[0]/3 142 142 143 143 … … 151 151 vertices = [[0,1,2]] 152 152 153 mesh = Mesh(points, vertices) 153 mesh = Mesh(points, vertices) 154 154 centroid = mesh.centroid_coordinates[0] 155 155 156 156 157 157 #Test that the centroid is located 2/3 of the way 158 158 #from each vertex to the midpoint of the opposite side 159 159 160 160 V = mesh.get_vertex_coordinates() 161 161 162 162 x0 = V[0,0] 163 163 y0 = V[0,1] … … 166 166 x2 = V[0,4] 167 167 y2 = V[0,5] 168 168 169 169 m0 = [(x1 + x2)/2, (y1 + y2)/2] 170 170 m1 = [(x0 + x2)/2, (y0 + y2)/2] … … 176 176 # 177 177 d0 = distance(centroid, [x1, y1]) 178 d1 = distance(m1, [x1, y1]) 178 d1 = distance(m1, [x1, y1]) 179 179 assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3) 180 180 181 181 d0 = distance(centroid, [x2, y2]) 182 182 d1 = distance(m2, [x2, y2]) 183 assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3) 184 183 assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3) 184 185 185 #Radius 186 186 d0 = distance(centroid, m0) 187 187 d1 = distance(centroid, m1) 188 d2 = distance(centroid, m2) 188 d2 = distance(centroid, m2) 189 189 assert mesh.radii[0] == min(d0, d1, d2) 190 190 … … 197 197 vertices = [[0,3,2], [2,3,1], [1,3,0]] 198 198 new_mesh = Mesh(points, vertices) 199 199 200 200 assert new_mesh.areas[0] == new_mesh.areas[1] 201 assert new_mesh.areas[1] == new_mesh.areas[2] 202 assert new_mesh.areas[1] == new_mesh.areas[2] 203 204 assert new_mesh.areas[1] == mesh.areas[0]/3 205 201 assert new_mesh.areas[1] == new_mesh.areas[2] 202 assert new_mesh.areas[1] == new_mesh.areas[2] 203 204 assert new_mesh.areas[1] == mesh.areas[0]/3 205 206 206 207 207 #Test that points are arranged in a counter clock wise order … … 217 217 e = [2.0, 2.0] 218 218 points = [a, b, c, e] 219 vertices = [ [1,0,2], [1,2,3] ] #bac, bce 219 vertices = [ [1,0,2], [1,2,3] ] #bac, bce 220 220 mesh = Mesh(points, vertices) 221 221 222 222 assert mesh.areas[0] == 2.0 223 223 224 assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3]) 224 assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3]) 225 225 226 226 227 227 #Test that points are arranged in a counter clock wise order 228 mesh.check_integrity() 229 230 231 228 mesh.check_integrity() 229 230 231 232 232 def test_more_triangles(self): 233 233 234 234 a = [0.0, 0.0] 235 235 b = [0.0, 2.0] … … 240 240 241 241 points = [a, b, c, d, e, f] 242 #bac, bce, ecf, dbe, daf, dae 242 #bac, bce, ecf, dbe, daf, dae 243 243 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]] 244 mesh = Mesh(points, vertices) 244 mesh = Mesh(points, vertices) 245 245 246 246 #Test that points are arranged in a counter clock wise order 247 mesh.check_integrity() 247 mesh.check_integrity() 248 248 249 249 assert mesh.areas[0] == 2.0 … … 252 252 assert mesh.areas[3] == 2.0 253 253 assert mesh.areas[4] == 8.0 254 assert mesh.areas[5] == 4.0 255 254 assert mesh.areas[5] == 4.0 255 256 256 assert mesh.edgelengths[1,0] == 2.0 257 257 assert mesh.edgelengths[1,1] == 2.0 … … 259 259 260 260 assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3]) 261 assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3]) 261 assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3]) 262 262 assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3]) 263 assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3]) 263 assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3]) 264 264 265 265 … … 278 278 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 279 279 mesh = Mesh(points, vertices) 280 280 281 281 mesh.check_integrity() 282 282 … … 284 284 T = mesh 285 285 tid = 0 286 assert T.number_of_boundaries[tid] == 2 286 assert T.number_of_boundaries[tid] == 2 287 287 assert T.neighbours[tid, 0] < 0 #Opposite point b (0,2) 288 288 assert T.neighbours[tid, 1] == 1 #Opposite point a (0,0) 289 assert T.neighbours[tid, 2] < 0 #Opposite point c (2,0) 289 assert T.neighbours[tid, 2] < 0 #Opposite point c (2,0) 290 290 291 291 tid = 1 292 assert T.number_of_boundaries[tid] == 0 292 assert T.number_of_boundaries[tid] == 0 293 293 assert T.neighbours[tid, 0] == 2 #Opposite point b (0,2) 294 294 assert T.neighbours[tid, 1] == 3 #Opposite point c (2,0) … … 296 296 297 297 tid = 2 298 assert T.number_of_boundaries[tid] == 2 298 assert T.number_of_boundaries[tid] == 2 299 299 assert T.neighbours[tid, 0] < 0 #Opposite point e (2,2) 300 300 assert T.neighbours[tid, 1] < 0 #Opposite point c (2,0) 301 assert T.neighbours[tid, 2] == 1 #Opposite point f (4,0) 301 assert T.neighbours[tid, 2] == 1 #Opposite point f (4,0) 302 302 303 303 tid = 3 304 assert T.number_of_boundaries[tid] == 2 304 assert T.number_of_boundaries[tid] == 2 305 305 assert T.neighbours[tid, 0] == 1 #Opposite point d (0,4) 306 306 assert T.neighbours[tid, 1] < 0 #Opposite point b (0,3) … … 311 311 assert T.neighbour_edges[tid, 0] < 0 #Opposite point b (0,2) 312 312 assert T.neighbour_edges[tid, 1] == 2 #Opposite point a (0,0) 313 assert T.neighbour_edges[tid, 2] < 0 #Opposite point c (2,0) 314 315 tid = 1 313 assert T.neighbour_edges[tid, 2] < 0 #Opposite point c (2,0) 314 315 tid = 1 316 316 assert T.neighbour_edges[tid, 0] == 2 #Opposite point b (0,2) 317 317 assert T.neighbour_edges[tid, 1] == 0 #Opposite point c (2,0) 318 318 assert T.neighbour_edges[tid, 2] == 1 #Opposite point e (2,2) 319 319 320 tid = 2 320 tid = 2 321 321 assert T.neighbour_edges[tid, 0] < 0 #Opposite point e (2,2) 322 322 assert T.neighbour_edges[tid, 1] < 0 #Opposite point c (2,0) 323 assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0) 324 325 tid = 3 323 assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0) 324 325 tid = 3 326 326 assert T.neighbour_edges[tid, 0] == 1 #Opposite point d (0,4) 327 327 assert T.neighbour_edges[tid, 1] < 0 #Opposite point b (0,3) 328 328 assert T.neighbour_edges[tid, 2] < 0 #Opposite point e (2,2) 329 329 330 330 331 331 def test_rectangular_mesh_basic(self): 332 332 M=1 … … 336 336 mesh = Mesh(points, vertices, boundary) 337 337 338 #Test that points are arranged in a counter clock wise order 338 #Test that points are arranged in a counter clock wise order 339 339 mesh.check_integrity() 340 340 … … 343 343 points, vertices, boundary = rectangular(M, N) 344 344 mesh = Mesh(points, vertices, boundary) 345 346 #Test that points are arranged in a counter clock wise order 345 346 #Test that points are arranged in a counter clock wise order 347 347 mesh.check_integrity() 348 348 349 349 #assert mesh.boundary[(7,1)] == 2 # top 350 350 assert mesh.boundary[(7,1)] == 'top' # top 351 assert mesh.boundary[(3,1)] == 'top' # top 352 353 351 assert mesh.boundary[(3,1)] == 'top' # top 352 353 354 354 def test_boundary_tags(self): 355 355 356 356 357 357 points, vertices, boundary = rectangular(4, 4) 358 358 mesh = Mesh(points, vertices, boundary) 359 360 361 #Test that points are arranged in a counter clock wise order 359 360 361 #Test that points are arranged in a counter clock wise order 362 362 mesh.check_integrity() 363 363 … … 369 369 370 370 for k in [24, 26, 28, 30]: 371 assert mesh.boundary[(k,2)] == 'right' 372 371 assert mesh.boundary[(k,2)] == 'right' 372 373 373 for k in [7, 15, 23, 31]: 374 374 assert mesh.boundary[(k,1)] == 'top' … … 377 377 378 378 379 379 380 380 def test_rectangular_mesh(self): 381 381 M=4 … … 383 383 len1 = 100.0 384 384 len2 = 17.0 385 385 386 386 points, vertices, boundary = rectangular(M, N, len1, len2) 387 387 mesh = Mesh(points, vertices, boundary) 388 388 389 assert len(mesh) == 2*M*N 389 assert len(mesh) == 2*M*N 390 390 391 391 for i in range(len(mesh)): … … 396 396 assert mesh.edgelengths[i, 1] == len1/M #x direction 397 397 assert mesh.edgelengths[i, 2] == len2/N #y direction 398 399 #Test that points are arranged in a counter clock wise order 400 mesh.check_integrity() 398 399 #Test that points are arranged in a counter clock wise order 400 mesh.check_integrity() 401 401 402 402 … … 404 404 #Check that integers don't cause trouble 405 405 N = 16 406 406 407 407 points, vertices, boundary = rectangular(2*N, N, len1=10, len2=10) 408 408 mesh = Mesh(points, vertices, boundary) … … 422 422 #bac, bce, ecf, dbe 423 423 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 424 mesh = Mesh(points, vertices) 424 mesh = Mesh(points, vertices) 425 425 mesh.check_integrity() 426 426 … … 428 428 T = mesh 429 429 tid = 0 430 assert T.number_of_boundaries[tid] == 2 431 assert T.surrogate_neighbours[tid, 0] == tid 432 assert T.surrogate_neighbours[tid, 1] == 1 433 assert T.surrogate_neighbours[tid, 2] == tid 430 assert T.number_of_boundaries[tid] == 2 431 assert T.surrogate_neighbours[tid, 0] == tid 432 assert T.surrogate_neighbours[tid, 1] == 1 433 assert T.surrogate_neighbours[tid, 2] == tid 434 434 435 435 tid = 1 436 assert T.number_of_boundaries[tid] == 0 436 assert T.number_of_boundaries[tid] == 0 437 437 assert T.surrogate_neighbours[tid, 0] == 2 438 438 assert T.surrogate_neighbours[tid, 1] == 3 … … 440 440 441 441 tid = 2 442 assert T.number_of_boundaries[tid] == 2 443 assert T.surrogate_neighbours[tid, 0] == tid 444 assert T.surrogate_neighbours[tid, 1] == tid 442 assert T.number_of_boundaries[tid] == 2 443 assert T.surrogate_neighbours[tid, 0] == tid 444 assert T.surrogate_neighbours[tid, 1] == tid 445 445 assert T.surrogate_neighbours[tid, 2] == 1 446 446 447 447 tid = 3 448 assert T.number_of_boundaries[tid] == 2 449 assert T.surrogate_neighbours[tid, 0] == 1 448 assert T.number_of_boundaries[tid] == 2 449 assert T.surrogate_neighbours[tid, 0] == 1 450 450 assert T.surrogate_neighbours[tid, 1] == tid 451 451 assert T.surrogate_neighbours[tid, 2] == tid … … 470 470 (2, 1): 'Fourth', 471 471 (3, 1): 'Fifth', 472 (3, 2): 'Sixth'} 473 474 475 mesh = Mesh(points, vertices, boundary) 472 (3, 2): 'Sixth'} 473 474 475 mesh = Mesh(points, vertices, boundary) 476 476 mesh.check_integrity() 477 477 478 478 479 479 #Check enumeration 480 #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments): 480 #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments): 481 481 # b = -k-1 482 482 # assert mesh.neighbours[vol_id, edge_id] == b … … 502 502 (2, 1): 'Fourth', 503 503 #(3, 1): 'Fifth', #Skip this 504 (3, 2): 'Sixth'} 505 506 507 mesh = Mesh(points, vertices, boundary) 504 (3, 2): 'Sixth'} 505 506 507 mesh = Mesh(points, vertices, boundary) 508 508 mesh.check_integrity() 509 509 … … 513 513 514 514 #Check enumeration 515 #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments): 515 #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments): 516 516 # b = -k-1 517 517 # assert mesh.neighbours[vol_id, edge_id] == b … … 535 535 (2, 1): 'Fourth', 536 536 #(3, 1): 'Fifth', #Skip this 537 (3, 2): 'Sixth'} 538 539 540 mesh = Mesh(points, vertices) #, boundary) 537 (3, 2): 'Sixth'} 538 539 540 mesh = Mesh(points, vertices) #, boundary) 541 541 mesh.check_integrity() 542 542 … … 545 545 assert mesh.boundary[ (0, 2) ] == default_boundary_tag 546 546 assert mesh.boundary[ (2, 0) ] == default_boundary_tag 547 assert mesh.boundary[ (2, 1) ] == default_boundary_tag 547 assert mesh.boundary[ (2, 1) ] == default_boundary_tag 548 548 assert mesh.boundary[ (3, 1) ] == default_boundary_tag 549 assert mesh.boundary[ (3, 2) ] == default_boundary_tag 549 assert mesh.boundary[ (3, 2) ] == default_boundary_tag 550 550 551 551 552 552 #Check enumeration 553 #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments): 553 #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments): 554 554 # b = -k-1 555 555 # assert mesh.neighbours[vol_id, edge_id] == b … … 575 575 #Too few points 576 576 try: 577 mesh = Mesh([points[0]], vertices) 577 mesh = Mesh([points[0]], vertices) 578 578 except AssertionError: 579 579 pass … … 583 583 #Too few points - 1 element 584 584 try: 585 mesh = Mesh([points[0]], [vertices[0]]) 586 except AssertionError: 587 pass 588 else: 589 raise 'Should have raised an exception' 590 591 #Wrong dimension of vertices 592 try: 593 mesh = Mesh(points, vertices[0]) 585 mesh = Mesh([points[0]], [vertices[0]]) 594 586 except AssertionError: 595 587 pass … … 597 589 raise 'Should have raised an exception' 598 590 591 #Wrong dimension of vertices 592 try: 593 mesh = Mesh(points, vertices[0]) 594 except AssertionError: 595 pass 596 else: 597 raise 'Should have raised an exception' 598 599 599 #Unsubscriptable coordinates object raises exception 600 600 try: 601 mesh = Mesh(points[0], [vertices[0]]) 601 mesh = Mesh(points[0], [vertices[0]]) 602 602 except AssertionError: 603 603 pass … … 610 610 #Not specifying all boundary tags 611 611 #try: 612 # mesh = Mesh(points, vertices, {(3,0): 'x'}) 612 # mesh = Mesh(points, vertices, {(3,0): 'x'}) 613 613 #except AssertionError: 614 614 # pass … … 618 618 #Specifying wrong non existing segment 619 619 try: 620 mesh = Mesh(points, vertices, {(5,0): 'x'}) 620 mesh = Mesh(points, vertices, {(5,0): 'x'}) 621 621 except AssertionError: 622 622 pass … … 634 634 from shallow_water import Domain 635 635 from Numeric import zeros, Float 636 636 637 637 #Create basic mesh 638 638 points, vertices, boundary = rectangular(1, 3) … … 648 648 'all':[0,1,2,3,4,5]}) 649 649 650 650 651 651 def test_boundary_polygon(self): 652 652 from mesh_factory import rectangular … … 654 654 from Numeric import zeros, Float 655 655 from util import inside_polygon 656 656 657 657 #Create basic mesh 658 658 points, vertices, boundary = rectangular(2, 2) 659 659 mesh = Mesh(points, vertices, boundary) 660 660 661 661 662 662 P = mesh.get_boundary_polygon() … … 668 668 for p in points: 669 669 #print p, P 670 assert inside_polygon(p, P) 670 assert inside_polygon(p, P) 671 671 672 672 … … 675 675 from Numeric import zeros, Float 676 676 from util import inside_polygon 677 677 678 678 #Points 679 679 a = [0.0, 0.0] #0 … … 695 695 696 696 mesh = Mesh(points, vertices) 697 698 mesh.check_integrity() 699 697 698 mesh.check_integrity() 699 700 700 P = mesh.get_boundary_polygon() 701 701 … … 705 705 for p in points: 706 706 #print p, P 707 assert inside_polygon(p, P) 707 assert inside_polygon(p, P) 708 708 709 709 … … 711 711 """Same as II but vertices ordered differently 712 712 """ 713 713 714 714 from mesh import Mesh 715 715 from Numeric import zeros, Float 716 716 from util import inside_polygon 717 717 718 718 #Points 719 719 a = [0.0, 0.0] #0 … … 738 738 mesh.check_integrity() 739 739 740 740 741 741 P = mesh.get_boundary_polygon() 742 742 … … 745 745 746 746 for p in points: 747 assert inside_polygon(p, P) 748 749 750 747 assert inside_polygon(p, P) 748 749 750 751 751 #------------------------------------------------------------- 752 752 if __name__ == "__main__": 753 suite = unittest.makeSuite(Test Case,'test')753 suite = unittest.makeSuite(Test_Mesh,'test') 754 754 runner = unittest.TextTestRunner() 755 755 runner.run(suite) 756 756 757 758 759 760 761 757 758 759 760 -
inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py
r969 r1018 23 23 Transmissive_boundary 24 24 25 26 class Test Case(unittest.TestCase):25 26 class Test_pmesh2domain(unittest.TestCase): 27 27 28 28 def setUp(self): 29 29 pass 30 30 31 31 def tearDown(self): 32 32 pass 33 33 34 34 def test_pmesh2Domain(self): 35 35 import os 36 36 import tempfile 37 37 38 38 fileName = tempfile.mktemp(".tsh") 39 39 file = open(fileName,"w") … … 66 66 0 # <x> <y> <attribute>...Mesh Regions... \n") 67 67 file.close() 68 68 69 69 tags = {} 70 70 b1 = Dirichlet_boundary(conserved_quantities = array([0.0])) … … 74 74 tags["2"] = b2 75 75 tags["3"] = b3 76 76 77 77 domain = pmesh_to_domain_instance(fileName, Domain) 78 78 os.remove(fileName) 79 79 80 ## check the quantities 80 ## check the quantities 81 81 #print domain.quantities['elevation'].vertex_values 82 82 answer = [[0., 8., 0.], 83 83 [0., 10., 8.]] 84 84 assert allclose(domain.quantities['elevation'].vertex_values, 85 answer) 86 85 answer) 86 87 87 #print domain.quantities['stage'].vertex_values 88 88 answer = [[0., 12., 10.], … … 90 90 assert allclose(domain.quantities['stage'].vertex_values, 91 91 answer) 92 92 93 93 #print domain.quantities['friction'].vertex_values 94 94 answer = [[0.01, 0.04, 0.03], … … 100 100 assert allclose(domain.tagged_elements['dsg'][0],0) 101 101 assert allclose(domain.tagged_elements['ole nielsen'][0],1) 102 102 103 103 self.failUnless( domain.boundary[(1, 0)] == '1', 104 104 "test_tags_to_boundaries failed. Single boundary wasn't added.") … … 123 123 meshDict['triangle_tags'] = [6.6,6.6] 124 124 meshDict['triangle_neighbors'] = [[-1,-1,1],[-1,0,-1]] 125 meshDict['segments'] = [[1,0],[0,2],[2,3],[3,1]] 125 meshDict['segments'] = [[1,0],[0,2],[2,3],[3,1]] 126 126 meshDict['segment_tags'] = [2,3,1,1] 127 127 128 128 domain = Domain.pmesh_dictionary_to_domain(meshDict) 129 129 … … 143 143 boundary_value = Boundary_value.instances[id] 144 144 boundary_obj = boundary_value.boundary_object 145 145 146 146 self.failUnless( boundary_obj == b1, 147 147 "test_tags_to_boundaries failed. Single boundary wasn't added.") 148 148 149 149 inverted_id = Volume.instances[0].neighbours[1] 150 150 id = -(inverted_id+1) 151 151 boundary_value = Boundary_value.instances[id] 152 boundary_obj = boundary_value.boundary_object 152 boundary_obj = boundary_value.boundary_object 153 153 self.failUnless( boundary_obj == b3, 154 154 "test_tags_to_boundaries failed. Single boundary wasn't added.") … … 156 156 #------------------------------------------------------------- 157 157 if __name__ == "__main__": 158 suite = unittest.makeSuite(Test Case, 'test')158 suite = unittest.makeSuite(Test_pmesh2domain, 'test') 159 159 runner = unittest.TextTestRunner() 160 160 runner.run(suite) 161 162 163 164 161 165 162 166 163 164 165 -
inundation/ga/storm_surge/pyvolution/test_quad.py
r707 r1018 9 9 #------------------------------------------------------------- 10 10 11 class Test Case(unittest.TestCase):11 class Test_Quad(unittest.TestCase): 12 12 13 13 def setUp(self): 14 14 self.cell = Cell(100, 140, 0, 40, 'cell') 15 15 16 16 a = [3, 107] 17 17 b = [5, 107] … … 24 24 25 25 points = [a, b, c, d, e, f, g, h] 26 #bac, bce, ecf, dbe, daf, dae 27 vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]] 26 #bac, bce, ecf, dbe, daf, dae 27 vertices = [[1,0,2], [1,3,4], [1,2,3], [5,4,7], [4,6,7]] 28 28 29 29 mesh = Mesh(points, vertices) 30 30 self.mesh = mesh 31 31 Cell.initialise(mesh) 32 32 33 33 def tearDown(self): 34 34 pass … … 37 37 self.cell.insert(0) 38 38 self.cell.insert(1) 39 39 40 40 result = self.cell.retrieve() 41 41 assert type(result) in [types.ListType,types.TupleType],\ 42 'should be a list' 42 'should be a list' 43 43 self.assertEqual(len(result),2) 44 44 45 45 def test_add_points_2_cellII(self): 46 46 self.cell.insert([0,1,2,3,4,5,6,7]) 47 47 48 48 result = self.cell.retrieve() 49 49 assert type(result) in [types.ListType,types.TupleType],\ 50 'should be a list' 50 'should be a list' 51 51 self.assertEqual(len(result),8) 52 52 … … 55 55 self.cell.insert([0,1,2,3,4,5,6,7]) 56 56 self.cell.split(4) 57 57 58 58 result = self.cell.search(x = 1, y = 101) 59 59 assert type(result) in [types.ListType,types.TupleType],\ 60 'should be a list' 60 'should be a list' 61 61 self.assertEqual(result, [0,1,2,3]) 62 63 62 63 64 64 def test_clear_1(self): 65 65 self.cell.insert([0,1,2,3,4,5,6,7]) 66 assert self.cell.count() == 8 66 assert self.cell.count() == 8 67 67 self.cell.clear() 68 68 69 69 #This one actually revealed a bug :-) 70 assert self.cell.count() == 0 71 70 assert self.cell.count() == 0 71 72 72 def test_clear_2(self): 73 73 self.cell.insert([0,1,2,3,4,5,6,7]) 74 assert self.cell.count() == 8 74 assert self.cell.count() == 8 75 75 self.cell.split(2) 76 assert self.cell.count() == 8 77 76 assert self.cell.count() == 8 77 78 78 self.cell.clear() 79 assert self.cell.count() == 0 80 81 82 79 assert self.cell.count() == 0 80 81 82 83 83 def test_split(self): 84 self.cell.insert([0,1,2,3,4,5,6,7], split = False) 85 84 self.cell.insert([0,1,2,3,4,5,6,7], split = False) 85 86 86 #No children yet 87 assert self.cell.children is None 87 assert self.cell.children is None 88 88 assert self.cell.count() == 8 89 89 … … 92 92 #self.cell.show() 93 93 #self.cell.show_all() 94 95 94 95 96 96 #Now there are children 97 assert self.cell.children is not None 98 assert self.cell.count() == 8 99 100 101 97 assert self.cell.children is not None 98 assert self.cell.count() == 8 99 100 101 102 102 def test_collapse(self): 103 self.cell.insert([0,1,2,3,4,5,6,7], split = False) 104 103 self.cell.insert([0,1,2,3,4,5,6,7], split = False) 104 105 105 #Split maximally 106 106 self.cell.split(1) 107 107 108 108 #Now there are children 109 assert self.cell.children is not None 110 assert self.cell.count() == 8 111 109 assert self.cell.children is not None 110 assert self.cell.count() == 8 111 112 112 #Collapse 113 113 self.cell.collapse(8) 114 114 115 115 #No children 116 assert self.cell.children is None 117 assert self.cell.count() == 8 116 assert self.cell.children is None 117 assert self.cell.count() == 8 118 118 119 119 def test_build_quadtree(self): … … 123 123 124 124 #Q.show() 125 125 126 126 result = Q.search(3, 105) 127 127 assert type(result) in [types.ListType,types.TupleType],\ … … 130 130 self.assertEqual(result, [0,1,2,3]) 131 131 132 132 133 133 def test_build_quadtreeII(self): 134 134 135 135 self.cell = Cell(100, 140, 0, 40, 'cell') 136 136 137 137 p0 = [34.6292076111,-7999.92529297] 138 138 p1 = [8000.0, 7999.0] … … 141 141 142 142 points = [p0,p1,p2, p3] 143 #bac, bce, ecf, dbe, daf, dae 143 #bac, bce, ecf, dbe, daf, dae 144 144 vertices = [[0,1,2],[0,2,3]] 145 145 … … 148 148 #This was causing round off error 149 149 Q = build_quadtree(mesh) 150 150 151 151 #------------------------------------------------------------- 152 152 if __name__ == "__main__": 153 153 154 mysuite = unittest.makeSuite(Test Case,'test')154 mysuite = unittest.makeSuite(Test_Sparse,'test') 155 155 runner = unittest.TextTestRunner() 156 156 runner.run(mysuite) 157 -
inundation/ga/storm_surge/pyvolution/test_quantity.py
r1011 r1018 9 9 from Numeric import allclose, array, ones, Float 10 10 11 12 class Test Case(unittest.TestCase):11 12 class Test_Quantity(unittest.TestCase): 13 13 def setUp(self): 14 14 from domain import Domain 15 15 16 16 a = [0.0, 0.0] 17 17 b = [0.0, 2.0] … … 26 26 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 27 27 28 self.mesh1 = Domain(points[:3], [elements[0]]) 28 self.mesh1 = Domain(points[:3], [elements[0]]) 29 29 self.mesh1.check_integrity() 30 31 self.mesh4 = Domain(points, elements) 30 31 self.mesh4 = Domain(points, elements) 32 32 self.mesh4.check_integrity() 33 33 34 34 def tearDown(self): 35 35 pass … … 38 38 39 39 def test_creation(self): 40 40 41 41 quantity = Quantity(self.mesh1, [[1,2,3]]) 42 42 assert allclose(quantity.vertex_values, [[1.,2.,3.]]) … … 55 55 pass 56 56 except: 57 raise 'Should have raised "mising mesh object" error' 58 57 raise 'Should have raised "mising mesh object" error' 58 59 59 60 60 def test_creation_zeros(self): 61 61 62 62 quantity = Quantity(self.mesh1) 63 63 assert allclose(quantity.vertex_values, [[0.,0.,0.]]) … … 66 66 quantity = Quantity(self.mesh4) 67 67 assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.], 68 [0.,0.,0.], [0.,0.,0.]]) 69 70 68 [0.,0.,0.], [0.,0.,0.]]) 69 70 71 71 def test_interpolation(self): 72 72 quantity = Quantity(self.mesh1, [[1,2,3]]) 73 73 assert allclose(quantity.centroid_values, [2.0]) #Centroid 74 74 75 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 75 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 76 76 77 77 … … 88 88 [7.5, 0.5, 1.], 89 89 [-5, -2.5, 7.5]]) 90 91 #print quantity.edge_values 90 91 #print quantity.edge_values 92 92 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 93 93 [5., 5., 5.], … … 120 120 quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 121 121 assert allclose(quantity.vertex_values, 122 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 122 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 123 123 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 124 124 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], … … 135 135 location = 'edges') 136 136 assert allclose(quantity.edge_values, 137 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 137 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 138 138 139 139 #Test exceptions … … 174 174 [3, 3, 3], 175 175 [3, 3, 3], 176 [3, 3, 3]]) 176 [3, 3, 3]]) 177 177 178 178 … … 186 186 #print "quantity.vertex_values",quantity.vertex_values 187 187 assert allclose(quantity.vertex_values, 188 [[2,0,2], [2,2,4], [4,2,4], [4,2,4]]) 188 [[2,0,2], [2,2,4], [4,2,4], [4,2,4]]) 189 189 assert allclose(quantity.centroid_values, 190 190 [4.0/3, 8.0/3, 10.0/3, 10.0/3]) 191 191 assert allclose(quantity.edge_values, 192 [[1,2,1], [3,3,2], [3,4,3], [3,4,3]]) 193 194 192 [[1,2,1], [3,3,2], [3,4,3], [3,4,3]]) 193 194 195 195 quantity.set_values(f, location = 'centroids') 196 196 assert allclose(quantity.centroid_values, … … 201 201 quantity = Quantity(self.mesh4) 202 202 203 #Try constants first 203 #Try constants first 204 204 const = 5 205 quantity.set_values(const, location = 'vertices') 205 quantity.set_values(const, location = 'vertices') 206 206 #print 'Q', quantity.get_integral() 207 208 assert allclose(quantity.get_integral(), 207 208 assert allclose(quantity.get_integral(), 209 209 self.mesh4.get_area() * const) 210 210 211 211 #Try with a linear function 212 212 def f(x, y): … … 214 214 215 215 quantity.set_values(f, location = 'vertices') 216 217 216 217 218 218 ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 219 220 assert allclose (quantity.get_integral(), ref_integral) 219 220 assert allclose (quantity.get_integral(), ref_integral) 221 221 222 222 … … 256 256 #Centroid 257 257 assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) 258 258 259 259 assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 260 260 [3., 2.5, 1.5], … … 268 268 quantity = Conserved_quantity(self.mesh4) 269 269 270 #Set up for a gradient of (3,0) at mid triangle 270 #Set up for a gradient of (3,0) at mid triangle 271 271 quantity.set_values([2.0, 4.0, 8.0, 2.0], 272 272 location = 'centroids') … … 294 294 [0., 0., 6.]]) 295 295 296 296 297 297 def test_second_order_extrapolation2(self): 298 quantity = Conserved_quantity(self.mesh4) 298 quantity = Conserved_quantity(self.mesh4) 299 299 300 300 #Set up for a gradient of (3,1), f(x) = 3x+y 301 301 quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3], 302 302 location = 'centroids') 303 303 304 304 #Gradients 305 305 a, b = quantity.compute_gradients() … … 314 314 315 315 quantity.extrapolate_second_order() 316 316 317 317 #print quantity.vertex_values 318 318 assert allclose(quantity.vertex_values[1,0], 2.0) 319 assert allclose(quantity.vertex_values[1,1], 6.0) 319 assert allclose(quantity.vertex_values[1,1], 6.0) 320 320 assert allclose(quantity.vertex_values[1,2], 8.0) 321 321 … … 323 323 324 324 def test_first_order_extrapolator(self): 325 quantity = Conserved_quantity(self.mesh4) 325 quantity = Conserved_quantity(self.mesh4) 326 326 327 327 #Test centroids … … 338 338 339 339 def test_second_order_extrapolator(self): 340 quantity = Conserved_quantity(self.mesh4) 341 342 #Set up for a gradient of (3,0) at mid triangle 340 quantity = Conserved_quantity(self.mesh4) 341 342 #Set up for a gradient of (3,0) at mid triangle 343 343 quantity.set_values([2.0, 4.0, 8.0, 2.0], 344 344 location = 'centroids') 345 345 346 346 347 347 … … 353 353 assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 354 354 assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1] 355 355 356 356 assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 357 357 assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 358 358 359 359 assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 360 360 assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1] 361 361 362 362 363 363 #Assert that quantities are conserved 364 364 from Numeric import sum … … 366 366 assert allclose (quantity.centroid_values[k], 367 367 sum(quantity.vertex_values[k,:])/3) 368 369 370 371 368 369 370 371 372 372 373 373 def test_limiter(self): … … 384 384 assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 385 385 assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1] 386 386 387 387 assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 388 388 assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 389 389 390 390 assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 391 391 assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1] 392 392 393 393 394 394 395 395 #Assert that quantities are conserved 396 396 from Numeric import sum … … 398 398 assert allclose (quantity.centroid_values[k], 399 399 sum(quantity.vertex_values[k,:])/3) 400 400 401 401 402 402 def test_limiter2(self): … … 418 418 quantity.limit() 419 419 420 420 421 421 assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9]) 422 422 423 423 424 424 #Assert that quantities are conserved … … 429 429 430 430 431 431 432 432 433 433 434 434 def test_distribute_first_order(self): 435 quantity = Conserved_quantity(self.mesh4) 435 quantity = Conserved_quantity(self.mesh4) 436 436 437 437 #Test centroids … … 444 444 445 445 #Interpolate 446 quantity.interpolate_from_vertices_to_edges() 446 quantity.interpolate_from_vertices_to_edges() 447 447 448 448 assert allclose(quantity.vertex_values, … … 453 453 454 454 def test_distribute_second_order(self): 455 quantity = Conserved_quantity(self.mesh4) 455 quantity = Conserved_quantity(self.mesh4) 456 456 457 457 #Test centroids … … 498 498 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 499 499 denom = ones(4, Float)-timestep*sem 500 500 501 501 x = array([1, 2, 3, 4])/denom 502 502 assert allclose( quantity.centroid_values, x) … … 512 512 #Set explicit_update 513 513 quantity.explicit_update = array( [4.,3.,2.,1.] ) 514 514 515 515 #Set semi implicit update 516 516 quantity.semi_implicit_update = array( [1.,1.,1.,1.] ) 517 517 518 518 #Update with given timestep 519 timestep = 0.1 519 timestep = 0.1 520 520 quantity.update(0.1) 521 521 … … 525 525 x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) 526 526 x /= denom 527 assert allclose( quantity.centroid_values, x) 528 529 530 531 532 #Test smoothing 527 assert allclose( quantity.centroid_values, x) 528 529 530 531 532 #Test smoothing 533 533 def test_smoothing(self): 534 534 … … 546 546 domain.reduction = mean 547 547 548 548 549 549 #Set some field values 550 domain.set_quantity('elevation', lambda x,y: x) 550 domain.set_quantity('elevation', lambda x,y: x) 551 551 domain.set_quantity('friction', 0.03) 552 552 553 553 554 ###################### 554 ###################### 555 555 # Boundary conditions 556 556 B = Transmissive_boundary(domain) 557 557 domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B}) 558 559 560 ###################### 558 559 560 ###################### 561 561 #Initial condition - with jumps 562 562 … … 566 566 h = 0.03 567 567 for i in range(stage.shape[0]): 568 if i % 2 == 0: 568 if i % 2 == 0: 569 569 stage[i,:] = bed[i,:] + h 570 570 else: 571 571 stage[i,:] = bed[i,:] 572 572 573 573 domain.set_quantity('stage', stage) 574 574 … … 579 579 Q = stage.vertex_values 580 580 581 582 assert A.shape[0] == 9 583 assert V.shape[0] == 8 584 assert V.shape[1] == 3 585 586 #First four points 587 assert allclose(A[0], (Q[0,2] + Q[1,1])/2) 581 582 assert A.shape[0] == 9 583 assert V.shape[0] == 8 584 assert V.shape[1] == 3 585 586 #First four points 587 assert allclose(A[0], (Q[0,2] + Q[1,1])/2) 588 588 assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3) 589 589 assert allclose(A[2], Q[3,0]) … … 593 593 assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\ 594 594 Q[5,0] + Q[6,2] + Q[7,1])/6) 595 595 596 596 597 597 #Check V 598 598 assert allclose(V[0,:], [3,4,0]) 599 599 assert allclose(V[1,:], [1,0,4]) 600 assert allclose(V[2,:], [4,5,1]) 600 assert allclose(V[2,:], [4,5,1]) 601 601 assert allclose(V[3,:], [2,1,5]) 602 assert allclose(V[4,:], [6,7,3]) 602 assert allclose(V[4,:], [6,7,3]) 603 603 assert allclose(V[5,:], [4,3,7]) 604 604 assert allclose(V[6,:], [7,8,4]) 605 assert allclose(V[7,:], [5,4,8]) 605 assert allclose(V[7,:], [5,4,8]) 606 606 607 607 #Get smoothed stage with XY … … 609 609 610 610 assert allclose(A, A1) 611 assert allclose(V, V1) 611 assert allclose(V, V1) 612 612 613 613 #Check XY … … 616 616 617 617 assert allclose(X[7], 1.0) 618 assert allclose(Y[7], 0.5) 618 assert allclose(Y[7], 0.5) 619 619 620 620 … … 628 628 from util import mean 629 629 630 630 631 631 #Create basic mesh 632 632 points, vertices, boundary = rectangular(2, 2) … … 637 637 domain.reduction = mean 638 638 639 639 640 640 #Set some field values 641 domain.set_quantity('elevation', lambda x,y: x) 641 domain.set_quantity('elevation', lambda x,y: x) 642 642 domain.set_quantity('friction', 0.03) 643 643 644 644 645 ###################### 645 ###################### 646 646 #Initial condition - with jumps 647 647 … … 651 651 h = 0.03 652 652 for i in range(stage.shape[0]): 653 if i % 2 == 0: 653 if i % 2 == 0: 654 654 stage[i,:] = bed[i,:] + h 655 655 else: 656 656 stage[i,:] = bed[i,:] 657 657 658 658 domain.set_quantity('stage', stage) 659 659 660 660 #Get stage 661 stage = domain.quantities['stage'] 661 stage = domain.quantities['stage'] 662 662 A, V = stage.get_vertex_values(xy=False, smooth=False) 663 663 Q = stage.vertex_values.flat … … 666 666 assert allclose(A[k], Q[k]) 667 667 668 668 669 669 for k in range(8): 670 670 assert V[k, 0] == 3*k 671 671 assert V[k, 1] == 3*k+1 672 assert V[k, 2] == 3*k+2 673 672 assert V[k, 2] == 3*k+2 673 674 674 675 675 676 676 X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False) 677 677 678 678 679 679 assert allclose(A, A1) 680 assert allclose(V, V1) 680 assert allclose(V, V1) 681 681 682 682 #Check XY … … 686 686 assert allclose(Y[4], 0.0) 687 687 assert allclose(X[12], 1.0) 688 assert allclose(Y[12], 0.0) 689 690 691 688 assert allclose(Y[12], 0.0) 689 690 691 692 692 def set_array_values_by_index(self): 693 693 … … 695 695 from shallow_water import Domain 696 696 from Numeric import zeros, Float 697 697 698 698 #Create basic mesh 699 699 points, vertices, boundary = rectangular(1, 1) … … 701 701 #Create shallow water domain 702 702 domain = Domain(points, vertices, boundary) 703 #print "domain.number_of_elements ",domain.number_of_elements 703 #print "domain.number_of_elements ",domain.number_of_elements 704 704 quantity = Quantity(domain,[[1,1,1],[2,2,2]]) 705 705 value = [7] … … 709 709 indexes = indexes) 710 710 #print "quantity.centroid_values",quantity.centroid_values 711 711 712 712 assert allclose(quantity.centroid_values, [1,7]) 713 713 … … 717 717 quantity.set_array_values([15,20,25], indexes = indexes) 718 718 assert allclose(quantity.centroid_values, [1,20]) 719 719 720 720 def test_setting_some_vertex_values(self): 721 721 """ … … 725 725 from shallow_water import Domain 726 726 from Numeric import zeros, Float 727 727 728 728 #Create basic mesh 729 729 points, vertices, boundary = rectangular(1, 3) … … 731 731 #Create shallow water domain 732 732 domain = Domain(points, vertices, boundary) 733 #print "domain.number_of_elements ",domain.number_of_elements 733 #print "domain.number_of_elements ",domain.number_of_elements 734 734 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 735 735 [4,4,4],[5,5,5],[6,6,6]]) … … 741 741 #print "quantity.centroid_values",quantity.centroid_values 742 742 assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 743 743 744 744 value = [[15,20,25]] 745 745 quantity.set_values(value, indexes = indexes) … … 758 758 assert allclose(quantity.centroid_values, [10,100,3,4,5,50]) 759 759 760 760 761 761 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 762 762 [4,4,4],[5,5,5],[6,6,6]]) … … 764 764 #this will be per unique vertex, indexing the vertices 765 765 #print "quantity.vertex_values",quantity.vertex_values 766 quantity.set_values(values, indexes = [0,1,5]) 766 quantity.set_values(values, indexes = [0,1,5]) 767 767 #print "quantity.vertex_values",quantity.vertex_values 768 768 assert allclose(quantity.vertex_values[0], [1,50,10]) … … 780 780 [4,4,4],[5,5,5],[6,6,6]] 781 781 quantity.set_values(values) 782 782 783 783 # testing the standard set values by vertex 784 784 # indexed by vertex_id in general_mesh.coordinates 785 785 values = [0,1,2,3,4,5,6,7] 786 786 787 787 quantity.set_values(values) 788 788 #print "1 quantity.vertex_values",quantity.vertex_values … … 793 793 [ 6., 7., 2.], 794 794 [ 3., 2., 7.]]) 795 795 796 796 def test_setting_unique_vertex_values(self): 797 797 """ … … 801 801 from shallow_water import Domain 802 802 from Numeric import zeros, Float 803 803 804 804 #Create basic mesh 805 805 points, vertices, boundary = rectangular(1, 3) … … 807 807 #Create shallow water domain 808 808 domain = Domain(points, vertices, boundary) 809 #print "domain.number_of_elements ",domain.number_of_elements 809 #print "domain.number_of_elements ",domain.number_of_elements 810 810 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 811 811 [4,4,4],[5,5,5]]) … … 819 819 assert allclose(quantity.vertex_values[1], [7,1,7]) 820 820 assert allclose(quantity.vertex_values[2], [7,2,7]) 821 822 821 822 823 823 def test_get_values(self): 824 824 """ … … 828 828 from shallow_water import Domain 829 829 from Numeric import zeros, Float 830 830 831 831 #Create basic mesh 832 832 points, vertices, boundary = rectangular(1, 3) … … 838 838 #Create shallow water domain 839 839 domain = Domain(points, vertices, boundary) 840 #print "domain.number_of_elements ",domain.number_of_elements 840 #print "domain.number_of_elements ",domain.number_of_elements 841 841 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 842 842 [4,4,4],[5,5,5]]) 843 843 844 844 #print "quantity.get_values(location = 'unique vertices')", \ 845 845 # quantity.get_values(location = 'unique vertices') … … 848 848 # quantity.get_values(indexes=[0,1,2,3,4,5,6,7], \ 849 849 # location = 'unique vertices') 850 850 851 851 answer = [0.5,2,4,5,0,1,3,4.5] 852 852 assert allclose(answer, 853 853 quantity.get_values(location = 'unique vertices')) 854 854 855 855 indexes = [0,5,3] 856 856 answer = [0.5,1,5] … … 861 861 #print "quantity.get_values(location = 'centroids') ",\ 862 862 # quantity.get_values(location = 'centroids') 863 863 864 864 def test_getting_some_vertex_values(self): 865 865 """ … … 869 869 from shallow_water import Domain 870 870 from Numeric import zeros, Float 871 871 872 872 #Create basic mesh 873 873 points, vertices, boundary = rectangular(1, 3) … … 879 879 #Create shallow water domain 880 880 domain = Domain(points, vertices, boundary) 881 #print "domain.number_of_elements ",domain.number_of_elements 881 #print "domain.number_of_elements ",domain.number_of_elements 882 882 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 883 883 [4,4,4],[5,5,5],[6,6,6]]) … … 893 893 quantity.get_values(location = 'centroids')) 894 894 895 895 896 896 value = [[15,20,25]] 897 897 quantity.set_values(value, indexes = indexes) … … 902 902 quantity.get_values(location = 'edges')) 903 903 904 # get a subset of elements 904 # get a subset of elements 905 905 subset = quantity.get_values(location='centroids', indexes=[0,5]) 906 906 answer = [quantity.centroid_values[0],quantity.centroid_values[5]] … … 920 920 assert allclose(subset, answer) 921 921 922 923 922 923 924 924 #------------------------------------------------------------- 925 925 if __name__ == "__main__": 926 suite = unittest.makeSuite(Test Case,'test')926 suite = unittest.makeSuite(Test_Quantity,'test') 927 927 #print "restricted test" 928 #suite = unittest.makeSuite(Test Case,'test_set_vertex_values_subset')928 #suite = unittest.makeSuite(Test_Quantity,'test_set_vertex_values_subset') 929 929 runner = unittest.TextTestRunner() 930 930 runner.run(suite) 931 931 932 -
inundation/ga/storm_surge/pyvolution/test_region.py
r773 r1018 12 12 def add_x_y(x, y): 13 13 return x+y 14 15 class Test Case(unittest.TestCase):14 15 class Test_Region(unittest.TestCase): 16 16 def setUp(self): 17 17 pass 18 18 19 19 20 20 def tearDown(self): 21 21 pass … … 29 29 from shallow_water import Domain 30 30 from Numeric import zeros, Float 31 31 32 32 #Create basic mesh 33 33 points, vertices, boundary = rectangular(1, 3) 34 34 35 35 #Create shallow water domain 36 36 domain = Domain(points, vertices, boundary) … … 39 39 'all':[0,1,2,3,4,5]}) 40 40 41 41 42 42 #Set friction 43 43 manning = 0.07 … … 66 66 [ 11.0, 11.0, 11.0], 67 67 [ 11.0, 11.0, 11.0]]) 68 69 # trying a function 68 69 # trying a function 70 70 domain.set_region(Set_region('top', 'friction', add_x_y)) 71 71 #print domain.quantities['friction'].get_values() … … 77 77 [ 5./3, 2.0, 2./3], 78 78 [ 1.0, 2./3, 2.0]]) 79 79 80 80 domain.set_quantity('elevation', 10.0) 81 81 domain.set_quantity('stage', 10.0) … … 89 89 [ 11.0, 11.0, 11.0], 90 90 [ 11.0, 11.0, 11.0]]) 91 91 92 92 def test_unique_vertices(self): 93 93 """ … … 97 97 from shallow_water import Domain 98 98 from Numeric import zeros, Float 99 99 100 100 #Create basic mesh 101 101 points, vertices, boundary = rectangular(1, 3) 102 102 103 103 #Create shallow water domain 104 104 domain = Domain(points, vertices, boundary) … … 106 106 'top':[4,5], 107 107 'all':[0,1,2,3,4,5]}) 108 108 109 109 #Set friction 110 110 manning = 0.07 … … 122 122 [ 0.07, 0.07, 0.07]]) 123 123 124 124 125 125 def test_unique_verticesII(self): 126 126 """ … … 130 130 from shallow_water import Domain 131 131 from Numeric import zeros, Float 132 132 133 133 #Create basic mesh 134 134 points, vertices, boundary = rectangular(1, 3) 135 135 136 136 #Create shallow water domain 137 137 domain = Domain(points, vertices, boundary) … … 139 139 'top':[4,5], 140 140 'all':[0,1,2,3,4,5]}) 141 141 142 142 #Set friction 143 143 manning = 0.07 … … 145 145 146 146 domain.set_region(Add_value_to_region('bottom', 'friction', 1.0,initial_quantity='friction', location = 'unique vertices')) 147 147 148 148 #print domain.quantities['friction'].get_values() 149 149 assert allclose(domain.quantities['friction'].get_values(),\ … … 156 156 #------------------------------------------------------------- 157 157 if __name__ == "__main__": 158 suite = unittest.makeSuite(Test Case,'test')158 suite = unittest.makeSuite(Test_Region,'test') 159 159 runner = unittest.TextTestRunner() 160 160 runner.run(suite) -
inundation/ga/storm_surge/pyvolution/test_sparse.py
r599 r1018 5 5 6 6 from sparse import * 7 from Numeric import allclose, array, transpose, Float 8 9 class Test Case(unittest.TestCase):7 from Numeric import allclose, array, transpose, Float 8 9 class Test_Sparse(unittest.TestCase): 10 10 11 11 def setUp(self): 12 12 pass 13 13 14 14 def tearDown(self): 15 15 pass … … 17 17 def test_init1(Self): 18 18 """Test initialisation from dimensions 19 """ 19 """ 20 20 A = Sparse(3,3) 21 21 A[1,1] = 4 … … 26 26 assert A[i,j] == 4.0 27 27 else: 28 assert A[i,j] == 0.0 29 28 assert A[i,j] == 0.0 29 30 30 def test_init2(Self): 31 31 """Test initialisation from dense matrix 32 32 """ 33 33 34 34 A = Sparse(4,3) 35 35 A[1,1] = 4 … … 37 37 A[0,1] = 13 38 38 A[0,2] = -6 39 A[2,2] = 1 39 A[2,2] = 1 40 40 41 41 B = A.todense() 42 42 43 43 C = Sparse(B) 44 44 45 45 assert allclose(C.todense(), B) 46 46 47 47 48 48 def test_dense(self): 49 49 A = Sparse(4,3) 50 A[1,1] = 4 51 50 A[1,1] = 4 51 52 52 assert allclose(A.todense(), [[0,0,0], [0,4,0], [0,0,0], [0,0,0]]) 53 53 … … 59 59 60 60 A = Sparse(3,3) 61 A[1,1] = 4 61 A[1,1] = 4 62 62 A[1,1] = 0 63 63 … … 68 68 A[1,2] = 0 69 69 assert len(A) == 0 70 assert allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]]) 70 assert allclose(A.todense(), [[0,0,0], [0,0,0], [0,0,0]]) 71 71 72 72 def test_sparse_multiplication_vector(self): 73 73 A = Sparse(3,3) 74 74 75 75 A[0,0] = 3 76 76 A[1,1] = 2 … … 91 91 92 92 u = A*v[:,1] 93 assert allclose(u, [12,16,4]) 93 assert allclose(u, [12,16,4]) 94 94 95 95 96 96 def test_sparse_multiplication_matrix(self): 97 97 A = Sparse(3,3) 98 98 99 99 A[0,0] = 3 100 100 A[1,1] = 2 … … 106 106 107 107 u = A*v 108 assert allclose(u, [[6,12], [14,16], [4,4]]) 109 110 111 112 def test_sparse_transpose_multiplication(self): 113 A = Sparse(3,3) 114 108 assert allclose(u, [[6,12], [14,16], [4,4]]) 109 110 111 112 def test_sparse_transpose_multiplication(self): 113 A = Sparse(3,3) 114 115 115 A[0,0] = 3 116 116 A[1,1] = 2 … … 128 128 """Test method __rmul__ 129 129 """ 130 130 131 131 A = Sparse(3,3) 132 132 … … 140 140 141 141 B = A*3 142 assert allclose(B.todense(), 3*A.todense()) 142 assert allclose(B.todense(), 3*A.todense()) 143 143 144 144 try: … … 152 152 """ Test sparse addition with dok format 153 153 """ 154 154 155 155 A = Sparse(3,3) 156 156 … … 171 171 """ Test conversion to csr format 172 172 """ 173 173 174 174 A = Sparse(4,3) 175 175 … … 181 181 A[2,0] = 5 182 182 183 #print ' ' 183 #print ' ' 184 184 #print A.todense() 185 185 … … 197 197 198 198 assert allclose(B*C2, [[15.0, 30.0],[10.0, 20.0],[8.0, 16.0],[0.0, 0.0]]) 199 199 200 200 201 201 202 202 #------------------------------------------------------------- 203 203 if __name__ == "__main__": 204 suite = unittest.makeSuite(Test Case, 'test')204 suite = unittest.makeSuite(Test_Sparse, 'test') 205 205 runner = unittest.TextTestRunner() 206 206 runner.run(suite) 207 -
inundation/ga/storm_surge/pyvolution/test_util.py
r999 r1018 13 13 return x+y 14 14 15 class Test Case(unittest.TestCase):15 class Test_Util(unittest.TestCase): 16 16 def setUp(self): 17 17 pass … … 25 25 x1 = 1.0; y1 = 0.0; z1 = -1.0 26 26 x2 = 0.0; y2 = 1.0; z2 = 0.0 27 27 28 28 zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) 29 29 … … 44 44 q1 = 8.0+2.0/3 45 45 q2 = 2.0+8.0/3 46 47 #Gradient of fitted pwl surface 46 47 #Gradient of fitted pwl surface 48 48 a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) 49 49 50 50 assert abs(a - 3.0) < epsilon 51 51 assert abs(b - 1.0) < epsilon 52 52 53 53 54 54 … … 67 67 import util_ext 68 68 69 69 70 70 def test_gradient_C_extension(self): 71 71 from util_ext import gradient as gradient_c … … 78 78 q1 = 8.0+2.0/3 79 79 q2 = 2.0+8.0/3 80 81 #Gradient of fitted pwl surface 80 81 #Gradient of fitted pwl surface 82 82 a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2) 83 83 84 84 assert abs(a - 3.0) < epsilon 85 85 assert abs(b - 1.0) < epsilon 86 86 87 87 88 88 def test_gradient_C_extension3(self): 89 89 from util_ext import gradient as gradient_c 90 90 91 91 from RandomArray import uniform, seed 92 92 seed(17, 53) … … 96 96 q0 = uniform(0.0, 10.0, 4) 97 97 q1 = uniform(1.0, 3.0, 4) 98 q2 = uniform(7.0, 20.0, 4) 99 98 q2 = uniform(7.0, 20.0, 4) 99 100 100 101 101 for i in range(4): 102 102 #Gradient of fitted pwl surface 103 from util import gradient_python 103 from util import gradient_python 104 104 a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, 105 105 q0[i], q1[i], q2[i]) … … 112 112 assert abs(a - a_ref) < epsilon 113 113 assert abs(b - b_ref) < epsilon 114 115 116 117 #Geometric 114 115 116 117 #Geometric 118 118 #def test_distance(self): 119 119 # from util import distance# … … 127 127 # 128 128 # self.failUnless( distance([9,8],[4,2]) == distance([4,2],[9,8]), 129 # 'distance is wrong!') 129 # 'distance is wrong!') 130 130 131 131 def test_angle(self): … … 136 136 137 137 def test_anglediff(self): 138 from util import anglediff 138 from util import anglediff 139 139 140 140 assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0) … … 163 163 fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600))) 164 164 t += dt 165 165 166 166 fid.close() 167 167 … … 174 174 175 175 #Exact linear intpolation 176 assert allclose(q[0], 2*t) 176 assert allclose(q[0], 2*t) 177 177 if i%6 == 0: 178 178 assert allclose(q[1], t**2) 179 assert allclose(q[2], sin(t*pi/600)) 179 assert allclose(q[2], sin(t*pi/600)) 180 180 181 181 #Check non-exact 182 182 183 183 t = 90 #Halfway between 60 and 120 184 q = F(t) 184 q = F(t) 185 185 assert allclose( (120**2 + 60**2)/2, q[1] ) 186 186 assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] ) … … 188 188 189 189 t = 100 #Two thirds of the way between between 60 and 120 190 q = F(t) 190 q = F(t) 191 191 assert allclose( 2*120**2/3 + 60**2/3, q[1] ) 192 192 assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] ) 193 193 194 194 os.remove(filename) 195 195 196 196 197 197 … … 215 215 216 216 #Nice test that may render some of the others redundant. 217 217 218 218 import os, time 219 219 from config import time_format … … 231 231 rectangular(4, 4, 15, 30, origin = (0, -20)) 232 232 233 233 234 234 #print 'Number of elements', len(vertices) 235 235 domain = Domain(points, vertices, boundary) … … 246 246 domain.starttime = start 247 247 248 248 249 249 #Store structure 250 250 domain.initialise_storage() … … 262 262 263 263 f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600) 264 domain.set_quantity('ymomentum', f3) 264 domain.set_quantity('ymomentum', f3) 265 265 266 266 #Store and advance time … … 275 275 276 276 #Set domain.starttime to too early 277 domain.starttime = start - 1 277 domain.starttime = start - 1 278 278 279 279 #Create file function … … 286 286 287 287 #Check that domain.starttime isn't updated if later 288 domain.starttime = start + 1 288 domain.starttime = start + 1 289 289 F = file_function(filename + '.sww', domain, 290 290 quantities = domain.conserved_quantities, 291 291 interpolation_points = interpolation_points) 292 assert allclose(domain.starttime, start+1) 292 assert allclose(domain.starttime, start+1) 293 293 294 294 domain.starttime = start 295 295 296 296 297 297 298 298 #Check linear interpolation in time 299 299 #for id in range(len(interpolation_points)): 300 300 for id in [1]: 301 301 x = interpolation_points[id][0] 302 y = interpolation_points[id][1] 303 302 y = interpolation_points[id][1] 303 304 304 for i in range(20): 305 305 t = i*10 306 306 k = i%6 307 307 308 308 if k == 0: 309 q0 = F(t, point_id=id) 310 q1 = F(t+60, point_id=id) 311 312 309 q0 = F(t, point_id=id) 310 q1 = F(t+60, point_id=id) 311 312 313 313 q = F(t, point_id=id) 314 314 assert allclose(q, (k*q1 + (6-k)*q0)/6) 315 316 315 316 317 317 #Another check of linear interpolation in time 318 318 for id in range(len(interpolation_points)): 319 319 q60 = F(60, point_id=id) 320 q120 = F(120, point_id=id) 321 320 q120 = F(120, point_id=id) 321 322 322 t = 90 #Halfway between 60 and 120 323 323 q = F(t,point_id=id) … … 333 333 #than file end time 334 334 delta = 23 335 domain.starttime = start + delta 335 domain.starttime = start + delta 336 336 F = file_function(filename + '.sww', domain, 337 337 quantities = domain.conserved_quantities, 338 338 interpolation_points = interpolation_points) 339 assert allclose(domain.starttime, start+delta) 340 341 342 343 339 assert allclose(domain.starttime, start+delta) 340 341 342 343 344 344 #Now try interpolation with delta offset 345 345 for id in [1]: 346 346 x = interpolation_points[id][0] 347 y = interpolation_points[id][1] 348 347 y = interpolation_points[id][1] 348 349 349 for i in range(20): 350 350 t = i*10 351 351 k = i%6 352 352 353 353 if k == 0: 354 q0 = F(t-delta, point_id=id) 355 q1 = F(t+60-delta, point_id=id) 356 354 q0 = F(t-delta, point_id=id) 355 q1 = F(t+60-delta, point_id=id) 356 357 357 q = F(t-delta, point_id=id) 358 358 assert allclose(q, (k*q1 + (6-k)*q0)/6) 359 360 359 360 361 361 os.remove(filename + '.sww') 362 362 363 363 364 364 … … 385 385 fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600))) 386 386 t += dt 387 387 388 388 fid.close() 389 389 … … 391 391 b = [4.0, 0.0] 392 392 c = [0.0, 3.0] 393 393 394 394 points = [a, b, c] 395 395 vertices = [[0,1,2]] 396 domain = Domain(points, vertices) 396 domain = Domain(points, vertices) 397 397 398 398 #Check that domain.starttime is updated if non-existing 399 399 F = file_function(filename, domain) 400 assert allclose(domain.starttime, start) 400 assert allclose(domain.starttime, start) 401 401 402 402 #Check that domain.starttime is updated if too early 403 domain.starttime = start - 1 403 domain.starttime = start - 1 404 404 F = file_function(filename, domain) 405 405 assert allclose(domain.starttime, start) 406 406 407 407 #Check that domain.starttime isn't updated if later 408 domain.starttime = start + 1 408 domain.starttime = start + 1 409 409 F = file_function(filename, domain) 410 assert allclose(domain.starttime, start+1) 410 assert allclose(domain.starttime, start+1) 411 411 412 412 domain.starttime = start 413 414 413 414 415 415 #Now try interpolation 416 416 for i in range(20): … … 419 419 420 420 #Exact linear intpolation 421 assert allclose(q[0], 2*t) 421 assert allclose(q[0], 2*t) 422 422 if i%6 == 0: 423 423 assert allclose(q[1], t**2) 424 assert allclose(q[2], sin(t*pi/600)) 424 assert allclose(q[2], sin(t*pi/600)) 425 425 426 426 #Check non-exact 427 427 428 428 t = 90 #Halfway between 60 and 120 429 q = F(t) 429 q = F(t) 430 430 assert allclose( (120**2 + 60**2)/2, q[1] ) 431 431 assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] ) … … 433 433 434 434 t = 100 #Two thirds of the way between between 60 and 120 435 q = F(t) 435 q = F(t) 436 436 assert allclose( 2*120**2/3 + 60**2/3, q[1] ) 437 437 assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] ) 438 438 439 439 os.remove(filename) 440 440 441 441 442 442 def test_file_function_time_with_domain_different_start(self): … … 464 464 fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600))) 465 465 t += dt 466 466 467 467 fid.close() 468 468 … … 470 470 b = [4.0, 0.0] 471 471 c = [0.0, 3.0] 472 472 473 473 points = [a, b, c] 474 474 vertices = [[0,1,2]] 475 domain = Domain(points, vertices) 475 domain = Domain(points, vertices) 476 476 477 477 #Check that domain.starttime isn't updated if later than file starttime but earlier 478 478 #than file end time 479 479 delta = 23 480 domain.starttime = start + delta 480 domain.starttime = start + delta 481 481 F = file_function(filename, domain) 482 assert allclose(domain.starttime, start+delta) 483 484 485 486 487 #Now try interpolation with delta offset 482 assert allclose(domain.starttime, start+delta) 483 484 485 486 487 #Now try interpolation with delta offset 488 488 for i in range(20): 489 489 t = i*10 … … 491 491 492 492 #Exact linear intpolation 493 assert allclose(q[0], 2*t) 493 assert allclose(q[0], 2*t) 494 494 if i%6 == 0: 495 495 assert allclose(q[1], t**2) 496 assert allclose(q[2], sin(t*pi/600)) 496 assert allclose(q[2], sin(t*pi/600)) 497 497 498 498 #Check non-exact 499 499 500 500 t = 90 #Halfway between 60 and 120 501 q = F(t-delta) 501 q = F(t-delta) 502 502 assert allclose( (120**2 + 60**2)/2, q[1] ) 503 503 assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] ) … … 505 505 506 506 t = 100 #Two thirds of the way between between 60 and 120 507 q = F(t-delta) 507 q = F(t-delta) 508 508 assert allclose( 2*120**2/3 + 60**2/3, q[1] ) 509 509 assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] ) 510 510 511 511 os.remove(filename) 512 512 … … 517 517 """ 518 518 import time 519 520 #Create sww file of simple propagation from left to right 519 520 #Create sww file of simple propagation from left to right 521 521 #through rectangular domain 522 522 from shallow_water import Domain, Dirichlet_boundary 523 523 from mesh_factory import rectangular 524 from Numeric import take, concatenate, reshape 524 from Numeric import take, concatenate, reshape 525 525 526 526 #Create basic mesh and shallow water domain 527 527 points, vertices, boundary = rectangular(3, 3) 528 528 domain1 = Domain(points, vertices, boundary) 529 529 530 530 from util import mean 531 domain1.reduction = mean 531 domain1.reduction = mean 532 532 domain1.smooth = True #NOTE: Mimic sww output where each vertex has 533 533 # only one value. 534 534 535 535 domain1.default_order = 2 536 domain1.store = True 536 domain1.store = True 537 537 domain1.set_datadir('.') 538 538 domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self))) 539 539 domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] 540 540 541 541 #Bed-slope, friction and IC at vertices (and interpolated elsewhere) 542 542 domain1.set_quantity('elevation', 0) 543 543 domain1.set_quantity('friction', 0) 544 domain1.set_quantity('stage', 0) 544 domain1.set_quantity('stage', 0) 545 545 546 546 # Boundary conditions 547 B0 = Dirichlet_boundary([0,0,0]) 548 B6 = Dirichlet_boundary([0.6,0,0]) 547 B0 = Dirichlet_boundary([0,0,0]) 548 B6 = Dirichlet_boundary([0.6,0,0]) 549 549 domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0}) 550 550 domain1.check_integrity() … … 560 560 from Scientific.IO.NetCDF import NetCDFFile 561 561 filename = domain1.get_name() + '.' + domain1.format 562 fid = NetCDFFile(filename) 562 fid = NetCDFFile(filename) 563 563 564 564 x = fid.variables['x'][:] 565 y = fid.variables['y'][:] 565 y = fid.variables['y'][:] 566 566 stage = fid.variables['stage'][:] 567 567 xmomentum = fid.variables['xmomentum'][:] … … 581 581 #this timestep are 582 582 r0 = (D[0] + D[1])/2 583 r1 = (D[1] + D[2])/2 584 r2 = (D[2] + D[3])/2 583 r1 = (D[1] + D[2])/2 584 r2 = (D[2] + D[3])/2 585 585 586 586 #And the midpoints are found now 587 587 Dx = take(reshape(x, (16,1)), [0,5,10,15]) 588 Dy = take(reshape(y, (16,1)), [0,5,10,15]) 588 Dy = take(reshape(y, (16,1)), [0,5,10,15]) 589 589 590 590 diag = concatenate( (Dx, Dy), axis=1) … … 597 597 598 598 q = f(timestep/10., point_id=0); assert allclose(r0, q) 599 q = f(timestep/10., point_id=1); assert allclose(r1, q) 599 q = f(timestep/10., point_id=1); assert allclose(r1, q) 600 600 q = f(timestep/10., point_id=2); assert allclose(r2, q) 601 601 … … 613 613 #this timestep are 614 614 r0 = (D[0] + D[1])/2 615 r1 = (D[1] + D[2])/2 616 r2 = (D[2] + D[3])/2 615 r1 = (D[1] + D[2])/2 616 r2 = (D[2] + D[3])/2 617 617 618 618 #Let us see if the file function can find the correct 619 #values 619 #values 620 620 q = f(0, point_id=0); assert allclose(r0, q) 621 621 q = f(0, point_id=1); assert allclose(r1, q) … … 626 626 #Now do it again for a timestep in the middle 627 627 628 timestep = 33 628 timestep = 33 629 629 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1)) 630 630 d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) … … 635 635 #this timestep are 636 636 r0 = (D[0] + D[1])/2 637 r1 = (D[1] + D[2])/2 638 r2 = (D[2] + D[3])/2 637 r1 = (D[1] + D[2])/2 638 r2 = (D[2] + D[3])/2 639 639 640 640 q = f(timestep/10., point_id=0); assert allclose(r0, q) 641 q = f(timestep/10., point_id=1); assert allclose(r1, q) 641 q = f(timestep/10., point_id=1); assert allclose(r1, q) 642 642 q = f(timestep/10., point_id=2); assert allclose(r2, q) 643 643 … … 656 656 #this timestep are 657 657 r0_0 = (D[0] + D[1])/2 658 r1_0 = (D[1] + D[2])/2 659 r2_0 = (D[2] + D[3])/2 658 r1_0 = (D[1] + D[2])/2 659 r2_0 = (D[2] + D[3])/2 660 660 661 661 # … … 669 669 #this timestep are 670 670 r0_1 = (D[0] + D[1])/2 671 r1_1 = (D[1] + D[2])/2 672 r2_1 = (D[2] + D[3])/2 671 r1_1 = (D[1] + D[2])/2 672 r2_1 = (D[2] + D[3])/2 673 673 674 674 # The reference values are 675 675 r0 = (r0_0 + r0_1)/2 676 676 r1 = (r1_0 + r1_1)/2 677 r2 = (r2_0 + r2_1)/2 677 r2 = (r2_0 + r2_1)/2 678 678 679 679 q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q) 680 q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q) 680 q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q) 681 681 q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q) 682 682 … … 688 688 r0 = (r0_0 + 2*r0_1)/3 689 689 r1 = (r1_0 + 2*r1_1)/3 690 r2 = (r2_0 + 2*r2_1)/3 690 r2 = (r2_0 + 2*r2_1)/3 691 691 692 692 #And the file function gives 693 693 q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q) 694 q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q) 694 q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q) 695 695 q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q) 696 696 697 697 fid.close() 698 698 import os 699 699 os.remove(filename) 700 700 701 701 702 702 def test_xya_ascii(self): … … 704 704 FN = 'xyatest' + str(time.time()) + '.xya' 705 705 fid = open(FN, 'w') 706 fid.write(' %s %s %s\n' %('a1', 'a2', 'a3') ) 706 fid.write(' %s %s %s\n' %('a1', 'a2', 'a3') ) 707 707 fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) ) 708 708 fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) ) 709 fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) ) 709 fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) ) 710 710 fid.close() 711 711 712 712 points, attributes = read_xya(FN, format = 'asc') 713 713 … … 715 715 assert allclose(attributes['a1'], [10,30,40.2]) 716 716 assert allclose(attributes['a2'], [20,20,40.3]) 717 assert allclose(attributes['a3'], [30,10,40.4]) 718 717 assert allclose(attributes['a3'], [30,10,40.4]) 718 719 719 os.remove(FN) 720 720 721 721 def test_xya_ascii_w_names(self): 722 722 import time, os 723 723 FN = 'xyatest' + str(time.time()) + '.xya' 724 724 fid = open(FN, 'w') 725 fid.write(' %s %s %s\n' %('a1', 'a2', 'a3') ) 725 fid.write(' %s %s %s\n' %('a1', 'a2', 'a3') ) 726 726 fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) ) 727 727 fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) ) 728 fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) ) 728 fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) ) 729 729 fid.close() 730 730 731 731 points, attributes = read_xya(FN, format = 'asc') 732 732 … … 735 735 assert allclose(attributes['a1'], [10,30,40.2]) 736 736 assert allclose(attributes['a2'], [20,20,40.3]) 737 assert allclose(attributes['a3'], [30,10,40.4]) 738 739 737 assert allclose(attributes['a3'], [30,10,40.4]) 738 739 740 740 os.remove(FN) 741 741 … … 743 743 744 744 745 #Polygon stuff 745 #Polygon stuff 746 746 def test_polygon_function_constants(self): 747 747 p1 = [[0,0], [10,0], [10,10], [0,10]] 748 748 p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] 749 750 f = Polygon_function( [(p1, 1.0)] ) 749 750 f = Polygon_function( [(p1, 1.0)] ) 751 751 z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1 752 752 assert allclose(z, [1,1,0,0]) 753 754 753 754 755 755 f = Polygon_function( [(p2, 2.0)] ) 756 756 z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2 757 assert allclose(z, [2,0,0,2]) 758 759 760 #Combined 757 assert allclose(z, [2,0,0,2]) 758 759 760 #Combined 761 761 f = Polygon_function( [(p1, 1.0), (p2, 2.0)] ) 762 z = f([5, 5, 27, 35], [5, 9, 8, -5]) 763 assert allclose(z, [2,1,0,2]) 764 765 762 z = f([5, 5, 27, 35], [5, 9, 8, -5]) 763 assert allclose(z, [2,1,0,2]) 764 765 766 766 def test_polygon_function_callable(self): 767 """Check that values passed into Polygon_function can be callable 767 """Check that values passed into Polygon_function can be callable 768 768 themselves. 769 769 """ 770 770 p1 = [[0,0], [10,0], [10,10], [0,10]] 771 771 p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] 772 773 f = Polygon_function( [(p1, test_function)] ) 772 773 f = Polygon_function( [(p1, test_function)] ) 774 774 z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1 775 775 assert allclose(z, [10,14,0,0]) 776 777 #Combined 776 777 #Combined 778 778 f = Polygon_function( [(p1, test_function), (p2, 2.0)] ) 779 z = f([5, 5, 27, 35], [5, 9, 8, -5]) 780 assert allclose(z, [2,14,0,2]) 781 782 783 #Combined w default 779 z = f([5, 5, 27, 35], [5, 9, 8, -5]) 780 assert allclose(z, [2,14,0,2]) 781 782 783 #Combined w default 784 784 f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14) 785 z = f([5, 5, 27, 35], [5, 9, 8, -5]) 786 assert allclose(z, [2,14,3.14,2]) 787 788 789 #Combined w default func 790 f = Polygon_function( [(p1, test_function), (p2, 2.0)], 785 z = f([5, 5, 27, 35], [5, 9, 8, -5]) 786 assert allclose(z, [2,14,3.14,2]) 787 788 789 #Combined w default func 790 f = Polygon_function( [(p1, test_function), (p2, 2.0)], 791 791 default = test_function) 792 z = f([5, 5, 27, 35], [5, 9, 8, -5]) 793 assert allclose(z, [2,14,35,2]) 794 795 792 z = f([5, 5, 27, 35], [5, 9, 8, -5]) 793 assert allclose(z, [2,14,35,2]) 794 795 796 796 def test_point_on_line(self): 797 798 #Endpoints first 799 assert point_on_line( 0, 0, 0,0, 1,0 ) 800 assert point_on_line( 1, 0, 0,0, 1,0 ) 801 802 #Then points on line 803 assert point_on_line( 0.5, 0, 0,0, 1,0 ) 804 assert point_on_line( 0, 0.5, 0,1, 0,0 ) 805 assert point_on_line( 1, 0.5, 1,1, 1,0 ) 806 assert point_on_line( 0.5, 0.5, 0,0, 1,1 ) 807 808 #Then points not on line 809 assert not point_on_line( 0.5, 0, 0,1, 1,1 ) 810 assert not point_on_line( 0, 0.5, 0,0, 1,1 ) 811 797 798 #Endpoints first 799 assert point_on_line( 0, 0, 0,0, 1,0 ) 800 assert point_on_line( 1, 0, 0,0, 1,0 ) 801 802 #Then points on line 803 assert point_on_line( 0.5, 0, 0,0, 1,0 ) 804 assert point_on_line( 0, 0.5, 0,1, 0,0 ) 805 assert point_on_line( 1, 0.5, 1,1, 1,0 ) 806 assert point_on_line( 0.5, 0.5, 0,0, 1,1 ) 807 808 #Then points not on line 809 assert not point_on_line( 0.5, 0, 0,1, 1,1 ) 810 assert not point_on_line( 0, 0.5, 0,0, 1,1 ) 811 812 812 #From real example that failed 813 813 assert not point_on_line( 40,50, 40,20, 40,40 ) 814 815 814 815 816 816 #From real example that failed 817 assert not point_on_line( 40,19, 40,20, 40,40 ) 818 819 820 821 817 assert not point_on_line( 40,19, 40,20, 40,40 ) 818 819 820 821 822 822 def test_inside_polygon_main(self): 823 823 824 824 825 825 #Simplest case: Polygon is the unit square 826 826 polygon = [[0,0], [1,0], [1,1], [0,1]] 827 827 828 828 assert inside_polygon( (0.5, 0.5), polygon ) 829 assert not inside_polygon( (0.5, 1.5), polygon ) 830 assert not inside_polygon( (0.5, -0.5), polygon ) 831 assert not inside_polygon( (-0.5, 0.5), polygon ) 832 assert not inside_polygon( (1.5, 0.5), polygon ) 833 829 assert not inside_polygon( (0.5, 1.5), polygon ) 830 assert not inside_polygon( (0.5, -0.5), polygon ) 831 assert not inside_polygon( (-0.5, 0.5), polygon ) 832 assert not inside_polygon( (1.5, 0.5), polygon ) 833 834 834 #Try point on borders 835 assert inside_polygon( (1., 0.5), polygon, closed=True) 836 assert inside_polygon( (0.5, 1), polygon, closed=True) 837 assert inside_polygon( (0., 0.5), polygon, closed=True) 838 assert inside_polygon( (0.5, 0.), polygon, closed=True) 839 840 assert not inside_polygon( (0.5, 1), polygon, closed=False) 841 assert not inside_polygon( (0., 0.5), polygon, closed=False) 842 assert not inside_polygon( (0.5, 0.), polygon, closed=False) 843 assert not inside_polygon( (1., 0.5), polygon, closed=False) 844 845 846 835 assert inside_polygon( (1., 0.5), polygon, closed=True) 836 assert inside_polygon( (0.5, 1), polygon, closed=True) 837 assert inside_polygon( (0., 0.5), polygon, closed=True) 838 assert inside_polygon( (0.5, 0.), polygon, closed=True) 839 840 assert not inside_polygon( (0.5, 1), polygon, closed=False) 841 assert not inside_polygon( (0., 0.5), polygon, closed=False) 842 assert not inside_polygon( (0.5, 0.), polygon, closed=False) 843 assert not inside_polygon( (1., 0.5), polygon, closed=False) 844 845 846 847 847 #From real example (that failed) 848 polygon = [[20,20], [40,20], [40,40], [20,40]] 848 polygon = [[20,20], [40,20], [40,40], [20,40]] 849 849 points = [ [40, 50] ] 850 850 res = inside_polygon(points, polygon) 851 851 assert len(res) == 0 852 853 polygon = [[20,20], [40,20], [40,40], [20,40]] 852 853 polygon = [[20,20], [40,20], [40,40], [20,40]] 854 854 points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ] 855 855 res = inside_polygon(points, polygon) 856 856 assert len(res) == 2 857 857 assert allclose(res, [0,1]) 858 859 860 858 859 860 861 861 #More convoluted and non convex polygon 862 862 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] 863 863 assert inside_polygon( (0.5, 0.5), polygon ) 864 assert inside_polygon( (1, -0.5), polygon ) 864 assert inside_polygon( (1, -0.5), polygon ) 865 865 assert inside_polygon( (1.5, 0), polygon ) 866 867 assert not inside_polygon( (0.5, 1.5), polygon ) 868 assert not inside_polygon( (0.5, -0.5), polygon ) 869 870 866 867 assert not inside_polygon( (0.5, 1.5), polygon ) 868 assert not inside_polygon( (0.5, -0.5), polygon ) 869 870 871 871 #Very convoluted polygon 872 872 polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] 873 assert inside_polygon( (5, 5), polygon ) 874 assert inside_polygon( (17, 7), polygon ) 875 assert inside_polygon( (27, 2), polygon ) 876 assert inside_polygon( (35, -5), polygon ) 873 assert inside_polygon( (5, 5), polygon ) 874 assert inside_polygon( (17, 7), polygon ) 875 assert inside_polygon( (27, 2), polygon ) 876 assert inside_polygon( (35, -5), polygon ) 877 877 assert not inside_polygon( (15, 7), polygon ) 878 878 assert not inside_polygon( (24, 3), polygon ) 879 assert not inside_polygon( (25, -10), polygon ) 880 881 882 879 assert not inside_polygon( (25, -10), polygon ) 880 881 882 883 883 #Another combination (that failed) 884 884 polygon = [[0,0], [10,0], [10,10], [0,10]] 885 assert inside_polygon( (5, 5), polygon ) 886 assert inside_polygon( (7, 7), polygon ) 887 assert not inside_polygon( (-17, 7), polygon ) 888 assert not inside_polygon( (7, 17), polygon ) 889 assert not inside_polygon( (17, 7), polygon ) 890 assert not inside_polygon( (27, 8), polygon ) 891 assert not inside_polygon( (35, -5), polygon ) 892 893 894 895 896 #Multiple polygons 897 885 assert inside_polygon( (5, 5), polygon ) 886 assert inside_polygon( (7, 7), polygon ) 887 assert not inside_polygon( (-17, 7), polygon ) 888 assert not inside_polygon( (7, 17), polygon ) 889 assert not inside_polygon( (17, 7), polygon ) 890 assert not inside_polygon( (27, 8), polygon ) 891 assert not inside_polygon( (35, -5), polygon ) 892 893 894 895 896 #Multiple polygons 897 898 898 polygon = [[0,0], [1,0], [1,1], [0,1], [0,0], 899 899 [10,10], [11,10], [11,11], [10,11], [10,10]] 900 assert inside_polygon( (0.5, 0.5), polygon ) 901 assert inside_polygon( (10.5, 10.5), polygon ) 902 900 assert inside_polygon( (0.5, 0.5), polygon ) 901 assert inside_polygon( (10.5, 10.5), polygon ) 902 903 903 #FIXME: Fails if point is 5.5, 5.5 904 assert not inside_polygon( (0, 5.5), polygon ) 905 906 #Polygon with a hole 904 assert not inside_polygon( (0, 5.5), polygon ) 905 906 #Polygon with a hole 907 907 polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1], 908 908 [0,0], [1,0], [1,1], [0,1], [0,0]] 909 909 910 910 assert inside_polygon( (0, -0.5), polygon ) 911 assert not inside_polygon( (0.5, 0.5), polygon ) 912 913 def test_inside_polygon_vector_version(self): 911 assert not inside_polygon( (0.5, 0.5), polygon ) 912 913 def test_inside_polygon_vector_version(self): 914 914 #Now try the vector formulation returning indices 915 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] 915 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] 916 916 points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] 917 917 res = inside_polygon( points, polygon ) 918 918 919 919 assert allclose( res, [0,1,2] ) 920 920 … … 938 938 for point in points: 939 939 assert inside_polygon(point, polygon) 940 941 942 940 941 942 943 943 #------------------------------------------------------------- 944 944 if __name__ == "__main__": 945 945 #suite = unittest.makeSuite(TestCase,'test_inside_polygon_main') 946 suite = unittest.makeSuite(Test Case,'test')946 suite = unittest.makeSuite(Test_Util,'test') 947 947 runner = unittest.TextTestRunner() 948 948 runner.run(suite) 949 950 951 952 953 949 950 951 952
Note: See TracChangeset
for help on using the changeset viewer.
