Changeset 4665
- Timestamp:
- Aug 8, 2007, 9:39:22 AM (17 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 3 edited
-
fit_interpolate/benchmark_least_squares.py (modified) (1 diff)
-
pmesh/mesh.py (modified) (22 diffs)
-
shallow_water/data_manager.py (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py
r4658 r4665 22 22 import tempfile 23 23 import profile , pstats 24 import tempfile25 24 26 25 from anuga.fit_interpolate.interpolate import Interpolate -
anuga_core/source/anuga/pmesh/mesh.py
r4663 r4665 28 28 29 29 #import load_mesh 30 from anuga.coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE 30 from anuga.coordinate_transforms.geo_reference import Geo_reference, \ 31 DEFAULT_ZONE 31 32 from anuga.utilities.polygon import point_in_polygon 32 33 from anuga.load_mesh.loadASCII import NOMAXAREA, export_mesh_file, \ … … 40 41 import kinds 41 42 except ImportError: 42 # Hand-built mockup of the things we need from the kinds package, since it 43 # Hand-built mockup of the things we need from the kinds package, since it 43 44 # was recently removed from the standard Numeric distro. Some users may 44 45 # not have it by default. … … 89 90 Returns the distance from this point to another 90 91 """ 91 SumOfSquares = ((self.x - OtherPoint.x)**2) + ((self.y - OtherPoint.y)**2) 92 SumOfSquares = ((self.x - OtherPoint.x)**2) + \ 93 ((self.y - OtherPoint.y)**2) 92 94 return math.sqrt(SumOfSquares) 93 95 … … 154 156 155 157 VERTEXSQUARESIDELENGTH = 6 156 def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0, yoffset =0, ): 158 def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0, 159 yoffset =0, ): 157 160 x = scale*(self.x + xoffset) 158 161 y = -1*scale*(self.y + yoffset) # - since for a canvas - is up … … 189 192 HOLECORNERLENGTH = 6 190 193 191 def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, yoffset =0, ): 194 def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, 195 yoffset =0, ): 192 196 x = scale*(self.x + xoffset) 193 197 y = -1*scale*(self.y + yoffset) # - since for a canvas - is up … … 260 264 return len(self.attributes)> 1 261 265 262 def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "black"): 266 def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, 267 colour = "black"): 263 268 """ 264 269 Draw a black cross, returning the objectID … … 301 306 """ 302 307 303 def __init__(self, vertex1, vertex2, vertex3, attribute = None, neighbors = None ): 308 def __init__(self, vertex1, vertex2, vertex3, attribute = None, 309 neighbors = None ): 304 310 """ 305 311 Vertices, the initial arguments, are listed in counterclockwise order. … … 347 353 else: 348 354 if offset == 1: 349 self.vertices = [self.vertices[1],self.vertices[2],self.vertices[0]] 350 self.neighbors = [self.neighbors[1],self.neighbors[2],self.neighbors[0]] 355 self.vertices = [self.vertices[1],self.vertices[2], 356 self.vertices[0]] 357 self.neighbors = [self.neighbors[1],self.neighbors[2], 358 self.neighbors[0]] 351 359 if offset == 2: 352 self.vertices = [self.vertices[2],self.vertices[0],self.vertices[1]] 353 self.neighbors = [self.neighbors[2],self.neighbors[0],self.neighbors[1]] 360 self.vertices = [self.vertices[2],self.vertices[0], 361 self.vertices[1]] 362 self.neighbors = [self.neighbors[2],self.neighbors[0], 363 self.neighbors[1]] 354 364 355 365 def rotate_longest_side(self): … … 395 405 return a+b+c 396 406 397 def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 =None):407 def setNeighbors(self,neighbor1=None, neighbor2=None, neighbor3=None): 398 408 """ 399 409 neighbor1 is the triangle opposite vertex1 and so on. … … 413 423 414 424 415 def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "green"): 425 def draw(self, canvas, tags, scale=1, xoffset=0, yoffset=0, 426 colour="green"): 416 427 """ 417 428 Draw a triangle, returning the objectID … … 448 459 449 460 450 def draw(self, canvas, tags,scale=1 , xoffset=0 , yoffset=0,colour='blue'):461 def draw(self, canvas, tags,scale=1, xoffset=0, yoffset=0,colour='blue'): 451 462 x=[] 452 463 y=[] … … 529 540 User attributes describe the mesh region/segments/vertices/attributes 530 541 531 mesh attributes describe the mesh that is produced eg triangles and vertices. 542 mesh attributes describe the mesh that is produced eg triangles and 543 vertices. 532 544 All point information is relative to the geo_reference passed in 533 545 … … 682 694 return self._addRegion(x, y) 683 695 696 def build_grid(self, vert_rows, vert_columns): 697 """ 698 Build a grid with vert_rows number of vertex rows and 699 vert_columns number if vertex columns 700 701 Grid spacing of 1, the origin is the lower left hand corner. 702 703 FIXME(DSG-DSG) no test. 704 """ 705 706 for i in range(vert_rows): 707 for j in range(vert_columns): 708 self.addUserVertex(j,i) 709 self.auto_segment() 710 self.generateMesh(mode = "Q", minAngle=20.0) 711 684 712 # Depreciated 685 713 def addRegionEN(self, x,y): … … 1141 1169 segconverter) 1142 1170 #print "processed gen",generatedMesh['generatedsegmentmarkerlist'] 1143 #print "generatedtriangleattributelist",generatedMesh['generatedtriangleattributelist']1144 1171 generatedMesh['generatedtriangleattributelist'] = \ 1145 1172 region_ints2strings(generatedMesh['generatedtriangleattributelist'], … … 1670 1697 for vertex in self.getUserVertices(): 1671 1698 points.append((vertex.x,vertex.y)) 1672 self.shape = anuga.alpha_shape.alpha_shape.Alpha_Shape(points, alpha = alpha) 1699 self.shape = anuga.alpha_shape.alpha_shape.Alpha_Shape(points, 1700 alpha=alpha) 1673 1701 1674 1702 … … 2683 2711 tag = parent_segment.tag 2684 2712 offspring = [line] 2685 offspring.extend(self.subtract_line(user_line,line)) 2713 offspring.extend(self.subtract_line(user_line, 2714 line)) 2686 2715 affine_lines.remove(user_line) 2687 2716 for newline in offspring: … … 2695 2724 point_keys[point]=vert 2696 2725 line_vertices.append(vert) 2697 segment = Segment(line_vertices[0],line_vertices[1],tag) 2726 segment = Segment(line_vertices[0], 2727 line_vertices[1],tag) 2698 2728 userSegments[newline]=segment 2699 2729 affine_lines.append(newline) … … 2714 2744 point_keys[point]=vert 2715 2745 line_vertices.append(vert) 2716 segment = Segment(line_vertices[0],line_vertices[1],tag) 2746 segment = Segment(line_vertices[0], 2747 line_vertices[1],tag) 2717 2748 userSegments[line]=segment 2718 2749 affine_lines.append(line) … … 2819 2850 return ((x1-x0)**2-(y1-y0)**2)**0.5 2820 2851 2821 def threshold(self,setName,min=None,max=None,attribute_name ='elevation'):2852 def threshold(self,setName,min=None,max=None,attribute_name='elevation'): 2822 2853 """ 2823 2854 threshold using d … … 3316 3347 """ 3317 3348 3318 new_triangle1 = Triangle(new_point,triangle.vertices[0],triangle.vertices[1],attribute = triangle.attribute, neighbors = None) 3319 new_triangle2 = Triangle(new_point,triangle.vertices[2],triangle.vertices[0],attribute = triangle.attribute, neighbors = None) 3349 new_triangle1 = Triangle(new_point,triangle.vertices[0], 3350 triangle.vertices[1], 3351 attribute = triangle.attribute, 3352 neighbors = None) 3353 new_triangle2 = Triangle(new_point,triangle.vertices[2], 3354 triangle.vertices[0], 3355 attribute = triangle.attribute, 3356 neighbors = None) 3320 3357 3321 3358 new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2) … … 3358 3395 if not self.neighbour is None: 3359 3396 self.neighbour.rotate_longest_side() 3360 self.next_case = self.get_next_case(mesh,self.neighbour,self.current_triangle) 3397 self.next_case = self.get_next_case(mesh,self.neighbour, 3398 self.current_triangle) 3361 3399 if self.case == 'vertex': 3362 3400 self.new_point=self.old_point … … 3393 3431 attributes.append(att) 3394 3432 new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])] 3395 newVertex = Vertex(new_coordinate[0],new_coordinate[1], attributes = attributes) 3433 newVertex = Vertex(new_coordinate[0],new_coordinate[1], 3434 attributes = attributes) 3396 3435 mesh.maxVertexIndex+=1 3397 3436 newVertex.index = mesh.maxVertexIndex -
anuga_core/source/anuga/shallow_water/data_manager.py
r4663 r4665 2175 2175 ascid.close() 2176 2176 fid.close() 2177 2178 2177 return basename_out 2179 2178 … … 4907 4906 If you get both, do a convertion from the old to the new. 4908 4907 4909 If you only get new_origin, the points are absolute, convert to relative 4908 If you only get new_origin, the points are absolute, 4909 convert to relative 4910 4910 4911 4911 if you only get the current_origin the points are relative, store … … 4928 4928 4929 4929 number_of_points = len(points_utm) 4930 volumes = array(volumes) 4930 volumes = array(volumes) 4931 points_utm = array(points_utm) 4931 4932 4932 4933 # given the two geo_refs and the points, do the stuff … … 4951 4952 4952 4953 # This will put the geo ref in the middle 4953 #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.) 4954 #geo_ref=Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.) 4955 4954 4956 x = points[:,0] 4955 4957 y = points[:,1] … … 5122 5124 #Title 5123 5125 fid.write('time' +d+ 'HA depth m'+d+ \ 5124 5126 'UA momentum East x m/sec' +d+ 'VA momentum North y m/sec' \ 5125 5127 + "\n") 5126 5128 for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']): … … 5231 5233 if extra_info <>'': 5232 5234 extra_info = '_'+str(extra_info) 5233 screen_output_name = dir_name + "screen_output%s%s%s.txt" %(myid,numprocs,extra_info) 5234 screen_error_name = dir_name + "screen_error%s%s%s.txt" %(myid,numprocs,extra_info) 5235 screen_output_name = dir_name + "screen_output%s%s%s.txt" %(myid, 5236 numprocs, 5237 extra_info) 5238 screen_error_name = dir_name + "screen_error%s%s%s.txt" %(myid, 5239 numprocs, 5240 extra_info) 5235 5241 print screen_output_name 5236 5242 #used to catch screen output to file … … 5323 5329 might be a better way to do this using CSV module Writer and writeDict 5324 5330 5325 writes file to "output_dir" unless "completed" is in kwargs, then it writes to5326 "file_name" kwargs5331 writes file to "output_dir" unless "completed" is in kwargs, then 5332 it writes to "file_name" kwargs 5327 5333 5328 5334 """ … … 5401 5407 fid.close() 5402 5408 else: 5403 #backup plan, if header is different and has completed will append info to 5409 #backup plan, 5410 # if header is different and has completed will append info to 5404 5411 #end of details_temp.cvs file in output directory 5405 5412 file = str(kwargs['output_dir'])+'detail_temp.csv' … … 5431 5438 verbose=False) 5432 5439 5433 filename is a NetCDF sww file containing ANUGA model output. 5434 Optional arguments polygon and time_interval restricts the maximum runup calculation 5440 filename is a NetCDF sww file containing ANUGA model output. 5441 Optional arguments polygon and time_interval restricts the maximum 5442 runup calculation 5435 5443 to a points that lie within the specified polygon and time interval. 5436 5444 … … 5465 5473 5466 5474 filename is a NetCDF sww file containing ANUGA model output. 5467 Optional arguments polygon and time_interval restricts the maximum runup calculation 5475 Optional arguments polygon and time_interval restricts the maximum 5476 runup calculation 5468 5477 to a points that lie within the specified polygon and time interval. 5469 5478 … … 5495 5504 5496 5505 5497 Algorithm is as in get_maximum_inundation_elevation from shallow_water_domain 5506 Algorithm is as in get_maximum_inundation_elevation from 5507 shallow_water_domain 5498 5508 except that this function works with the sww file and computes the maximal 5499 5509 runup height over multiple timesteps. 5500 5510 5501 Optional arguments polygon and time_interval restricts the maximum runup calculation 5502 to a points that lie within the specified polygon and time interval. Polygon is 5511 Optional arguments polygon and time_interval restricts the 5512 maximum runup calculation 5513 to a points that lie within the specified polygon and time interval. 5514 Polygon is 5503 5515 assumed to be in (absolute) UTM coordinates in the same zone as domain. 5504 5516 … … 5561 5573 5562 5574 5563 # Here's where one could convert nodal information to centroid information 5575 # Here's where one could convert nodal information to centroid 5576 # information 5564 5577 # but is probably something we need to write in C. 5565 5578 # Here's a Python thought which is NOT finished!!! … … 5574 5587 msg = 'polygon must be a sequence of points.' 5575 5588 assert len(polygon[0]) == 2, msg 5576 # FIXME (Ole): Make a generic polygon input check in polygon.py and call it here 5589 # FIXME (Ole): Make a generic polygon input check in polygon.py 5590 # and call it here 5577 5591 5578 5592 points = concatenate((x[:,NewAxis], y[:,NewAxis]), axis=1) … … 5600 5614 msg = 'time_interval must be a sequence of length 2.' 5601 5615 assert len(time_interval) == 2, msg 5602 msg = 'time_interval %s must not be decreasing.' %(time_interval) 5616 msg = 'time_interval %s must not be decreasing.' %(time_interval) 5603 5617 assert time_interval[1] >= time_interval[0], msg 5604 5618 … … 5633 5647 for i in timesteps: 5634 5648 if use_centroid_values is True: 5635 stage_i = get_centroid_values(stage[i,:], volumes) 5649 stage_i = get_centroid_values(stage[i,:], volumes) 5636 5650 else: 5637 5651 stage_i = stage[i,:]
Note: See TracChangeset
for help on using the changeset viewer.
