Changeset 5858 for anuga_work/development/anuga_1d/quantity.py
- Timestamp:
- Oct 22, 2008, 7:36:23 PM (16 years ago)
- File:
-
- 1 edited
-
anuga_work/development/anuga_1d/quantity.py (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/quantity.py
r5844 r5858 16 16 17 17 18 18 19 class Quantity: 19 20 21 20 22 def __init__(self, domain, vertex_values=None): 21 22 #from domain import Domain 23 #from domain_order2 import Domain 23 #Initialise Quantity using optional vertex values. 24 24 25 from domain import Domain 25 26 from Numeric import array, zeros, Float … … 78 79 self.beta = domain.beta 79 80 80 #Methods for operator overloading 81 81 82 def __len__(self): 83 """ 84 Returns number of intervals. 85 """ 82 86 return self.centroid_values.shape[0] 83 87 84 88 def interpolate(self): 85 """Compute interpolated values at centroid 89 """ 90 Compute interpolated values at centroid 86 91 Pre-condition: vertex_values have been set 87 92 """ … … 104 109 In case of location == 'centroid' the dimension values must 105 110 be a list of a Numerical array of length N, N being the number 106 of elements in the mesh. Otherwise it must be of dimension Nx 3111 of elements in the mesh. Otherwise it must be of dimension Nx2 107 112 108 113 The values will be stored in elements following their … … 117 122 """ 118 123 119 if location not in ['vertices', 'centroids']: #, 'edges']:124 if location not in ['vertices', 'centroids']: 120 125 msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location 121 126 raise msg … … 143 148 #Intialise centroid and edge values 144 149 self.interpolate() 150 145 151 146 152 … … 217 223 return X, Compatible list, Numeric array (see below) 218 224 location: Where values are to be stored. 219 Permissible options are: vertices, edges,centroid225 Permissible options are: vertices, centroid 220 226 and unique vertices. Default is 'vertices' 221 227 … … 233 239 from Numeric import take 234 240 235 #if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:236 241 if location not in ['vertices', 'centroids', 'unique vertices']: 237 242 msg = 'Invalid location: %s' %location … … 247 252 indices = range(len(self)) 248 253 return take(self.centroid_values,indices) 249 #elif location == 'edges':250 # if (indices == None):251 # indices = range(len(self))252 # return take(self.edge_values,indices)253 254 elif location == 'unique vertices': 254 255 if (indices == None): … … 257 258 #Go through list of unique vertices 258 259 for unique_vert_id in indices: 259 triangles = self.domain.vertexlist[unique_vert_id]260 cells = self.domain.vertexlist[unique_vert_id] 260 261 261 262 #In case there are unused points 262 if triangles is None:263 msg = 'Unique vertex not associated with triangles'263 if cells is None: 264 msg = 'Unique vertex not associated with cells' 264 265 raise msg 265 266 266 # Go through all triangle, vertex pairs267 # Go through all cells, vertex pairs 267 268 # Average the values 268 269 sum = 0 269 for triangle_id, vertex_id in triangles:270 sum += self.vertex_values[ triangle_id, vertex_id]271 vert_values.append(sum/len( triangles))270 for cell_id, vertex_id in cells: 271 sum += self.vertex_values[cell_id, vertex_id] 272 vert_values.append(sum/len(cells)) 272 273 return Numeric.array(vert_values) 273 274 else: … … 276 277 return take(self.vertex_values,indices) 277 278 278 #Method for outputting model results 279 #FIXME: Split up into geometric and numeric stuff. 280 #FIXME: Geometric (X,Y,V) should live in mesh.py 281 #FIXME: STill remember to move XY to mesh 279 282 280 def get_vertex_values(self, 283 #xy=True,284 281 x=True, 285 282 smooth = None, … … 289 286 290 287 The vertex values are returned as one sequence in the 1D float array A. 291 If requested the coordinates will be returned in 1D arrays X and Y.288 If requested the coordinates will be returned in 1D arrays X. 292 289 293 290 The connectivity is represented as an integer array, V, of dimension 294 M x 3, where M is the number of volumes. Each row has threeindices295 into the X, Y, A arrays defining the triangle.291 M x 2, where M is the number of volumes. Each row has two indices 292 into the X, A arrays defining the element. 296 293 297 294 if smooth is True, vertex values corresponding to one common … … 302 299 If no smoothings is required, vertex coordinates and values will 303 300 be aggregated as a concatenation of values at 304 vertices 0, vertices 1 and vertices 2301 vertices 0, vertices 1 305 302 306 303 307 304 Calling convention 308 if x yis True:309 X, Y,A,V = get_vertex_values305 if x is True: 306 X,A,V = get_vertex_values 310 307 else: 311 308 A,V = get_vertex_values … … 351 348 A[k] = reduction(contributions) 352 349 353 354 #if xy is True:355 350 if x is True: 356 351 #X = self.domain.coordinates[:,0].astype(precision) … … 380 375 381 376 #Do vertex coordinates 382 #if xy is True:383 377 if x is True: 384 C= self.domain.get_vertex_coordinates()385 386 X = C[:,0:6:2].copy()387 Y = C[:,1:6:2].copy()388 389 return X.flat, Y.flat,A, V378 X = self.domain.get_vertex_coordinates() 379 380 #X = C[:,0:6:2].copy() 381 #Y = C[:,1:6:2].copy() 382 383 return X.flat, A, V 390 384 else: 391 385 return A, V … … 412 406 N = self.centroid_values.shape[0] 413 407 414 #print "pre update",self.centroid_values415 #print "timestep",timestep416 #print min(self.domain.areas)417 #print "explicit_update", self.explicit_update418 419 408 #Explicit updates 420 409 self.centroid_values += timestep*self.explicit_update 421 422 #print "post update",self.centroid_values423 410 424 411 #Semi implicit updates … … 438 425 """ 439 426 427 #print 'compute_gradient' 440 428 441 429 from Numeric import array, zeros, Float … … 457 445 k0 = k 458 446 k1 = k+1 447 k2 = k+2 459 448 460 449 q0 = Q[k0] 461 450 q1 = Q[k1] 451 q2 = Q[k2] 462 452 463 453 x0 = X[k0] #V0 centroid 464 454 x1 = X[k1] #V1 centroid 455 x2 = X[k2] 465 456 466 457 #Gradient 467 G[k] = (q1 - q0)/(x1 - x0) 458 #G[k] = (q1 - q0)/(x1 - x0) 459 460 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 461 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 468 462 469 463 elif k == N-1: … … 472 466 k0 = k 473 467 k1 = k-1 468 k2 = k-2 474 469 475 470 q0 = Q[k0] 476 471 q1 = Q[k1] 472 q2 = Q[k2] 477 473 478 474 x0 = X[k0] #V0 centroid 479 475 x1 = X[k1] #V1 centroid 476 x2 = X[k2] 480 477 481 478 #Gradient 482 G[k] = (q1 - q0)/(x1 - x0) 479 #G[k] = (q1 - q0)/(x1 - x0) 480 481 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 482 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 483 484 ## q0 = Q[k0] 485 ## q1 = Q[k1] 486 ## 487 ## x0 = X[k0] #V0 centroid 488 ## x1 = X[k1] #V1 centroid 489 ## 490 ## #Gradient 491 ## G[k] = (q1 - q0)/(x1 - x0) 483 492 484 493 else: … … 507 516 """ 508 517 518 #print 'compute_minmod_gradients' 519 509 520 from Numeric import array, zeros, Float,sign 510 521 … … 533 544 k0 = k 534 545 k1 = k+1 546 k2 = k+2 535 547 536 548 q0 = Q[k0] 537 549 q1 = Q[k1] 550 q2 = Q[k2] 538 551 539 552 x0 = X[k0] #V0 centroid 540 553 x1 = X[k1] #V1 centroid 554 x2 = X[k2] 541 555 542 556 #Gradient 543 G[k] = (q1 - q0)/(x1 - x0) 557 #G[k] = (q1 - q0)/(x1 - x0) 558 559 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 560 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 544 561 545 562 elif k == N-1: … … 548 565 k0 = k 549 566 k1 = k-1 567 k2 = k-2 550 568 551 569 q0 = Q[k0] 552 570 q1 = Q[k1] 571 q2 = Q[k2] 553 572 554 573 x0 = X[k0] #V0 centroid 555 574 x1 = X[k1] #V1 centroid 575 x2 = X[k2] 556 576 557 577 #Gradient 558 G[k] = (q1 - q0)/(x1 - x0) 578 #G[k] = (q1 - q0)/(x1 - x0) 579 580 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 581 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 582 583 ## #Get data 584 ## k0 = k 585 ## k1 = k-1 586 ## 587 ## q0 = Q[k0] 588 ## q1 = Q[k1] 589 ## 590 ## x0 = X[k0] #V0 centroid 591 ## x1 = X[k1] #V1 centroid 592 ## 593 ## #Gradient 594 ## G[k] = (q1 - q0)/(x1 - x0) 559 595 560 596 elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2): … … 732 768 #print limiter 733 769 770 #print 'limit_range' 734 771 N = self.domain.number_of_elements 735 772 qc = self.centroid_values … … 866 903 # Call correct module function 867 904 # (either from this module or C-extension) 868 #saxpy_centroid_values(self,a,b)869 905 self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f') 870 906 … … 883 919 884 920 921 ## 922 ##def newLinePlot(title='Simple Plot'): 923 ## import Gnuplot 924 ## g = Gnuplot.Gnuplot() 925 ## g.title(title) 926 ## g('set data style linespoints') 927 ## g.xlabel('x') 928 ## g.ylabel('y') 929 ## return g 930 ## 931 ##def linePlot(g,x,y): 932 ## import Gnuplot 933 ## g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat)) 885 934 886 935 def newLinePlot(title='Simple Plot'): 887 import Gnuplot 888 g = Gnuplot.Gnuplot() 936 import pylab as g 937 g.ion() 938 g.hold(False) 889 939 g.title(title) 890 g('set data style linespoints')891 940 g.xlabel('x') 892 941 g.ylabel('y') 893 return g 894 895 def linePlot(g,x,y): 896 import Gnuplot 897 g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat)) 898 899 900 901 902 942 943 944 def linePlot(x,y): 945 import pylab as g 946 g.plot(x.flat,y.flat) 947 948 949 def closePlots(): 950 import pylab as g 951 g.close('all') 952 903 953 if __name__ == "__main__": 904 954 #from domain import Domain … … 950 1000 Qc[1,0] = 3 951 1001 952 Q1.beta = 1.0953 1002 Q1.extrapolate_second_order() 954 Q1.limit_minmod() 955 956 g1 = newLinePlot('plot 1') 957 linePlot(g1,Xc,Qc) 1003 #Q1.limit_minmod() 1004 1005 newLinePlot('plots') 1006 linePlot(Xc,Qc) 1007 raw_input('press return') 958 1008 959 1009 points2 = arange(10) … … 964 1014 Xc = Q2.domain.vertices 965 1015 Qc = Q2.vertex_values 966 967 g2 = newLinePlot('plot 2') 968 linePlot(g2,Xc,Qc) 969 1016 linePlot(Xc,Qc) 1017 raw_input('press return') 970 1018 971 1019 972 1020 Q2.extrapolate_second_order() 973 Q2.limit_minmod()1021 #Q2.limit_minmod() 974 1022 Xc = Q2.domain.vertices 975 1023 Qc = Q2.vertex_values 976 977 1024 print Q2.centroid_values 978 1025 print Qc 979 raw_input('press_return') 980 981 g3 = newLinePlot('plot 3') 982 linePlot(g3,Xc,Qc) 1026 linePlot(Xc,Qc) 983 1027 raw_input('press return') 984 1028 985 1029 986 1030 for i in range(10): 1031 import pylab as g 1032 g.hold(True) 987 1033 fun = FunClass(i/10.0) 988 Q2.set_values(fun,'vertices') 1034 Q2.set_values(fun,'centroids') 1035 Q2.extrapolate_second_order() 1036 #Q2.limit_minmod() 989 1037 Qc = Q2.vertex_values 990 linePlot(g3,Xc,Qc) 991 992 raw_input('press return') 993 1038 linePlot(Xc,Qc) 1039 raw_input('press return') 1040 1041 raw_input('press return to quit') 1042 closePlots()
Note: See TracChangeset
for help on using the changeset viewer.
