Changeset 7276
- Timestamp:
- Jun 30, 2009, 2:07:41 PM (14 years ago)
- Files:
-
- 1 deleted
- 201 edited
- 181 copied
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r7053 r7276 26 26 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ 27 27 import Transmissive_boundary 28 29 28 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain 30 29 from anuga.abstract_2d_finite_volumes.region\ … … 34 33 from quantity import Quantity 35 34 36 import Numericas num35 import numpy as num 37 36 38 37 … … 95 94 """ 96 95 97 98 96 number_of_full_nodes=None 99 97 number_of_full_triangles=None … … 193 191 buffer_shape = self.full_send_dict[key][0].shape[0] 194 192 self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys), 195 num. Float))193 num.float)) 196 194 197 195 for key in self.ghost_recv_dict: 198 196 buffer_shape = self.ghost_recv_dict[key][0].shape[0] 199 197 self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys), 200 num. Float))198 num.float)) 201 199 202 200 # Setup cell full flag … … 205 203 N = len(self) #number_of_elements 206 204 self.number_of_elements = N 207 self.tri_full_flag = num.ones(N, num. Int)205 self.tri_full_flag = num.ones(N, num.int) 208 206 for i in self.ghost_recv_dict.keys(): 209 207 for id in self.ghost_recv_dict[i][0]: … … 267 265 # (boolean) array, to be used during the flux calculation. 268 266 N = len(self) # Number_of_triangles 269 self.already_computed_flux = num.zeros((N, 3), num. Int)267 self.already_computed_flux = num.zeros((N, 3), num.int) 270 268 271 269 # Storage for maximal speeds computed for each triangle by 272 270 # compute_fluxes. 273 271 # This is used for diagnostics only (reset at every yieldstep) 274 self.max_speed = num.zeros(N, num. Float)272 self.max_speed = num.zeros(N, num.float) 275 273 276 274 if mesh_filename is not None: … … 391 389 raise Exception, msg 392 390 393 q = num.zeros(len(self.conserved_quantities), num. Float)391 q = num.zeros(len(self.conserved_quantities), num.float) 394 392 395 393 for i, name in enumerate(self.conserved_quantities): … … 463 461 464 462 name: Name of quantity 465 value: Compatible list, Numeric array, const or function (see below)463 value: Compatible list, numeric array, const or function (see below) 466 464 467 465 The values will be stored in elements following their internal ordering. … … 895 893 896 894 msg += '------------------------------------------------\n' 897 msg += ' Speeds in [%f, %f]\n' % ( min(self.max_speed),898 max(self.max_speed))895 msg += ' Speeds in [%f, %f]\n' % (num.min(self.max_speed), 896 num.max(self.max_speed)) 899 897 msg += ' Histogram:\n' 900 898 … … 908 906 else: 909 907 # Closed upper interval 910 hi = max(self.max_speed)908 hi = num.max(self.max_speed) 911 909 msg += ' [%f, %f]: %d\n' % (lo, hi, count) 912 910 913 N = len(self.max_speed )911 N = len(self.max_speed.flat) 914 912 if N > 10: 915 913 msg += ' Percentiles (10%):\n' … … 1373 1371 self.number_of_steps = 0 1374 1372 self.number_of_first_order_steps = 0 1375 self.max_speed = num.zeros(N, num. Float)1373 self.max_speed = num.zeros(N, num.float) 1376 1374 1377 1375 ## … … 1642 1640 # momentum 1643 1641 if self.protect_against_isolated_degenerate_timesteps is True and \ 1644 self.max_speed> 10.0: # FIXME (Ole): Make this configurable1642 num.max(self.max_speed) > 10.0: # FIXME (Ole): Make this configurable 1645 1643 1646 1644 # Setup 10 bins for speed histogram … … 1657 1655 1658 1656 # Find triangles in last bin 1659 # FIXME - speed up using Numeric1657 # FIXME - speed up using numeric package 1660 1658 d = 0 1661 1659 for i in range(self.number_of_full_triangles): … … 1763 1761 # We have a list with ghosts expecting updates 1764 1762 1765 1766 from Numeric import take,put1767 1768 1769 1763 #Update of ghost cells 1770 1764 iproc = self.processor … … 1779 1773 for i, q in enumerate(self.conserved_quantities): 1780 1774 Q_cv = self.quantities[q].centroid_values 1781 put(Q_cv, Idg, take(Q_cv, Idf))1775 num.put(Q_cv, Idg, num.take(Q_cv, Idf, axis=0)) 1782 1776 1783 1777 -
anuga_core/source/anuga/abstract_2d_finite_volumes/ermapper_grids.py
r6161 r7276 2 2 3 3 # from os import open, write, read 4 import Numericas num5 6 celltype_map = {'IEEE4ByteReal': num. Float32, 'IEEE8ByteReal': num.Float64}4 import numpy as num 5 6 celltype_map = {'IEEE4ByteReal': num.float32, 'IEEE8ByteReal': num.float64} 7 7 8 8 … … 11 11 write_ermapper_grid(ofile, data, header = {}): 12 12 13 Function to write a 2D Numeric array to an ERMapper grid. There are a series of conventions adopted within13 Function to write a 2D numeric array to an ERMapper grid. There are a series of conventions adopted within 14 14 this code, specifically: 15 15 1) The registration coordinate for the data is the SW (or lower-left) corner of the data 16 16 2) The registration coordinates refer to cell centres 17 3) The data is a 2D Numeric array with the NW-most data in element (0,0) and the SE-most data in element (N,M)17 3) The data is a 2D numeric array with the NW-most data in element (0,0) and the SE-most data in element (N,M) 18 18 where N is the last line and M is the last column 19 19 4) There has been no testng of the use of a rotated grid. Best to keep data in an NS orientation … … 163 163 return header 164 164 165 def write_ermapper_data(grid, ofile, data_format = num.Float32):165 def write_ermapper_data(grid, ofile, data_format=num.float32): 166 166 167 167 … … 193 193 194 194 195 def read_ermapper_data(ifile, data_format = num. Float32):195 def read_ermapper_data(ifile, data_format = num.float32): 196 196 # open input file in a binary format and read the input string 197 197 fid = open(ifile,'rb') … … 199 199 fid.close() 200 200 201 # convert input string to required format (Note default format is num. Float32)201 # convert input string to required format (Note default format is num.float32) 202 202 grid_as_float = num.fromstring(input_string,data_format) 203 203 return grid_as_float -
anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
r6654 r7276 1 import Numeric as num 1 import copy 2 import numpy as num 2 3 3 4 from anuga.coordinate_transforms.geo_reference import Geo_reference … … 10 11 which is represented as a coordinate set (x,y). 11 12 Vertices from different triangles can point to the same node. 12 The nodes are implemented as an Nx2 Numeric array containing the13 The nodes are implemented as an Nx2 numeric array containing the 13 14 x and y coordinates. 14 15 … … 19 20 where 20 21 21 nodes is either a list of 2-tuples or an Nx2 Numeric array of22 nodes is either a list of 2-tuples or an Nx2 numeric array of 22 23 floats representing all x, y coordinates in the mesh. 23 24 24 triangles is either a list of 3-tuples or an Mx3 Numeric array of25 triangles is either a list of 3-tuples or an Mx3 numeric array of 25 26 integers representing indices of all vertices in the mesh. 26 27 Each vertex is identified by its index i in [0, N-1]. … … 37 38 triangles = [ [1,0,2], [1,2,3] ] # bac, bce 38 39 39 # Create mesh with two triangles: bac and bce 40 # Create mesh with two triangles: bac and bce 40 41 mesh = Mesh(nodes, triangles) 41 42 … … 56 57 """ 57 58 58 #FIXME: It would be a good idea to use geospatial data as an alternative 59 #input 60 def __init__(self, nodes, triangles, 61 geo_reference=None, 59 # FIXME: It would be a good idea to use geospatial data as an alternative 60 # input 61 def __init__(self, 62 nodes, 63 triangles, 64 geo_reference=None, 62 65 number_of_full_nodes=None, 63 number_of_full_triangles=None, 66 number_of_full_triangles=None, 64 67 verbose=False): 65 68 """Build triangular 2d mesh from nodes and triangle information 66 69 67 70 Input: 68 71 69 72 nodes: x,y coordinates represented as a sequence of 2-tuples or 70 a Nx2 Numeric array of floats.71 72 triangles: sequence of 3-tuples or Mx3 Numeric array of73 a Nx2 numeric array of floats. 74 75 triangles: sequence of 3-tuples or Mx3 numeric array of 73 76 non-negative integers representing indices into 74 77 the nodes array. 75 78 76 79 georeference (optional): If specified coordinates are 77 80 assumed to be relative to this origin. … … 83 86 In this case it is usefull to specify the number of real (called full) 84 87 nodes and triangles. If omitted they will default to all. 85 86 """ 87 88 if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' 89 90 self.triangles = num.array(triangles, num.Int) 91 self.nodes = num.array(nodes, num.Float) 92 93 94 # Register number of elements and nodes 88 89 """ 90 91 if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' 92 93 self.triangles = num.array(triangles, num.int) 94 self.nodes = num.array(nodes, num.float) 95 96 # Register number of elements and nodes 95 97 self.number_of_triangles = N = self.triangles.shape[0] 96 self.number_of_nodes = self.nodes.shape[0] 97 98 98 self.number_of_nodes = self.nodes.shape[0] 99 99 100 100 if number_of_full_nodes is None: … … 102 102 else: 103 103 assert int(number_of_full_nodes) 104 self.number_of_full_nodes = number_of_full_nodes 105 104 self.number_of_full_nodes = number_of_full_nodes 106 105 107 106 if number_of_full_triangles is None: 108 self.number_of_full_triangles = self.number_of_triangles 107 self.number_of_full_triangles = self.number_of_triangles 109 108 else: 110 assert int(number_of_full_triangles) 109 assert int(number_of_full_triangles) 111 110 self.number_of_full_triangles = number_of_full_triangles 112 113 114 #print self.number_of_full_nodes, self.number_of_nodes115 #print self.number_of_full_triangles, self.number_of_triangles116 117 118 111 119 112 # FIXME: this stores a geo_reference, but when coords are returned 120 113 # This geo_ref is not taken into account! 121 114 if geo_reference is None: 122 self.geo_reference = Geo_reference() #Use defaults115 self.geo_reference = Geo_reference() # Use defaults 123 116 else: 124 117 self.geo_reference = geo_reference 125 118 126 119 # Input checks 127 msg = 'Triangles must an Mx3 Numeric array or a sequence of 3-tuples. '128 msg += 'The supplied array has the shape: %s'\129 % str(self.triangles.shape)120 msg = ('Triangles must an Mx3 numeric array or a sequence of 3-tuples. ' 121 'The supplied array has the shape: %s' 122 % str(self.triangles.shape)) 130 123 assert len(self.triangles.shape) == 2, msg 131 124 132 msg = 'Nodes must an Nx2 Numeric array or a sequence of 2-tuples' 133 msg += 'The supplied array has the shape: %s'\ 134 %str(self.nodes.shape) 125 msg = ('Nodes must an Nx2 numeric array or a sequence of 2-tuples' 126 'The supplied array has the shape: %s' % str(self.nodes.shape)) 135 127 assert len(self.nodes.shape) == 2, msg 136 128 137 129 msg = 'Vertex indices reference non-existing coordinate sets' 138 assert max(self.triangles.flat) < self.nodes.shape[0], msg 139 130 assert num.max(self.triangles) < self.nodes.shape[0], msg 140 131 141 132 # FIXME: Maybe move to statistics? 142 133 # Or use with get_extent 143 xy_extent = [ min(self.nodes[:,0]), min(self.nodes[:,1]) , 144 max(self.nodes[:,0]), max(self.nodes[:,1]) ] 145 146 self.xy_extent = num.array(xy_extent, num.Float) 147 134 xy_extent = [min(self.nodes[:,0]), min(self.nodes[:,1]), 135 max(self.nodes[:,0]), max(self.nodes[:,1])] 136 137 self.xy_extent = num.array(xy_extent, num.float) 148 138 149 139 # Allocate space for geometric quantities 150 self.normals = num.zeros((N, 6), num. Float)151 self.areas = num.zeros(N, num. Float)152 self.edgelengths = num.zeros((N, 3), num. Float)140 self.normals = num.zeros((N, 6), num.float) 141 self.areas = num.zeros(N, num.float) 142 self.edgelengths = num.zeros((N, 3), num.float) 153 143 154 144 # Get x,y coordinates for all triangles and store 155 145 self.vertex_coordinates = V = self.compute_vertex_coordinates() 156 146 157 158 147 # Initialise each triangle 159 148 if verbose: 160 print 'General_mesh: Computing areas, normals and edgeleng hts'161 149 print 'General_mesh: Computing areas, normals and edgelengths' 150 162 151 for i in range(N): 163 if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' % (i, N)164 152 if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' % (i, N) 153 165 154 x0, y0 = V[3*i, :] 166 155 x1, y1 = V[3*i+1, :] 167 x2, y2 = V[3*i+2, :] 156 x2, y2 = V[3*i+2, :] 168 157 169 158 # Area 170 self.areas[i] = abs((x1*y0-x0*y1) +(x2*y1-x1*y2)+(x0*y2-x2*y0))/2171 172 msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' % (x0,y0,x1,y1,x2,y2)173 msg += ' is degenerate: area == %f' % self.areas[i]159 self.areas[i] = abs((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2 160 161 msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' % (x0,y0,x1,y1,x2,y2) 162 msg += ' is degenerate: area == %f' % self.areas[i] 174 163 assert self.areas[i] > 0.0, msg 175 176 164 177 165 # Normals … … 184 172 # the first vertex, etc) 185 173 # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle 186 187 n0 = num.array([x2 - x1, y2 - y1], num.Float) 174 n0 = num.array([x2-x1, y2-y1], num.float) 188 175 l0 = num.sqrt(num.sum(n0**2)) 189 176 190 n1 = num.array([x0 - x2, y0 - y2], num.Float)177 n1 = num.array([x0-x2, y0-y2], num.float) 191 178 l1 = num.sqrt(num.sum(n1**2)) 192 179 193 n2 = num.array([x1 - x0, y1 - y0], num.Float)180 n2 = num.array([x1-x0, y1-y0], num.float) 194 181 l2 = num.sqrt(num.sum(n2**2)) 195 182 … … 207 194 self.edgelengths[i, :] = [l0, l1, l2] 208 195 209 210 # Build structure listing which trianglse belong to which node. 211 if verbose: print 'Building inverted triangle structure' 196 # Build structure listing which triangles belong to which node. 197 if verbose: print 'Building inverted triangle structure' 212 198 self.build_inverted_triangle_structure() 213 214 215 199 216 200 def __len__(self): 217 201 return self.number_of_triangles 218 219 202 220 203 def __repr__(self): 221 return 'Mesh: %d vertices, %d triangles'\222 %(self.nodes.shape[0], len(self))204 return ('Mesh: %d vertices, %d triangles' 205 % (self.nodes.shape[0], len(self))) 223 206 224 207 def get_normals(self): 225 208 """Return all normal vectors. 209 226 210 Return normal vectors for all triangles as an Nx6 array 227 211 (ordered as x0, y0, x1, y1, x2, y2 for each triangle) 228 212 """ 213 229 214 return self.normals 230 231 215 232 216 def get_normal(self, i, j): 233 217 """Return normal vector j of the i'th triangle. 218 234 219 Return value is the numeric array slice [x, y] 235 220 """ 221 236 222 return self.normals[i, 2*j:2*j+2] 237 223 … … 245 231 def get_number_of_nodes(self): 246 232 return self.number_of_nodes 247 233 248 234 def get_nodes(self, absolute=False): 249 235 """Return all nodes in mesh. … … 256 242 are to be made absolute by taking georeference into account 257 243 Default is False as many parts of ANUGA expects relative coordinates. 258 (To see which, switch to default absolute=True and run tests). 244 (To see which, switch to default absolute=True and run tests). 259 245 """ 260 246 … … 264 250 if not self.geo_reference.is_absolute(): 265 251 V = self.geo_reference.get_absolute(V) 266 252 267 253 return V 268 269 def get_node(self, i, 270 absolute=False): 254 255 def get_node(self, i, absolute=False): 271 256 """Return node coordinates for triangle i. 272 257 … … 274 259 are to be made absolute by taking georeference into account 275 260 Default is False as many parts of ANUGA expects relative coordinates. 276 (To see which, switch to default absolute=True and run tests). 277 """ 278 279 261 (To see which, switch to default absolute=True and run tests). 262 263 Note: This method returns a modified _copy_ of the nodes slice if 264 absolute is True. If absolute is False, just return the slice. 265 This is related to the ensure_numeric() returning a copy problem. 266 """ 267 280 268 V = self.nodes[i,:] 281 269 if absolute is True: 282 270 if not self.geo_reference.is_absolute(): 283 return V + num.array([self.geo_reference.get_xllcorner(), 284 self.geo_reference.get_yllcorner()], num.Float) 285 else: 286 return V 287 else: 288 return V 289 290 291 292 def get_vertex_coordinates(self, 293 triangle_id=None, 294 absolute=False): 295 """Return vertex coordinates for all triangles. 296 271 # get a copy so as not to modify the internal self.nodes array 272 V = copy.copy(V) 273 V += num.array([self.geo_reference.get_xllcorner(), 274 self.geo_reference.get_yllcorner()], num.float) 275 return V 276 277 def get_vertex_coordinates(self, triangle_id=None, absolute=False): 278 """Return vertex coordinates for all triangles. 279 297 280 Return all vertex coordinates for all triangles as a 3*M x 2 array 298 281 where the jth vertex of the ith triangle is located in row 3*i+j and … … 301 284 if triangle_id is specified (an integer) the 3 vertex coordinates 302 285 for triangle_id are returned. 303 286 304 287 Boolean keyword argument absolute determines whether coordinates 305 288 are to be made absolute by taking georeference into account 306 289 Default is False as many parts of ANUGA expects relative coordinates. 307 290 """ 308 291 309 292 V = self.vertex_coordinates 310 293 311 if triangle_id is None: 294 if triangle_id is None: 312 295 if absolute is True: 313 296 if not self.geo_reference.is_absolute(): 314 297 V = self.geo_reference.get_absolute(V) 315 316 298 return V 317 299 else: … … 320 302 assert int(i) == i, msg 321 303 assert 0 <= i < self.number_of_triangles 322 323 i3 = 3*i 304 305 i3 = 3*i 324 306 if absolute is True and not self.geo_reference.is_absolute(): 325 307 offset=num.array([self.geo_reference.get_xllcorner(), 326 self.geo_reference.get_yllcorner()], num. Float)308 self.geo_reference.get_yllcorner()], num.float) 327 309 return num.array([V[i3,:]+offset, 328 310 V[i3+1,:]+offset, 329 V[i3+2,:]+offset], num. Float)311 V[i3+2,:]+offset], num.float) 330 312 else: 331 return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.Float) 332 333 313 return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.float) 334 314 335 315 def get_vertex_coordinate(self, i, j, absolute=False): … … 340 320 msg = 'vertex id j must be an integer in [0,1,2]' 341 321 assert j in [0,1,2], msg 342 343 V = self.get_vertex_coordinates(triangle_id=i, 344 absolute=absolute) 322 323 V = self.get_vertex_coordinates(triangle_id=i, absolute=absolute) 345 324 return V[j,:] 346 347 348 325 349 326 def compute_vertex_coordinates(self): … … 356 333 357 334 M = self.number_of_triangles 358 vertex_coordinates = num.zeros((3*M, 2), num. Float)335 vertex_coordinates = num.zeros((3*M, 2), num.float) 359 336 360 337 for i in range(M): … … 365 342 return vertex_coordinates 366 343 367 368 369 344 def get_triangles(self, indices=None): 370 345 """Get mesh triangles. … … 382 357 if indices is None: 383 358 return self.triangles 384 #indices = range(M) 385 386 return num.take(self.triangles, indices) 387 388 359 360 return num.take(self.triangles, indices, axis=0) 389 361 390 362 def get_disconnected_triangles(self): … … 405 377 406 378 The triangles created will have the format 407 408 379 [[0,1,2], 409 380 [3,4,5], 410 381 [6,7,8], 411 382 ... 412 [3*M-3 3*M-2 3*M-1]] 383 [3*M-3 3*M-2 3*M-1]] 413 384 """ 414 385 415 386 M = len(self) # Number of triangles 416 387 K = 3*M # Total number of unique vertices 417 418 T = num.reshape(num.arange(K, typecode=num.Int), (M,3)) 419 420 return T 421 422 423 424 def get_unique_vertices(self, indices=None): 425 """FIXME(Ole): This function needs a docstring 426 """ 388 return num.reshape(num.arange(K, dtype=num.int), (M,3)) 389 390 def get_unique_vertices(self, indices=None): 391 """FIXME(Ole): This function needs a docstring""" 392 427 393 triangles = self.get_triangles(indices=indices) 428 394 unique_verts = {} 429 395 for triangle in triangles: 430 unique_verts[triangle[0]] = 0431 unique_verts[triangle[1]] = 0432 unique_verts[triangle[2]] = 0396 unique_verts[triangle[0]] = 0 397 unique_verts[triangle[1]] = 0 398 unique_verts[triangle[2]] = 0 433 399 return unique_verts.keys() 434 435 400 436 401 def get_triangles_and_vertices_per_node(self, node=None): … … 447 412 # Get index for this node 448 413 first = num.sum(self.number_of_triangles_per_node[:node]) 449 414 450 415 # Get number of triangles for this node 451 416 count = self.number_of_triangles_per_node[node] … … 459 424 triangle_list.append( (volume_id, vertex_id) ) 460 425 461 triangle_list = num.array(triangle_list, num. Int) #array default#426 triangle_list = num.array(triangle_list, num.int) #array default# 462 427 else: 463 428 # Get info for all nodes recursively. … … 465 430 # working directly with the inverted triangle structure 466 431 for i in range(self.number_of_full_nodes): 467 468 432 L = self.get_triangles_and_vertices_per_node(node=i) 469 470 433 triangle_list.append(L) 471 434 472 435 return triangle_list 473 474 436 475 437 def build_inverted_triangle_structure(self): … … 481 443 listing for each node how many triangles use it. N is the number of 482 444 nodes in mesh. 483 445 484 446 vertex_value_indices: An array of length M listing indices into 485 447 triangles ordered by node number. The (triangle_id, vertex_id) … … 489 451 used to average vertex values efficiently. 490 452 491 492 453 Example: 493 494 454 a = [0.0, 0.0] # node 0 495 455 b = [0.0, 2.0] # node 1 … … 498 458 e = [2.0, 2.0] # node 4 499 459 f = [4.0, 0.0] # node 5 500 501 460 nodes = array([a, b, c, d, e, f]) 502 503 #bac, bce, ecf, dbe, daf, dae 504 triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 505 506 507 For this structure 508 461 462 # bac, bce, ecf, dbe 463 triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 464 465 For this structure: 509 466 number_of_triangles_per_node = [1 3 3 1 3 1] 510 467 which means that node a has 1 triangle associated with it, node b 511 468 has 3, node has 3 and so on. 512 469 513 470 vertex_value_indices = [ 1 0 3 10 2 4 7 9 5 6 11 8] 514 471 which reflects the fact that … … 524 481 and by triangle 2, vertex 0 (index = 6) 525 482 and by triangle 3, vertex 2 (index = 11) 526 node 5 is used by triangle 2, vertex 2 (index = 8) 527 483 node 5 is used by triangle 2, vertex 2 (index = 8) 528 484 529 485 Preconditions: … … 532 488 Postcondition: 533 489 self.number_of_triangles_per_node is built 534 self.vertex_value_indices is built 490 self.vertex_value_indices is built 535 491 """ 536 492 537 493 # Count number of triangles per node 538 number_of_triangles_per_node = num.zeros(self.number_of_full_nodes) 494 number_of_triangles_per_node = num.zeros(self.number_of_full_nodes, 495 num.int) #array default# 539 496 for volume_id, triangle in enumerate(self.get_triangles()): 540 497 for vertex_id in triangle: … … 543 500 # Allocate space for inverted structure 544 501 number_of_entries = num.sum(number_of_triangles_per_node) 545 vertex_value_indices = num.zeros(number_of_entries )502 vertex_value_indices = num.zeros(number_of_entries, num.int) #array default# 546 503 547 504 # Register (triangle, vertex) indices for each node 548 vertexlist = [None] *self.number_of_full_nodes505 vertexlist = [None] * self.number_of_full_nodes 549 506 for volume_id in range(self.number_of_full_triangles): 550 551 507 a = self.triangles[volume_id, 0] 552 508 b = self.triangles[volume_id, 1] 553 509 c = self.triangles[volume_id, 2] 554 510 555 for vertex_id, node_id in enumerate([a, b,c]):511 for vertex_id, node_id in enumerate([a, b, c]): 556 512 if vertexlist[node_id] is None: 557 513 vertexlist[node_id] = [] 558 559 vertexlist[node_id].append( (volume_id, vertex_id) ) 560 514 vertexlist[node_id].append((volume_id, vertex_id)) 561 515 562 516 # Build inverted triangle index array … … 566 520 for volume_id, vertex_id in vertices: 567 521 vertex_value_indices[k] = 3*volume_id + vertex_id 568 569 522 k += 1 570 523 … … 573 526 self.vertex_value_indices = vertex_value_indices 574 527 575 576 577 578 528 def get_extent(self, absolute=False): 579 529 """Return min and max of all x and y coordinates … … 582 532 are to be made absolute by taking georeference into account 583 533 """ 584 585 586 534 587 535 C = self.get_vertex_coordinates(absolute=absolute) … … 589 537 Y = C[:,1:6:2].copy() 590 538 591 xmin = min(X.flat)592 xmax = max(X.flat)593 ymin = min(Y.flat)594 ymax = max(Y.flat)595 #print "C",C 539 xmin = num.min(X) 540 xmax = num.max(X) 541 ymin = num.min(Y) 542 ymax = num.max(Y) 543 596 544 return xmin, xmax, ymin, ymax 597 545 598 546 def get_areas(self): 599 """Get areas of all individual triangles.600 """ 601 return self.areas 547 """Get areas of all individual triangles.""" 548 549 return self.areas 602 550 603 551 def get_area(self): 604 """Return total area of mesh 605 """ 552 """Return total area of mesh""" 606 553 607 554 return num.sum(self.areas) … … 609 556 def set_georeference(self, g): 610 557 self.geo_reference = g 611 558 612 559 def get_georeference(self): 613 560 return self.geo_reference 614 615 616 561 -
anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
r6928 r7276 8 8 from anuga.fit_interpolate.interpolate import Modeltime_too_early 9 9 10 import Numericas num10 import numpy as num 11 11 12 12 … … 69 69 raise Exception, msg 70 70 71 self.conserved_quantities=num.array(conserved_quantities, num. Float)71 self.conserved_quantities=num.array(conserved_quantities, num.float) 72 72 73 73 def __repr__(self): … … 129 129 130 130 try: 131 q = num.array(q, num. Float)131 q = num.array(q, num.float) 132 132 except: 133 133 msg = 'Return value from time boundary function could ' 134 msg += 'not be converted into a Numeric array of floats.\n'134 msg += 'not be converted into a numeric array of floats.\n' 135 135 msg += 'Specified function should return either list or array.\n' 136 136 msg += 'I got %s' %str(q) … … 241 241 242 242 if verbose: print 'Find midpoint coordinates of entire boundary' 243 self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num. Float)243 self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.float) 244 244 boundary_keys = domain.boundary.keys() 245 245 … … 261 261 262 262 # Compute midpoints 263 if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num. Float)264 if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num. Float)265 if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num. Float)263 if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float) 264 if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float) 265 if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float) 266 266 267 267 # Convert to absolute UTM coordinates … … 371 371 # See http://docs.python.org/lib/warning-filter.html 372 372 self.default_boundary_invoked = True 373 374 375 if res == NAN: 373 374 if num.any(res == NAN): 376 375 x,y=self.midpoint_coordinates[i,:] 377 376 msg = 'NAN value found in file_boundary at ' … … 423 422 424 423 This was added by Nils Goseberg et al in April 2009 425 426 424 """ 427 425 … … 429 427 use_cache=False, verbose=False): 430 428 import time 431 from Numeric import array, zeros, Float432 429 from anuga.config import time_format 433 430 from anuga.abstract_2d_finite_volumes.util import file_function … … 445 442 446 443 if verbose: print 'Find midpoint coordinates of entire boundary' 447 self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)444 self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.float) 448 445 boundary_keys = domain.boundary.keys() 449 450 446 451 447 xllcorner = domain.geo_reference.get_xllcorner() 452 448 yllcorner = domain.geo_reference.get_yllcorner() 453 454 449 455 450 # Make ordering unique #FIXME: should this happen in domain.py? … … 459 454 self.boundary_indices = {} 460 455 for i, (vol_id, edge_id) in enumerate(boundary_keys): 461 462 456 base_index = 3*vol_id 463 457 x0, y0 = V[base_index, :] … … 466 460 467 461 # Compute midpoints 468 if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])469 if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])470 if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])462 if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2]) 463 if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2]) 464 if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2]) 471 465 472 466 # Convert to absolute UTM coordinates … … 483 477 if verbose: print 'Initialise file_function' 484 478 self.F = file_function(filename, domain, 485 interpolation_points=self.midpoint_coordinates,479 interpolation_points=self.midpoint_coordinates, 486 480 time_thinning=time_thinning, 487 481 use_cache=use_cache, … … 497 491 # only part of mesh boundary is actually used with a given 498 492 # file boundary sww file. 499 if hasattr(self.F, 'indices_outside_mesh') and\500 len(self.F.indices_outside_mesh) > 0:493 if (hasattr(self.F, 'indices_outside_mesh') and 494 len(self.F.indices_outside_mesh)) > 0: 501 495 msg = 'WARNING: File_boundary has points outside the mesh ' 502 msg += 'given in %s. ' % filename496 msg += 'given in %s. ' % filename 503 497 msg += 'See warning message issued by Interpolation_function ' 504 498 msg += 'for details (should appear above somewhere if ' … … 510 504 #raise Exception(msg) 511 505 512 513 506 q = self.F(0, point_id=0) 514 507 515 508 d = len(domain.conserved_quantities) 516 msg = 'Values specified in file %s must be ' % filename517 msg += ' a list or an array of length %d' %d509 msg = 'Values specified in file %s must be ' % filename 510 msg += 'a list or an array of length %d' % d 518 511 assert len(q) == d, msg 519 512 … … 525 518 def evaluate(self, vol_id=None, edge_id=None): 526 519 """Return linearly interpolated values based on domain.time 527 at midpoint of segment defined by vol_id and edge_id.520 at midpoint of segment defined by vol_id and edge_id. 528 521 """ 529 q = self.domain.get_conserved_quantities(vol_id, edge = edge_id) 522 523 q = self.domain.get_conserved_quantities(vol_id, edge=edge_id) 530 524 t = self.domain.time 531 525 532 526 if vol_id is not None and edge_id is not None: 533 i = self.boundary_indices[ vol_id, edge_id]534 res = self.F(t, point_id =i)527 i = self.boundary_indices[vol_id, edge_id] 528 res = self.F(t, point_id=i) 535 529 536 530 if res == NAN: 537 x,y =self.midpoint_coordinates[i,:]531 x,y = self.midpoint_coordinates[i,:] 538 532 msg = 'NAN value found in file_boundary at ' 539 msg += 'point id #%d: (%.2f, %.2f).\n' % (i, x, y)540 541 if hasattr(self.F, 'indices_outside_mesh') and\542 len(self.F.indices_outside_mesh) > 0 :533 msg += 'point id #%d: (%.2f, %.2f).\n' % (i, x, y) 534 535 if (hasattr(self.F, 'indices_outside_mesh') and 536 len(self.F.indices_outside_mesh) > 0): 543 537 # Check if NAN point is due it being outside 544 538 # boundary defined in sww file. … … 546 540 if i in self.F.indices_outside_mesh: 547 541 msg += 'This point refers to one outside the ' 548 msg += 'mesh defined by the file %s.\n'\ 549 %self.F.filename 542 msg += 'mesh defined by the file %s.\n' % self.F.filename 550 543 msg += 'Make sure that the file covers ' 551 544 msg += 'the boundary segment it is assigned to ' … … 553 546 else: 554 547 msg += 'This point is inside the mesh defined ' 555 msg += 'the file %s.\n' % self.F.filename548 msg += 'the file %s.\n' % self.F.filename 556 549 msg += 'Check this file for NANs.' 557 550 raise Exception, msg -
anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py
r6717 r7276 3 3 """ 4 4 5 import Numericas num5 import numpy as num 6 6 7 7 … … 94 94 index = Index(n,m) 95 95 96 points = num.zeros((Np, 2), num. Float)96 points = num.zeros((Np, 2), num.float) 97 97 98 98 for i in range(m+1): … … 106 106 107 107 108 elements = num.zeros((Nt, 3), num. Int)108 elements = num.zeros((Nt, 3), num.int) 109 109 boundary = {} 110 110 nt = -1 … … 203 203 204 204 205 206 205 def rectangular_periodic(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)): 207 206 … … 260 259 E = EIndex(n,m) 261 260 262 points = num.zeros( (Np,2), num. Float)261 points = num.zeros( (Np,2), num.float) 263 262 264 263 for i in range(m+1): … … 272 271 273 272 274 elements = num.zeros( (Nt,3), num. Int)273 elements = num.zeros( (Nt,3), num.int) 275 274 boundary = {} 276 275 Idgl = [] … … 340 339 Idgr.extend(Idgl) 341 340 342 Idfl = num.array(Idfl, num. Int)343 Idgr = num.array(Idgr, num. Int)341 Idfl = num.array(Idfl, num.int) 342 Idgr = num.array(Idgr, num.int) 344 343 345 344 full_send_dict[processor] = [Idfl, Idfl] -
anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
r6654 r7276 12 12 from math import pi, sqrt 13 13 14 import Numericas num15 14 import numpy as num 15 16 16 17 17 class Mesh(General_mesh): … … 24 24 coordinate set. 25 25 26 Coordinate sets are implemented as an N x 2 Numeric array containing26 Coordinate sets are implemented as an N x 2 numeric array containing 27 27 x and y coordinates. 28 28 … … 33 33 where 34 34 35 coordinates is either a list of 2-tuples or an Mx2 Numeric array of35 coordinates is either a list of 2-tuples or an Mx2 numeric array of 36 36 floats representing all x, y coordinates in the mesh. 37 37 38 triangles is either a list of 3-tuples or an Nx3 Numeric array of38 triangles is either a list of 3-tuples or an Nx3 numeric array of 39 39 integers representing indices of all vertices in the mesh. 40 40 Each vertex is identified by its index i in [0, M-1]. … … 76 76 """ 77 77 Build triangles from x,y coordinates (sequence of 2-tuples or 78 Mx2 Numeric array of floats) and triangles (sequence of 3-tuples79 or Nx3 Numeric array of non-negative integers).78 Mx2 numeric array of floats) and triangles (sequence of 3-tuples 79 or Nx3 numeric array of non-negative integers). 80 80 """ 81 81 … … 91 91 verbose=verbose) 92 92 93 if verbose: print 'Initialising mesh' 93 if verbose: print 'Initialising mesh' 94 94 95 95 N = len(self) #Number_of_triangles … … 98 98 99 99 #Allocate space for geometric quantities 100 self.centroid_coordinates = num.zeros((N, 2), num. Float)101 102 self.radii = num.zeros(N, num. Float)103 104 self.neighbours = num.zeros((N, 3), num. Int)105 self.neighbour_edges = num.zeros((N, 3), num. Int)106 self.number_of_boundaries = num.zeros(N, num. Int)107 self.surrogate_neighbours = num.zeros((N, 3), num. Int)100 self.centroid_coordinates = num.zeros((N, 2), num.float) 101 102 self.radii = num.zeros(N, num.float) 103 104 self.neighbours = num.zeros((N, 3), num.int) 105 self.neighbour_edges = num.zeros((N, 3), num.int) 106 self.number_of_boundaries = num.zeros(N, num.int) 107 self.surrogate_neighbours = num.zeros((N, 3), num.int) 108 108 109 109 #Get x,y coordinates for all triangles and store … … 111 111 112 112 #Initialise each triangle 113 if verbose: print 'Mesh: Computing centroids and radii' 113 if verbose: print 'Mesh: Computing centroids and radii' 114 114 for i in range(N): 115 115 if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N) … … 117 117 x0, y0 = V[3*i, :] 118 118 x1, y1 = V[3*i+1, :] 119 x2, y2 = V[3*i+2, :] 119 x2, y2 = V[3*i+2, :] 120 120 121 121 #x0 = V[i, 0]; y0 = V[i, 1] … … 124 124 125 125 #Compute centroid 126 centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num. Float)126 centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num.float) 127 127 self.centroid_coordinates[i] = centroid 128 128 … … 133 133 134 134 #Midpoints 135 m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num. Float)136 m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num. Float)137 m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num. Float)135 m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float) 136 m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float) 137 m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float) 138 138 139 139 #The radius is the distance from the centroid of … … 166 166 167 167 #Build neighbour structure 168 if verbose: print 'Mesh: Building neigbour structure' 168 if verbose: print 'Mesh: Building neigbour structure' 169 169 self.build_neighbour_structure() 170 170 … … 178 178 179 179 #Build tagged element dictionary mapping (tag) to array of elements 180 if verbose: print 'Mesh: Building tagged elements dictionary' 180 if verbose: print 'Mesh: Building tagged elements dictionary' 181 181 self.build_tagged_elements_dictionary(tagged_elements) 182 182 … … 184 184 self.lone_vertices = [] 185 185 #Check that all vertices have been registered 186 for node, count in enumerate(self.number_of_triangles_per_node): 186 for node, count in enumerate(self.number_of_triangles_per_node): 187 187 #msg = 'Node %d does not belong to an element.' %node 188 188 #assert count > 0, msg 189 189 if count == 0: 190 190 self.lone_vertices.append(node) 191 191 192 192 #Update boundary indices FIXME: OBSOLETE 193 193 #self.build_boundary_structure() 194 194 195 195 #FIXME check integrity? 196 if verbose: print 'Mesh: Done' 196 if verbose: print 'Mesh: Done' 197 197 198 198 def __repr__(self): … … 394 394 #Check that all keys in given boundary exist 395 395 for tag in tagged_elements.keys(): 396 tagged_elements[tag] = num.array(tagged_elements[tag], num. Int)396 tagged_elements[tag] = num.array(tagged_elements[tag], num.int) 397 397 398 398 msg = 'Not all elements exist. ' … … 400 400 #print "tagged_elements", tagged_elements 401 401 self.tagged_elements = tagged_elements 402 402 403 403 def get_tagged_elements(self): 404 404 return self.tagged_elements … … 459 459 460 460 Using the mesh boundary, derive a bounding polygon for this mesh. 461 If multiple vertex values are present (vertices stored uniquely), 462 461 If multiple vertex values are present (vertices stored uniquely), 462 the algorithm will select the path that contains the entire mesh. 463 463 464 464 All points are in absolute UTM coordinates 465 465 """ 466 467 from anuga.utilities.numerical_tools import angle, ensure_numeric 468 466 467 from anuga.utilities.numerical_tools import angle, ensure_numeric 469 468 470 469 # Get mesh extent 471 470 xmin, xmax, ymin, ymax = self.get_extent(absolute=True) 472 471 pmin = ensure_numeric([xmin, ymin]) 473 pmax = ensure_numeric([xmax, ymax]) 474 472 pmax = ensure_numeric([xmax, ymax]) 475 473 476 474 # Assemble dictionary of boundary segments and choose starting point … … 478 476 inverse_segments = {} 479 477 p0 = None 480 mindist = num.sqrt(num.sum((pmax-pmin)**2)) # Start value across entire mesh 478 479 # Start value across entire mesh 480 mindist = num.sqrt(num.sum((pmax-pmin)**2)) 481 481 for i, edge_id in self.boundary.keys(): 482 482 # Find vertex ids for boundary segment … … 485 485 if edge_id == 2: a = 0; b = 1 486 486 487 A = self.get_vertex_coordinate(i, a, absolute=True) # Start 488 B = self.get_vertex_coordinate(i, b, absolute=True) # End 489 487 A = self.get_vertex_coordinate(i, a, absolute=True) # Start 488 B = self.get_vertex_coordinate(i, b, absolute=True) # End 490 489 491 490 # Take the point closest to pmin as starting point … … 503 502 p0 = B 504 503 505 506 504 # Sanity check 507 505 if p0 is None: 508 raise Exception('Impossible')509 506 msg = 'Impossible: p0 is None!?' 507 raise Exception, msg 510 508 511 509 # Register potential paths from A to B 512 510 if not segments.has_key(tuple(A)): 513 segments[tuple(A)] = [] # Empty list for candidate points 514 515 segments[tuple(A)].append(B) 516 511 segments[tuple(A)] = [] # Empty list for candidate points 512 513 segments[tuple(A)].append(B) 517 514 518 515 # Start with smallest point and follow boundary (counter clock wise) 519 516 polygon = [list(p0)]# Storage for final boundary polygon 520 point_registry = {} # Keep track of storage to avoid multiple runs 521 # around boundary. This will only be the case if 522 517 point_registry = {} # Keep track of storage to avoid multiple runs 518 # around boundary. This will only be the case if 519 # there are more than one candidate. 523 520 # FIXME (Ole): Perhaps we can do away with polygon 524 521 # and use only point_registry to save space. 525 522 526 point_registry[tuple(p0)] = 0 527 523 point_registry[tuple(p0)] = 0 524 528 525 while len(point_registry) < len(self.boundary): 529 530 526 candidate_list = segments[tuple(p0)] 531 527 if len(candidate_list) > 1: 532 # Multiple points detected (this will be the case for meshes 533 # with duplicate points as those used for discontinuous 534 535 # Take the candidate that is furthest to the clockwise 536 537 538 539 528 # Multiple points detected (this will be the case for meshes 529 # with duplicate points as those used for discontinuous 530 # triangles with vertices stored uniquely). 531 # Take the candidate that is furthest to the clockwise 532 # direction, as that will follow the boundary. 533 # 534 # This will also be the case for pathological triangles 535 # that have no neighbours. 540 536 541 537 if verbose: 542 print 'Point %s has multiple candidates: %s'\543 %(str(p0), candidate_list)538 print ('Point %s has multiple candidates: %s' 539 % (str(p0), candidate_list)) 544 540 545 541 # Check that previous are not in candidate list … … 548 544 549 545 # Choose vector against which all angles will be measured 550 if len(polygon) > 1: 551 v_prev = p0 - polygon[-2] # Vector that leads to p0552 546 if len(polygon) > 1: 547 v_prev = p0 - polygon[-2] # Vector that leads to p0 548 # from previous point 553 549 else: 554 # FIXME (Ole): What do we do if the first point has 555 550 # FIXME (Ole): What do we do if the first point has 551 # multiple candidates? 556 552 # Being the lower left corner, perhaps we can use the 557 # vector [1, 0], but I really don't know if this is 558 553 # vector [1, 0], but I really don't know if this is 554 # completely watertight. 559 555 v_prev = [1.0, 0.0] 560 561 562 # Choose candidate with minimum angle 556 557 # Choose candidate with minimum angle 563 558 minimum_angle = 2*pi 564 559 for pc in candidate_list: 565 566 vc = pc-p0 # Candidate vector (from p0 to candidate pt) 567 560 vc = pc-p0 # Candidate vector (from p0 to candidate pt) 561 568 562 # Angle between each candidate and the previous vector 569 563 # in [-pi, pi] 570 564 ac = angle(vc, v_prev) 571 565 if ac > pi: 572 # Give preference to angles on the right hand side 573 # of v_prev 566 # Give preference to angles on the right hand side 567 # of v_prev 574 568 # print 'pc = %s, changing angle from %f to %f'\ 575 569 # %(pc, ac*180/pi, (ac-2*pi)*180/pi) 576 570 ac = ac-2*pi 577 571 578 # Take the minimal angle corresponding to the 579 572 # Take the minimal angle corresponding to the 573 # rightmost vector 580 574 if ac < minimum_angle: 581 575 minimum_angle = ac 582 p1 = pc # Best candidate 583 576 p1 = pc # Best candidate 584 577 585 578 if verbose is True: 586 579 print ' Best candidate %s, angle %f'\ 587 580 %(p1, minimum_angle*180/pi) 588 589 581 else: 590 582 p1 = candidate_list[0] 591 583 592 593 584 if point_registry.has_key(tuple(p1)): 594 # We have reached a point already visited. 595 596 if num.allclose(p1, polygon[0]): 597 # If it is the initial point, the polygon is complete. 598 585 # We have reached a point already visited. 586 if num.allclose(p1, polygon[0]): 587 # If it is the initial point, the polygon is complete. 599 588 if verbose is True: 600 601 print polygon 602 589 print ' Stop criterion fulfilled at point %s' %p1 590 print polygon 591 603 592 # We have completed the boundary polygon - yeehaa 604 605 else: 606 607 # This would be a pathological triangle, but the 608 609 610 593 break 594 else: 595 # The point already visited is not the initial point 596 # This would be a pathological triangle, but the 597 # algorithm must be able to deal with this 598 pass 599 611 600 else: 612 601 # We are still finding new points on the boundary 613 602 point_registry[tuple(p1)] = len(point_registry) 614 615 polygon.append(list(p1)) # De-Numeric each point :-)603 604 polygon.append(list(p1)) # De-numeric each point :-) 616 605 p0 = p1 617 606 618 619 607 return polygon 620 621 608 622 609 def check_integrity(self): … … 641 628 x1, y1 = V[3*i+1, :] 642 629 x2, y2 = V[3*i+2, :] 643 630 644 631 # Check that area hasn't been compromised 645 632 area = self.areas[i] … … 678 665 679 666 x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product 680 667 681 668 msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v) 682 669 msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,)) 683 msg += ' Inner product is %e.' %x 670 msg += ' Inner product is %e.' %x 684 671 assert x < epsilon, msg 685 672 … … 690 677 for i in range(N): 691 678 # For each triangle 692 679 693 680 for k, neighbour_id in enumerate(self.neighbours[i,:]): 694 681 … … 752 739 # Node is lone - i.e. not part of the mesh 753 740 continue 754 741 755 742 k += 1 756 743 757 744 volume_id = index / 3 758 745 vertex_id = index % 3 759 746 760 747 msg = 'Triangle %d, vertex %d points to %d. Should have been %d'\ 761 748 %(volume_id, vertex_id, self.triangles[volume_id, vertex_id], current_node) 762 749 assert self.triangles[volume_id, vertex_id] == current_node, msg 763 750 764 751 if self.number_of_triangles_per_node[current_node] == k: 765 752 # Move on to next node … … 788 775 if not self.geo_reference.is_absolute(): 789 776 V = self.geo_reference.get_absolute(V) 790 777 791 778 return V 792 779 793 780 794 781 def get_radii(self): 795 782 """Return all radii. 796 783 Return radius of inscribed cirle for all triangles 797 784 """ 798 return self.radii 785 return self.radii 799 786 800 787 … … 826 813 str += ' Areas [m^2]:\n' 827 814 str += ' A in [%f, %f]\n' %(min(areas), max(areas)) 828 str += ' number of distinct areas: %d\n' %(len(areas)) 815 str += ' number of distinct areas: %d\n' %(len(areas)) 829 816 str += ' Histogram:\n' 830 817 … … 833 820 lo = hi 834 821 if i+1 < len(bins): 835 #Open upper interval 822 #Open upper interval 836 823 hi = bins[i+1] 837 str += ' [%f, %f[: %d\n' %(lo, hi, count) 824 str += ' [%f, %f[: %d\n' %(lo, hi, count) 838 825 else: 839 826 #Closed upper interval … … 849 836 k = 0 850 837 lower = min(areas) 851 for i, a in enumerate(areas): 852 if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas 838 for i, a in enumerate(areas): 839 if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas 853 840 str += ' %d triangles in [%f, %f]\n' %(i-k, lower, a) 854 841 lower = a 855 842 k = i 856 843 857 844 str += ' %d triangles in [%f, %f]\n'\ 858 %(N-k, lower, max(areas)) 859 860 845 %(N-k, lower, max(areas)) 846 847 861 848 str += ' Boundary:\n' 862 849 str += ' Number of boundary segments == %d\n' %(len(self.boundary)) 863 str += ' Boundary tags == %s\n' %self.get_boundary_tags() 850 str += ' Boundary tags == %s\n' %self.get_boundary_tags() 864 851 str += '------------------------------------------------\n' 865 852 866 853 867 854 return str 868 855 869 856 870 857 def get_triangle_containing_point(self, point): … … 874 861 875 862 """ 876 863 877 864 # FIXME(Ole): This function is currently brute force 878 865 # because I needed it for diagnostics. … … 884 871 polygon = self.get_boundary_polygon() 885 872 #print polygon, point 886 873 887 874 if is_outside_polygon(point, polygon): 888 875 msg = 'Point %s is outside mesh' %str(point) … … 899 886 if is_inside_polygon(point, poly, closed=True): 900 887 return i 901 888 902 889 return 903 890 … … 930 917 Resulting segments are unsorted 931 918 """ 932 919 933 920 V = self.get_vertex_coordinates() 934 921 N = len(self) 935 922 936 923 # Adjust polyline to mesh spatial origin 937 924 polyline = self.geo_reference.get_relative(polyline) … … 939 926 if use_cache is True: 940 927 segments = cache(get_intersecting_segments, 941 (V, N, polyline), 928 (V, N, polyline), 942 929 {'verbose': verbose}, 943 930 verbose=verbose) 944 else: 931 else: 945 932 segments = get_intersecting_segments(V, N, polyline, 946 933 verbose=verbose) 947 934 948 935 949 936 return segments 950 951 937 938 952 939 953 940 def get_triangle_neighbours(self, tri_id): … … 956 943 957 944 Negative returned triangle id's represent a boundary as a neighbour. 958 945 959 946 If the given triangle id is bad, return an empty list. 960 947 """ … … 964 951 except IndexError: 965 952 return [] 966 953 967 954 968 955 def get_interpolation_object(self): 969 956 """Get object I that will allow linear interpolation using this mesh 970 971 This is a time consuming process but it needs only to be 957 958 This is a time consuming process but it needs only to be 972 959 once for the mesh. 973 974 Interpolation can then be done using 975 976 result = I.interpolate_block(vertex_values, interpolation_points) 977 960 961 Interpolation can then be done using 962 963 result = I.interpolate_block(vertex_values, interpolation_points) 964 978 965 where vertex values have been obtained from a quantity using 979 966 vertex_values, triangles = self.get_vertex_values() … … 984 971 else: 985 972 from anuga.fit_interpolate.interpolate import Interpolate 986 987 # Get discontinuous mesh - this will match internal 973 974 # Get discontinuous mesh - this will match internal 988 975 # representation of vertex values 989 976 triangles = self.get_disconnected_triangles() … … 992 979 I = Interpolate(vertex_coordinates, triangles) 993 980 self.interpolation_object = I 994 995 return I 996 981 982 return I 983 997 984 998 985 class Triangle_intersection: 999 986 """Store information about line segments intersecting a triangle 1000 987 1001 988 Attributes are 1002 989 … … 1014 1001 length=None, 1015 1002 triangle_id=None): 1016 self.segment = segment 1003 self.segment = segment 1017 1004 self.normal = normal 1018 1005 self.length = length 1019 1006 self.triangle_id = triangle_id 1020 1007 1021 1008 1022 1009 def __repr__(self): … … 1027 1014 self.length, 1028 1015 self.triangle_id) 1029 1016 1030 1017 return s 1031 1018 1032 1019 1033 1020 1034 1021 def _get_intersecting_segments(V, N, line, … … 1051 1038 from anuga.utilities.polygon import intersection 1052 1039 from anuga.utilities.polygon import is_inside_polygon 1053 1040 1054 1041 msg = 'Line segment must contain exactly two points' 1055 1042 assert len(line) == 2, msg … … 1060 1047 eta0 = line[0][1] 1061 1048 1062 1049 1063 1050 # Check intersection with edge segments for all triangles 1064 1051 # FIXME (Ole): This should be implemented in C … … 1069 1056 x1, y1 = V[3*i+1, :] 1070 1057 x2, y2 = V[3*i+2, :] 1071 1058 1072 1059 1073 1060 edge_segments = [[[x0,y0], [x1, y1]], … … 1076 1063 1077 1064 # Find segments that are intersected by line 1078 1065 1079 1066 intersections = {} # Use dictionary to record points only once 1080 1067 for edge in edge_segments: … … 1082 1069 status, value = intersection(line, edge) 1083 1070 #if value is not None: print 'Triangle %d, Got' %i, status, value 1084 1071 1085 1072 if status == 1: 1086 1073 # Normal intersection of one edge or vertex 1087 intersections[tuple(value)] = i 1074 intersections[tuple(value)] = i 1088 1075 1089 1076 # Exclude singular intersections with vertices … … 1101 1088 # along this edge, it will be recorded here. 1102 1089 intersections[tuple(value[0,:])] = i 1103 intersections[tuple(value[1,:])] = i 1104 1105 1090 intersections[tuple(value[1,:])] = i 1091 1092 1106 1093 if len(intersections) == 1: 1107 1094 # Check if either line end point lies fully within this triangle … … 1113 1100 intersections[tuple(line[1])] = i 1114 1101 elif is_inside_polygon(line[0], poly, closed=False): 1115 intersections[tuple(line[0])] = i 1102 intersections[tuple(line[0])] = i 1116 1103 else: 1117 # Ignore situations where one vertex is touch, for instance 1104 # Ignore situations where one vertex is touch, for instance 1118 1105 continue 1119 1106 … … 1139 1126 1140 1127 # Distances from line origin to the two intersections 1141 z0 = num.array([x0 - xi0, y0 - eta0], num. Float)1142 z1 = num.array([x1 - xi0, y1 - eta0], num. Float)1143 d0 = num.sqrt(num.sum(z0**2)) 1128 z0 = num.array([x0 - xi0, y0 - eta0], num.float) 1129 z1 = num.array([x1 - xi0, y1 - eta0], num.float) 1130 d0 = num.sqrt(num.sum(z0**2)) 1144 1131 d1 = num.sqrt(num.sum(z1**2)) 1145 1132 1146 1133 if d1 < d0: 1147 1134 # Swap … … 1151 1138 1152 1139 # (x0,y0) is now the origin of the intersecting segment 1153 1140 1154 1141 1155 1142 # Normal direction: 1156 1143 # Right hand side relative to line direction 1157 vector = num.array([x1 - x0, y1 - y0], num. Float) # Segment vector1144 vector = num.array([x1 - x0, y1 - y0], num.float) # Segment vector 1158 1145 length = num.sqrt(num.sum(vector**2)) # Segment length 1159 normal = num.array([vector[1], -vector[0]], num. Float)/length1160 1161 1162 segment = ((x0,y0), (x1, y1)) 1146 normal = num.array([vector[1], -vector[0]], num.float)/length 1147 1148 1149 segment = ((x0,y0), (x1, y1)) 1163 1150 T = Triangle_intersection(segment=segment, 1164 1151 normal=normal, … … 1168 1155 1169 1156 # Add segment unless it was done earlier 1170 if not triangle_intersections.has_key(segment): 1157 if not triangle_intersections.has_key(segment): 1171 1158 triangle_intersections[segment] = T 1172 1159 1173 1160 1174 # Return segments as a list 1161 # Return segments as a list 1175 1162 return triangle_intersections.values() 1176 1163 1177 1164 1178 1165 def get_intersecting_segments(V, N, polyline, 1179 verbose=False): 1166 verbose=False): 1180 1167 """Internal function to find edges intersected by Polyline 1181 1168 1182 1169 Input: 1183 1170 V: Vertex coordinates as obtained by mesh.get_vertex_coordinates() 1184 1171 N: Number of triangles in mesh 1185 polyline - list of points forming a segmented line 1172 polyline - list of points forming a segmented line 1186 1173 verbose 1187 1174 Output: … … 1190 1177 This method is used by the public method 1191 1178 get_intersecting_segments(self, polyline) which also contains 1192 more documentation. 1179 more documentation. 1193 1180 """ 1194 1181 1195 1182 msg = 'Polyline must contain at least two points' 1196 1183 assert len(polyline) >= 2, msg 1197 1198 1184 1185 1199 1186 # For all segments in polyline 1200 1187 triangle_intersections = [] … … 1206 1193 print '(%.2f, %.2f) - (%.2f, %.2f)' %(point0[0], point0[1], 1207 1194 point1[0], point1[1]) 1208 1195 1209 1196 line = [point0, point1] 1210 1197 triangle_intersections += _get_intersecting_segments(V, N, line, … … 1214 1201 msg = 'No segments found' 1215 1202 assert len(triangle_intersections) > 0, msg 1216 1217 1203 1204 1218 1205 return triangle_intersections 1219 1220 1221 1222 1223 1206 1207 1208 1209 1210 1224 1211 def segment_midpoints(segments): 1225 1212 """Calculate midpoints of all segments 1226 1213 1227 1214 Inputs: 1228 1215 segments: List of instances of class Segment 1229 1216 1230 1217 Ouputs: 1231 1218 midpoints: List of points 1232 1219 """ 1233 1220 1234 1221 midpoints = [] 1235 1222 msg = 'Elements of input list to segment_midpoints must be of class Triangle_intersection' 1236 1223 for segment in segments: 1237 1224 assert isinstance(segment, Triangle_intersection), msg 1238 1239 midpoint = num.sum(num.array(segment.segment, num. Float))/21225 1226 midpoint = num.sum(num.array(segment.segment, num.float), axis=0)/2 1240 1227 midpoints.append(midpoint) 1241 1228 1242 1229 return midpoints 1243 1244 1245 1230 1231 1232 -
anuga_core/source/anuga/abstract_2d_finite_volumes/pmesh2domain.py
r6688 r7276 1 """Class pmesh2domain - Converting .tsh files to do amains1 """Class pmesh2domain - Converting .tsh files to domains 2 2 3 3 … … 8 8 9 9 import sys 10 11 import Numeric as num 12 13 14 def pmesh_instance_to_domain_instance(mesh, 15 DomainClass): 16 """ 17 Convert a pmesh instance/object into a domain instance. 18 19 Use pmesh_to_domain_instance to convert a mesh file to a domain instance. 20 """ 21 22 vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \ 23 ,tagged_elements_dict, geo_reference = \ 24 pmesh_to_domain(mesh_instance=mesh) 10 import numpy as num 11 12 13 ## 14 # @brief Convert a pmesh instance to a domain instance. 15 # @param mesh The pmesh instance to convert. 16 # @param DomainClass The class to instantiate and return. 17 # @return The converted pmesh instance (as a 'DomainClass' instance). 18 def pmesh_instance_to_domain_instance(mesh, DomainClass): 19 """Convert a pmesh instance/object into a domain instance. 20 21 Uses pmesh_to_domain_instance to convert a mesh file to a domain instance. 22 """ 23 24 (vertex_coordinates, vertices, tag_dict, vertex_quantity_dict, 25 tagged_elements_dict, geo_reference) = pmesh_to_domain(mesh_instance=mesh) 26 27 # NOTE(Ole): This import cannot be at the module level 28 # due to mutual dependency with domain.py 29 from anuga.abstract_2d_finite_volumes.domain import Domain 30 31 # ensure that the required 'DomainClass' actually is an instance of Domain 32 msg = ('The class %s is not a subclass of the generic domain class %s' 33 % (DomainClass, Domain)) 34 assert issubclass(DomainClass, Domain), msg 35 36 # instantiate the result class 37 result = DomainClass(coordinates=vertex_coordinates, 38 vertices=vertices, 39 boundary=tag_dict, 40 tagged_elements=tagged_elements_dict, 41 geo_reference=geo_reference) 42 43 # set the water stage to be the elevation 44 if (vertex_quantity_dict.has_key('elevation') and 45 not vertex_quantity_dict.has_key('stage')): 46 vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation'] 47 result.set_quantity_vertices_dict(vertex_quantity_dict) 48 49 return result 50 51 52 ## 53 # @brief Convert a mesh file to a Domain instance. 54 # @param file_name Name of the file to convert (TSH or MSH). 55 # @param DomainClass Class of return instance. 56 # @param use_cache True if caching is to be used. 57 # @param verbose True if this function is to be verbose. 58 # @return An instance of 'DomainClass' containing the file data. 59 def pmesh_to_domain_instance(file_name, DomainClass, use_cache=False, 60 verbose=False): 61 """Converts a mesh file(.tsh or .msh), to a Domain instance. 62 63 file_name is the name of the mesh file to convert, including the extension 64 65 DomainClass is the Class that will be returned. 66 It must be a subclass of Domain, with the same interface as domain. 67 68 use_cache: True means that caching is attempted for the computed domain. 69 """ 70 71 if use_cache is True: 72 from caching import cache 73 result = cache(_pmesh_to_domain_instance, (file_name, DomainClass), 74 dependencies=[file_name], verbose=verbose) 75 else: 76 result = apply(_pmesh_to_domain_instance, (file_name, DomainClass)) 77 78 return result 79 80 81 ## 82 # @brief Convert a mesh file to a Domain instance. 83 # @param file_name Name of the file to convert (TSH or MSH). 84 # @param DomainClass Class of return instance. 85 # @return The DomainClass instance containing the file data. 86 def _pmesh_to_domain_instance(file_name, DomainClass): 87 """Converts a mesh file(.tsh or .msh), to a Domain instance. 88 89 Internal function. See public interface pmesh_to_domain_instance for details 90 """ 91 92 (vertex_coordinates, vertices, tag_dict, vertex_quantity_dict, 93 tagged_elements_dict, geo_reference) = pmesh_to_domain(file_name=file_name) 25 94 26 95 # NOTE(Ole): This import cannot be at the module level due to mutual … … 28 97 from anuga.abstract_2d_finite_volumes.domain import Domain 29 98 30 31 32 33 msg = 'The class %s is not a subclass of the generic domain class %s'\ 34 %(DomainClass, Domain) 99 # ensure the required class is a subclass of Domain 100 msg = ('The class %s is not a subclass of the generic domain class %s' 101 % (DomainClass, Domain)) 35 102 assert issubclass(DomainClass, Domain), msg 36 37 103 38 104 domain = DomainClass(coordinates = vertex_coordinates, … … 42 108 geo_reference = geo_reference ) 43 109 44 # set the water stage to be the elevation 45 if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'): 46 vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation'] 47 48 domain.set_quantity_vertices_dict(vertex_quantity_dict) 49 #print "vertex_quantity_dict",vertex_quantity_dict 50 return domain 51 52 53 54 def pmesh_to_domain_instance(file_name, DomainClass, use_cache = False, verbose = False): 55 """ 56 Converts a mesh file(.tsh or .msh), to a Domain instance. 57 58 file_name is the name of the mesh file to convert, including the extension 59 60 DomainClass is the Class that will be returned. 61 It must be a subclass of Domain, with the same interface as domain. 62 63 use_cache: True means that caching is attempted for the computed domain. 64 """ 65 66 if use_cache is True: 67 from caching import cache 68 result = cache(_pmesh_to_domain_instance, (file_name, DomainClass), 69 dependencies = [file_name], 70 verbose = verbose) 71 72 else: 73 result = apply(_pmesh_to_domain_instance, (file_name, DomainClass)) 74 75 return result 76 77 78 79 80 def _pmesh_to_domain_instance(file_name, DomainClass): 81 """ 82 Converts a mesh file(.tsh or .msh), to a Domain instance. 83 84 Internal function. See public interface pmesh_to_domain_instance for details 85 """ 86 87 vertex_coordinates, vertices, tag_dict, vertex_quantity_dict, \ 88 tagged_elements_dict, geo_reference = \ 89 pmesh_to_domain(file_name=file_name) 90 91 92 # NOTE(Ole): This import cannot be at the module level due to mutual 93 # dependency with domain.py 94 from anuga.abstract_2d_finite_volumes.domain import Domain 95 96 97 msg = 'The class %s is not a subclass of the generic domain class %s'\ 98 %(DomainClass, Domain) 99 assert issubclass(DomainClass, Domain), msg 100 101 102 103 domain = DomainClass(coordinates = vertex_coordinates, 104 vertices = vertices, 105 boundary = tag_dict, 106 tagged_elements = tagged_elements_dict, 107 geo_reference = geo_reference ) 108 109 110 111 #FIXME (Ole): Is this really the right place to apply the a default 112 #value specific to the shallow water wave equation? 113 #The 'assert' above indicates that any subclass of Domain is acceptable. 114 #Suggestion - module shallow_water.py will eventually take care of this 115 #(when I get around to it) so it should be removed from here. 110 # FIXME (Ole): Is this really the right place to apply a default 111 # value specific to the shallow water wave equation? 112 # The 'assert' above indicates that any subclass of Domain is acceptable. 113 # Suggestion - module shallow_water.py will eventually take care of this 114 # (when I get around to it) so it should be removed from here. 116 115 117 116 # This doesn't work on the domain instance. 118 117 # This is still needed so -ve elevations don't cuase 'lakes' 119 118 # The fixme we discussed was to only create a quantity when its values 120 # are set.119 # are set. 121 120 # I think that's the way to go still 122 121 123 122 # set the water stage to be the elevation 124 if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'): 123 if (vertex_quantity_dict.has_key('elevation') and 124 not vertex_quantity_dict.has_key('stage')): 125 125 vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation'] 126 127 126 domain.set_quantity_vertices_dict(vertex_quantity_dict) 128 #print "vertex_quantity_dict",vertex_quantity_dict 127 129 128 return domain 130 129 131 130 132 def pmesh_to_domain(file_name=None, 133 mesh_instance=None, 134 use_cache=False, 131 ## 132 # @brief Convert pmesh file/instance to list(s) that can instantiate a Domain. 133 # @param file_name Path to file to convert. 134 # @param mesh_instance Instance to convert. 135 # @param use_cache True if we are to cache. 136 # @param verbose True if this function is to be verbose. 137 # @return ?? 138 def pmesh_to_domain(file_name=None, mesh_instance=None, use_cache=False, 135 139 verbose=False): 136 """ 137 Convert a pmesh file or a pmesh mesh instance to a bunch of lists 140 """Convert a pmesh file or a pmesh mesh instance to a bunch of lists 138 141 that can be used to instanciate a domain object. 139 142 … … 144 147 from caching import cache 145 148 result = cache(_pmesh_to_domain, (file_name, mesh_instance), 146 dependencies = [file_name], 147 verbose = verbose) 149 dependencies=[file_name], verbose=verbose) 148 150 149 151 else: … … 153 155 154 156 155 def _pmesh_to_domain(file_name=None, 156 mesh_instance=None, 157 use_cache=False, 157 ## 158 # @brief Convert pmesh file/instance to list(s) that can instantiate a Domain. 159 # @param file_name Path to file to convert. 160 # @param mesh_instance Instance to convert. 161 # @param use_cache True if we are to cache. 162 # @param verbose True if this function is to be verbose. 163 # @return ?? 164 def _pmesh_to_domain(file_name=None, mesh_instance=None, use_cache=False, 158 165 verbose=False): 159 """ 160 Convert a pmesh file or a pmesh mesh instance to a bunch of lists 161 that can be used to instanciate a domain object. 166 """Convert a pmesh file or a pmesh mesh instance to a bunch of lists 167 that can be used to instantiate a domain object. 162 168 """ 163 169 164 170 from load_mesh.loadASCII import import_mesh_file 165 171 172 # get data from mesh instance or file 166 173 if file_name is None: 167 174 mesh_dict = mesh_instance.Mesh2IODict() 168 175 else: 169 176 mesh_dict = import_mesh_file(file_name) 170 #print "mesh_dict",mesh_dict 177 178 # extract required data from the mesh dictionary 171 179 vertex_coordinates = mesh_dict['vertices'] 172 180 volumes = mesh_dict['triangles'] 173 181 vertex_quantity_dict = {} 182 183 # num.transpose(None) gives scalar array of value None 174 184 point_atts = mesh_dict['vertex_attributes'] 175 point_titles = mesh_dict['vertex_attribute_titles'] 176 geo_reference = mesh_dict['geo_reference'] 185 186 point_titles = mesh_dict['vertex_attribute_titles'] 187 geo_reference = mesh_dict['geo_reference'] 177 188 if point_atts is not None: 178 point_atts = num.transpose(point_atts) 179 for quantity, value_vector in map 189 point_atts = num.transpose(point_atts) 190 for quantity, value_vector in map(None, point_titles, point_atts): 180 191 vertex_quantity_dict[quantity] = value_vector 181 192 tag_dict = pmesh_dict_to_tag_dict(mesh_dict) 182 193 tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict) 183 return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference 184 194 195 return (vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, 196 tagged_elements_dict, geo_reference) 185 197 186 198 187 199 def build_tagged_elements_dictionary(mesh_dict): 188 200 """Build the dictionary of element tags. 201 189 202 tagged_elements is a dictionary of element arrays, 190 keyed by tag: 191 { (tag): [e1, e2, e3..] }192 """ 203 keyed by tag: { (tag): [e1, e2, e3..] } 204 """ 205 193 206 tri_atts = mesh_dict['triangle_tags'] 194 207 tagged_elements = {} … … 199 212 tagged_elements.setdefault(tri_atts[tri_att_index], 200 213 []).append(tri_att_index) 214 201 215 return tagged_elements 216 202 217 203 218 def pmesh_dict_to_tag_dict(mesh_dict): … … 205 220 to a dictionary of tags, indexed with volume id and face number. 206 221 """ 222 207 223 triangles = mesh_dict['triangles'] 208 224 sides = calc_sides(triangles) 209 225 tag_dict = {} 210 for seg, tag in map(None, mesh_dict['segments'], 211 mesh_dict['segment_tags']): 226 for seg, tag in map(None, mesh_dict['segments'], mesh_dict['segment_tags']): 212 227 v1 = int(seg[0]) 213 228 v2 = int(seg[1]) … … 223 238 224 239 def calc_sides(triangles): 225 #Build dictionary mapping from sides (2-tuple of points) 226 #to left hand side neighbouring triangle 240 '''Build dictionary mapping from sides (2-tuple of points) 241 to left hand side neighbouring triangle 242 ''' 243 227 244 sides = {} 245 228 246 for id, triangle in enumerate(triangles): 229 247 a = int(triangle[0]) 230 248 b = int(triangle[1]) 231 249 c = int(triangle[2]) 250 232 251 sides[a,b] = (id, 2) #(id, face) 233 252 sides[b,c] = (id, 0) #(id, face) 234 253 sides[c,a] = (id, 1) #(id, face) 254 235 255 return sides 256 -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
r7044 r7276 14 14 Otherwise raise an exception 15 15 """ 16 17 from types import FloatType, IntType, LongType, NoneType 16 18 17 19 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar … … 23 25 from anuga.caching import cache 24 26 25 import Numeric as num 27 import anuga.utilities.numerical_tools as aunt 28 29 import numpy as num 26 30 27 31 … … 43 47 if vertex_values is None: 44 48 N = len(domain) # number_of_elements 45 self.vertex_values = num.zeros((N, 3), num. Float)49 self.vertex_values = num.zeros((N, 3), num.float) 46 50 else: 47 self.vertex_values = num.array(vertex_values, num. Float)51 self.vertex_values = num.array(vertex_values, num.float) 48 52 49 53 N, V = self.vertex_values.shape … … 57 61 58 62 # Allocate space for other quantities 59 self.centroid_values = num.zeros(N, num. Float)60 self.edge_values = num.zeros((N, 3), num. Float)63 self.centroid_values = num.zeros(N, num.float) 64 self.edge_values = num.zeros((N, 3), num.float) 61 65 62 66 # Allocate space for Gradient 63 self.x_gradient = num.zeros(N, num. Float)64 self.y_gradient = num.zeros(N, num. Float)67 self.x_gradient = num.zeros(N, num.float) 68 self.y_gradient = num.zeros(N, num.float) 65 69 66 70 # Allocate space for Limiter Phi 67 self.phi = num.zeros(N, num. Float)71 self.phi = num.zeros(N, num.float) 68 72 69 73 # Intialise centroid and edge_values … … 72 76 # Allocate space for boundary values 73 77 L = len(domain.boundary) 74 self.boundary_values = num.zeros(L, num. Float)78 self.boundary_values = num.zeros(L, num.float) 75 79 76 80 # Allocate space for updates of conserved quantities by … … 78 82 79 83 # Allocate space for update fields 80 self.explicit_update = num.zeros(N, num. Float )81 self.semi_implicit_update = num.zeros(N, num. Float )82 self.centroid_backup_values = num.zeros(N, num. Float)84 self.explicit_update = num.zeros(N, num.float ) 85 self.semi_implicit_update = num.zeros(N, num.float ) 86 self.centroid_backup_values = num.zeros(N, num.float) 83 87 84 88 self.set_beta(1.0) … … 316 320 317 321 numeric: 318 Compatible list, Numeric array (see below) or constant.322 Compatible list, numeric array (see below) or constant. 319 323 If callable it will treated as a function (see below) 320 324 If instance of another Quantity it will be treated as such. … … 357 361 358 362 In case of location == 'centroids' the dimension values must 359 be a list of a Numerical array of length N,363 be a list of a numerical array of length N, 360 364 N being the number of elements. 361 365 Otherwise it must be of dimension Nx3 … … 402 406 403 407 from anuga.geospatial_data.geospatial_data import Geospatial_data 404 from types import FloatType, IntType, LongType, ListType, NoneType405 408 406 409 # Treat special case: Polygon situation … … 457 460 raise Exception, msg 458 461 459 msg = 'Indices must be a list or None'460 assert type(indices) in [ListType, NoneType, num.ArrayType], msg462 msg = 'Indices must be a list, array or None' 463 assert isinstance(indices, (NoneType, list, num.ndarray)), msg 461 464 462 465 # Determine which 'set_values_from_...' to use 463 466 if numeric is not None: 464 if type(numeric) in [FloatType, IntType, LongType]: 465 self.set_values_from_constant(numeric, location, 466 indices, verbose) 467 elif type(numeric) in [num.ArrayType, ListType]: 467 if isinstance(numeric, (list, num.ndarray)): 468 468 self.set_values_from_array(numeric, location, indices, 469 469 use_cache=use_cache, verbose=verbose) … … 479 479 indices, verbose=verbose, 480 480 use_cache=use_cache) 481 else: 482 msg = 'Illegal type for argument numeric: %s' % str(numeric) 483 raise msg 481 else: # see if it's coercible to a float (float, int or long, etc) 482 try: 483 numeric = float(numeric) 484 except ValueError: 485 msg = ("Illegal type for variable 'numeric': %s" 486 % type(numeric)) 487 raise Exception, msg 488 self.set_values_from_constant(numeric, location, 489 indices, verbose) 484 490 elif quantity is not None: 485 491 self.set_values_from_quantity(quantity, location, indices, verbose) … … 515 521 self.extrapolate_first_order() 516 522 517 518 519 523 ############################################################################ 520 524 # Specific internal functions for setting values based on type … … 583 587 """Set values for quantity 584 588 585 values: Numeric array589 values: numeric array 586 590 location: Where values are to be stored. 587 591 Permissible options are: vertices, centroid, unique vertices … … 593 597 594 598 In case of location == 'centroid' the dimension values must 595 be a list of a Numerical array of length N, N being the number599 be a list of a numerical array of length N, N being the number 596 600 of elements. 597 601 … … 607 611 """ 608 612 609 values = num.array(values, num. Float)613 values = num.array(values, num.float) 610 614 611 615 if indices is not None: 612 indices = num.array(indices, num. Int)616 indices = num.array(indices, num.int) 613 617 msg = ('Number of values must match number of indices: You ' 614 618 'specified %d values and %d indices' … … 637 641 'Values array must be 1d') 638 642 639 self.set_vertex_values(values.flat , indices=indices,643 self.set_vertex_values(values.flatten(), indices=indices, 640 644 use_cache=use_cache, verbose=verbose) 641 645 else: … … 690 694 # @param use_cache ?? 691 695 # @param verbose True if this method is to be verbose. 692 def set_values_from_function(self, f, 693 location='vertices', 694 indices=None, 695 use_cache=False, 696 verbose=False): 696 def set_values_from_function(self, 697 f, 698 location='vertices', 699 indices=None, 700 use_cache=False, 701 verbose=False): 697 702 """Set values for quantity using specified function 698 703 … … 716 721 indices = range(len(self)) 717 722 718 V = num.take(self.domain.get_centroid_coordinates(), indices )723 V = num.take(self.domain.get_centroid_coordinates(), indices, axis=0) 719 724 x = V[:,0]; y = V[:,1] 720 725 if use_cache is True: … … 772 777 # @param verbose ?? 773 778 # @param use_cache ?? 774 def set_values_from_geospatial_data(self, geospatial_data, 775 alpha, 776 location, 777 indices, 778 verbose=False, 779 use_cache=False): 779 def set_values_from_geospatial_data(self, 780 geospatial_data, 781 alpha, 782 location, 783 indices, 784 verbose=False, 785 use_cache=False): 780 786 """Set values based on geo referenced geospatial data object.""" 781 787 … … 788 794 from anuga.coordinate_transforms.geo_reference import Geo_reference 789 795 790 points = ensure_numeric(points, num. Float)791 values = ensure_numeric(values, num. Float)796 points = ensure_numeric(points, num.float) 797 values = ensure_numeric(values, num.float) 792 798 793 799 if location != 'vertices': … … 829 835 # @param verbose True if this method is to be verbose. 830 836 # @param use_cache ?? 831 def set_values_from_points(self, points, 832 values, 833 alpha, 834 location, 835 indices, 836 data_georef=None, 837 verbose=False, 838 use_cache=False): 837 def set_values_from_points(self, 838 points, 839 values, 840 alpha, 841 location, 842 indices, 843 data_georef=None, 844 verbose=False, 845 use_cache=False): 839 846 """Set quantity values from arbitray data points using fit_interpolate.fit""" 840 847 … … 851 858 # @param use_cache 852 859 # @param max_read_lines 853 def set_values_from_file(self, filename, 854 attribute_name, 855 alpha, 856 location, 857 indices, 858 verbose=False, 859 use_cache=False, 860 max_read_lines=None): 860 def set_values_from_file(self, 861 filename, 862 attribute_name, 863 alpha, 864 location, 865 indices, 866 verbose=False, 867 use_cache=False, 868 max_read_lines=None): 861 869 """Set quantity based on arbitrary points in a points file using 862 870 attribute_name selects name of attribute present in file. … … 1061 1069 # @param use_cache ?? 1062 1070 # @param verbose True if this method is to be verbose. 1063 def get_interpolated_values(self, interpolation_points, 1064 use_cache=False, 1065 verbose=False): 1071 def get_interpolated_values(self, 1072 interpolation_points, 1073 use_cache=False, 1074 verbose=False): 1066 1075 """Get values at interpolation points 1067 1076 … … 1079 1088 1080 1089 # Ensure that interpolation points is either a list of 1081 # points, Nx2 array, or geospatial and convert to Numeric array1090 # points, Nx2 array, or geospatial and convert to numeric array 1082 1091 if isinstance(interpolation_points, Geospatial_data): 1083 1092 # Ensure interpolation points are in absolute UTM coordinates … … 1111 1120 # @param use_cache 1112 1121 # @param verbose True if this method is to be verbose. 1113 def get_values(self, 1122 def get_values(self, 1114 1123 interpolation_points=None, 1115 1124 location='vertices', … … 1119 1128 """Get values for quantity 1120 1129 1121 Extract values for quantity as a Numeric array.1130 Extract values for quantity as a numeric array. 1122 1131 1123 1132 Inputs: … … 1137 1146 1138 1147 In case of location == 'centroids' the dimension of returned 1139 values will be a list or a Numerical array of length N, N being1148 values will be a list or a numerical array of length N, N being 1140 1149 the number of elements. 1141 1150 … … 1174 1183 'edges', 'unique vertices']: 1175 1184 msg = 'Invalid location: %s' % location 1176 raise msg1185 raise Exception, msg 1177 1186 1178 1187 import types 1179 1188 1180 assert type(indices) in [types.ListType, types.NoneType, 1181 num.ArrayType], \ 1182 'Indices must be a list or None' 1189 msg = "'indices' must be a list, array or None" 1190 assert isinstance(indices, (NoneType, list, num.ndarray)), msg 1183 1191 1184 1192 if location == 'centroids': 1185 1193 if (indices == None): 1186 1194 indices = range(len(self)) 1187 return num.take(self.centroid_values, indices )1195 return num.take(self.centroid_values, indices, axis=0) 1188 1196 elif location == 'edges': 1189 1197 if (indices == None): 1190 1198 indices = range(len(self)) 1191 return num.take(self.edge_values, indices )1199 return num.take(self.edge_values, indices, axis=0) 1192 1200 elif location == 'unique vertices': 1193 1201 if (indices == None): … … 1211 1219 sum += self.vertex_values[triangle_id, vertex_id] 1212 1220 vert_values.append(sum / len(triangles)) 1213 return num.array(vert_values, num. Float)1221 return num.array(vert_values, num.float) 1214 1222 else: 1215 1223 if (indices is None): 1216 1224 indices = range(len(self)) 1217 return num.take(self.vertex_values, indices )1225 return num.take(self.vertex_values, indices, axis=0) 1218 1226 1219 1227 ## … … 1223 1231 # @param use_cache ?? 1224 1232 # @param verbose?? 1225 def set_vertex_values(self, A, 1226 indices=None, 1227 use_cache=False, 1228 verbose=False): 1233 def set_vertex_values(self, 1234 A, 1235 indices=None, 1236 use_cache=False, 1237 verbose=False): 1229 1238 """Set vertex values for all unique vertices based on input array A 1230 1239 which has one entry per unique vertex, i.e. one value for each row in … … 1237 1246 1238 1247 # Check that A can be converted to array and is of appropriate dim 1239 A = ensure_numeric(A, num. Float)1248 A = ensure_numeric(A, num.float) 1240 1249 assert len(A.shape) == 1 1241 1250 … … 1337 1346 1338 1347 if precision is None: 1339 precision = num. Float1348 precision = num.float 1340 1349 1341 1350 if smooth is True: … … 1343 1352 V = self.domain.get_triangles() 1344 1353 N = self.domain.number_of_full_nodes # Ignore ghost nodes if any 1345 A = num.zeros(N, num. Float)1354 A = num.zeros(N, num.float) 1346 1355 points = self.domain.get_nodes() 1347 1356 … … 1361 1370 if current_node == N: 1362 1371 msg = 'Current node exceeding number of nodes (%d) ' % N 1363 raise msg1372 raise Exception, msg 1364 1373 1365 1374 k += 1 … … 1382 1391 V = self.domain.get_disconnected_triangles() 1383 1392 points = self.domain.get_vertex_coordinates() 1384 A = self.vertex_values.flat .astype(precision)1393 A = self.vertex_values.flatten().astype(precision) 1385 1394 1386 1395 # Return -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
r7154 r7276 11 11 12 12 #include "Python.h" 13 #include " Numeric/arrayobject.h"13 #include "numpy/arrayobject.h" 14 14 #include "math.h" 15 15 … … 1028 1028 } 1029 1029 1030 // check that numpy array objects arrays are C contiguous memory 1031 CHECK_C_CONTIG(vertex_value_indices); 1032 CHECK_C_CONTIG(number_of_triangles_per_node); 1033 CHECK_C_CONTIG(vertex_values); 1034 CHECK_C_CONTIG(A); 1035 1030 1036 N = vertex_value_indices -> dimensions[0]; 1031 1037 // printf("Got parameters, N=%d\n", N); -
anuga_core/source/anuga/abstract_2d_finite_volumes/region.py
r6145 r7276 8 8 # FIXME (DSG-DSG) add better comments 9 9 10 import Numericas num10 import numpy as num 11 11 12 12 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py
r6717 r7276 7 7 from anuga.config import epsilon 8 8 9 import Numericas num9 import numpy as num 10 10 11 11 … … 65 65 66 66 67 assert domain.get_conserved_quantities(0, edge=1) == 0.67 assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.) 68 68 69 69 … … 345 345 346 346 347 A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], num. Float)348 B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], num. Float)347 A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], num.float) 348 B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], num.float) 349 349 350 350 #print A … … 549 549 [0,0,9], [-6, 3, 3]]) 550 550 551 assert num.allclose( domain.quantities['stage'].centroid_values,552 [2,5,3,0])551 assert num.allclose( domain.quantities['stage'].centroid_values, 552 [2,5,3,0] ) 553 553 554 554 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], … … 561 561 domain.distribute_to_vertices_and_edges() 562 562 563 # 563 #First order extrapolation 564 564 assert num.allclose( domain.quantities['stage'].vertex_values, 565 565 [[ 2., 2., 2.], … … 614 614 615 615 sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4]) 616 denom = num.ones(4, num. Float)-domain.timestep*sem616 denom = num.ones(4, num.float) - domain.timestep*sem 617 617 618 618 # x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) … … 669 669 domain.distribute_to_vertices_and_edges() 670 670 671 # 671 #First order extrapolation 672 672 assert num.allclose( domain.quantities['stage'].vertex_values, 673 673 [[ 2., 2., 2.], … … 766 766 767 767 #print domain.quantities['friction'].get_values() 768 msg = ("domain.quantities['friction'].get_values()=\n%s\n" 769 'should equal\n' 770 '[[ 0.09, 0.09, 0.09],\n' 771 ' [ 0.09, 0.09, 0.09],\n' 772 ' [ 0.07, 0.07, 0.07],\n' 773 ' [ 0.07, 0.07, 0.07],\n' 774 ' [ 1.0, 1.0, 1.0],\n' 775 ' [ 1.0, 1.0, 1.0]]' 776 % str(domain.quantities['friction'].get_values())) 768 777 assert num.allclose(domain.quantities['friction'].get_values(), 769 778 [[ 0.09, 0.09, 0.09], … … 772 781 [ 0.07, 0.07, 0.07], 773 782 [ 1.0, 1.0, 1.0], 774 [ 1.0, 1.0, 1.0]]) 783 [ 1.0, 1.0, 1.0]]), msg 775 784 776 785 domain.set_region([set_bottom_friction, set_top_friction]) … … 794 803 [ 11.0, 11.0, 11.0]]) 795 804 796 797 805 def test_rectangular_periodic_and_ghosts(self): 798 806 … … 909 917 910 918 #------------------------------------------------------------- 919 911 920 if __name__ == "__main__": 912 921 suite = unittest.makeSuite(Test_Domain,'test') -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_ermapper.py
r6145 r7276 7 7 from os import remove 8 8 9 import Numericas num9 import numpy as num 10 10 11 11 … … 45 45 46 46 # Write test data 47 ermapper_grids.write_ermapper_data(original_grid, filename, num. Float64)47 ermapper_grids.write_ermapper_data(original_grid, filename, num.float64) 48 48 49 49 # Read in the test data 50 new_grid = ermapper_grids.read_ermapper_data(filename, num. Float64)50 new_grid = ermapper_grids.read_ermapper_data(filename, num.float64) 51 51 52 52 # Check that the test data that has been read in matches the original data … … 183 183 184 184 #------------------------------------------------------------- 185 185 186 if __name__ == "__main__": 186 187 suite = unittest.makeSuite(Test_ERMapper,'test') 187 188 runner = unittest.TextTestRunner() 188 189 runner.run(suite) 189 190 191 192 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py
r6717 r7276 10 10 from anuga.coordinate_transforms.geo_reference import Geo_reference 11 11 12 import Numericas num12 import numpy as num 13 13 14 14 … … 53 53 e = [2.0, 2.0] 54 54 f = [4.0, 0.0] 55 56 55 nodes = num.array([a, b, c, d, e, f]) 57 56 58 57 nodes_absolute = geo.get_absolute(nodes) 59 58 60 #bac, bce, ecf, dbe, daf, dae 61 triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.Int) 62 63 domain = General_mesh(nodes, triangles, 64 geo_reference = geo) 65 verts = domain.get_vertex_coordinates(triangle_id=0) 66 self.assert_(num.allclose(num.array([b,a,c]), verts)) 59 # bac, bce, ecf, dbe 60 triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int) 61 62 domain = General_mesh(nodes, triangles, geo_reference=geo) 63 64 verts = domain.get_vertex_coordinates(triangle_id=0) # bac 65 msg = ("num.array([b,a,c])=\n%s\nshould be close to 'verts'=\n%s" 66 % (str(num.array([b,a,c])), str(verts))) 67 self.failUnless(num.allclose(num.array([b,a,c]), verts), msg) 68 67 69 verts = domain.get_vertex_coordinates(triangle_id=0) 68 self.assert_(num.allclose(num.array([b,a,c]), verts)) 70 msg = ("num.array([b,a,c])=\n%s\nshould be close to 'verts'=\n%s" 71 % (str(num.array([b,a,c])), str(verts))) 72 self.assert_(num.allclose(num.array([b,a,c]), verts), msg) 73 74 verts = domain.get_vertex_coordinates(triangle_id=0, absolute=True) 75 msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s" 76 % (str(num.array([nodes_absolute[1], 77 nodes_absolute[0], 78 nodes_absolute[2]])), 79 str(verts))) 80 self.assert_(num.allclose(num.array([nodes_absolute[1], 81 nodes_absolute[0], 82 nodes_absolute[2]]), 83 verts), msg) 84 69 85 verts = domain.get_vertex_coordinates(triangle_id=0, 70 86 absolute=True) 87 msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s" 88 % (str(num.array([nodes_absolute[1], 89 nodes_absolute[0], 90 nodes_absolute[2]])), 91 str(verts))) 71 92 self.assert_(num.allclose(num.array([nodes_absolute[1], 72 nodes_absolute[0], 73 nodes_absolute[2]]), verts)) 74 verts = domain.get_vertex_coordinates(triangle_id=0, 75 absolute=True) 76 self.assert_(num.allclose(num.array([nodes_absolute[1], 77 nodes_absolute[0], 78 nodes_absolute[2]]), verts)) 79 80 93 nodes_absolute[0], 94 nodes_absolute[2]]), 95 verts), msg) 81 96 82 97 def test_get_vertex_coordinates_triangle_id(self): … … 117 132 domain = General_mesh(nodes, triangles) 118 133 119 value = [7] 120 assert num.allclose(domain.get_triangles(), triangles) 134 msg = ("domain.get_triangles()=\n%s\nshould be the same as " 135 "'triangles'=\n%s" 136 % (str(domain.get_triangles()), str(triangles))) 137 assert num.allclose(domain.get_triangles(), triangles), msg 138 msg = ("domain.get_triangles([0,4])=\n%s\nshould be the same as " 139 "'[triangles[0], triangles[4]]' which is\n%s" 140 % (str(domain.get_triangles([0,4])), 141 str([triangles[0], triangles[4]]))) 121 142 assert num.allclose(domain.get_triangles([0,4]), 122 [triangles[0], triangles[4]])143 [triangles[0], triangles[4]]), msg 123 144 124 145 … … 271 292 nodes_absolute = geo.get_absolute(nodes) 272 293 273 # bac, bce, ecf, dbe, daf, dae294 # bac, bce, ecf, dbe 274 295 triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 275 296 276 domain = General_mesh(nodes, triangles, 277 geo_reference = geo) 297 domain = General_mesh(nodes, triangles, geo_reference = geo) 278 298 node = domain.get_node(2) 279 self.assertEqual(c, node) 299 msg = ('\nc=%s\nnode=%s' % (str(c), str(node))) 300 self.failUnless(num.alltrue(c == node), msg) 301 302 # repeat get_node(), see if result same 303 node = domain.get_node(2) 304 msg = ('\nc=%s\nnode=%s' % (str(c), str(node))) 305 self.failUnless(num.alltrue(c == node), msg) 280 306 281 307 node = domain.get_node(2, absolute=True) 282 self.assertEqual(nodes_absolute[2], node) 283 308 msg = ('\nnodes_absolute[2]=%s\nnode=%s' 309 % (str(nodes_absolute[2]), str(node))) 310 self.failUnless(num.alltrue(nodes_absolute[2] == node), msg) 311 312 # repeat get_node(2, absolute=True), see if result same 284 313 node = domain.get_node(2, absolute=True) 285 self.assertEqual(nodes_absolute[2], node) 314 msg = ('\nnodes_absolute[2]=%s\nnode=%s' 315 % (str(nodes_absolute[2]), str(node))) 316 self.failUnless(num.alltrue(nodes_absolute[2] == node), msg) 286 317 287 318 … … 321 352 self.failUnlessRaises(AssertionError, General_mesh, 322 353 nodes, triangles, geo_reference=geo) 323 324 325 326 #------------------------------------------------------------- 354 355 ################################################################################ 356 327 357 if __name__ == "__main__": 328 suite = unittest.makeSuite(Test_General_Mesh,'test') 329 #suite = unittest.makeSuite(Test_General_Mesh,'test_get_node') 358 suite = unittest.makeSuite(Test_General_Mesh, 'test') 330 359 runner = unittest.TextTestRunner() 331 360 runner.run(suite) -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py
r6195 r7276 8 8 from anuga.config import epsilon 9 9 10 import Numericas num10 import numpy as num 11 11 12 12 … … 270 270 271 271 272 273 274 275 272 #------------------------------------------------------------- 273 276 274 if __name__ == "__main__": 277 275 suite = unittest.makeSuite(Test_Generic_Boundary_Conditions, 'test') 278 276 runner = unittest.TextTestRunner() 279 277 runner.run(suite) 280 281 282 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_ghost.py
r6195 r7276 6 6 from anuga.abstract_2d_finite_volumes.domain import * 7 7 from anuga.config import epsilon 8 9 import numpy as num 8 10 9 11 … … 40 42 41 43 42 assert domain.get_conserved_quantities(0, edge=1) == 0. 43 44 44 assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.) 45 45 46 46 47 47 #------------------------------------------------------------- 48 48 49 if __name__ == "__main__": 49 50 suite = unittest.makeSuite(Test_Domain,'test') -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py
r6717 r7276 19 19 from anuga.utilities.numerical_tools import ensure_numeric 20 20 21 import Numericas num21 import numpy as num 22 22 23 23 … … 993 993 [ 75735.4765625 , 23762.00585938], 994 994 [ 52341.70703125, 38563.39453125]] 995 996 ##points = ensure_numeric(points, Int)/1000 # Simplify for ease of interpretation997 995 998 996 triangles = [[19, 0,15], … … 1097 1095 [ 35406.3359375 , 79332.9140625 ]] 1098 1096 1099 scaled_points = ensure_numeric(points, num. Int)/1000 # Simplify for ease of interpretation1097 scaled_points = ensure_numeric(points, num.int)/1000 # Simplify for ease of interpretation 1100 1098 1101 1099 triangles = [[ 0, 1, 2], … … 1176 1174 """ 1177 1175 1178 1179 1176 # test 1180 1177 a = [0.0, 0.0] … … 1183 1180 1184 1181 absolute_points = [a, b, c] 1185 vertices = [[0, 1,2]]1186 1187 geo = Geo_reference(56, 67,-56)1182 vertices = [[0, 1, 2]] 1183 1184 geo = Geo_reference(56, 67, -56) 1188 1185 1189 1186 relative_points = geo.change_points_geo_ref(absolute_points) 1190 1191 #print 'Relative', relative_points1192 #print 'Absolute', absolute_points1193 1187 1194 1188 mesh = Mesh(relative_points, vertices, geo_reference=geo) … … 1846 1840 1847 1841 1848 1849 1850 1842 #------------------------------------------------------------- 1843 1851 1844 if __name__ == "__main__": 1852 #suite = unittest.makeSuite(Test_Mesh,'test_two_triangles') 1853 #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_partially_coinciding') 1854 #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments7') 1855 suite = unittest.makeSuite(Test_Mesh,'test') 1845 suite = unittest.makeSuite(Test_Mesh, 'test') 1856 1846 runner = unittest.TextTestRunner() 1857 1847 runner.run(suite) 1858 1859 1860 1861 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py
r6191 r7276 14 14 from anuga.pmesh.mesh import importMeshFromFile 15 15 16 import Numericas num16 import numpy as num 17 17 18 18 … … 255 255 256 256 #------------------------------------------------------------- 257 257 258 if __name__ == "__main__": 258 259 suite = unittest.makeSuite(Test_pmesh2domain, 'test') 259 260 runner = unittest.TextTestRunner() 260 261 runner.run(suite) 261 262 263 264 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
r6195 r7276 15 15 from anuga.utilities.polygon import * 16 16 17 import Numericas num17 import numpy as num 18 18 19 19 … … 504 504 quantity.set_values(0.0) 505 505 quantity.set_values(3.14, polygon=polygon) 506 msg = ('quantity.vertex_values=\n%s\nshould be close to\n' 507 '[[0,0,0],\n' 508 ' [3.14,3.14,3.14],\n' 509 ' [3.14,3.14,3.14],\n' 510 ' [0,0,0]]' % str(quantity.vertex_values)) 506 511 assert num.allclose(quantity.vertex_values, 507 512 [[0,0,0], 508 513 [3.14,3.14,3.14], 509 514 [3.14,3.14,3.14], 510 [0,0,0]]) 515 [0,0,0]]), msg 511 516 512 517 … … 1765 1770 quantity = Quantity(self.mesh4) 1766 1771 1767 quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num. Float)1772 quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.float) 1768 1773 1769 1774 quantity.interpolate_from_vertices_to_edges() … … 1781 1786 [3., 2.5, 1.5], 1782 1787 [3.5, 4.5, 3.], 1783 [2.5, 3.5, 2]], num. Float)1788 [2.5, 3.5, 2]], num.float) 1784 1789 1785 1790 quantity.interpolate_from_edges_to_vertices() … … 1835 1840 1836 1841 sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4]) 1837 denom = num.ones(4, num. Float)-timestep*sem1842 denom = num.ones(4, num.float)-timestep*sem 1838 1843 1839 1844 x = num.array([1, 2, 3, 4])/denom … … 1859 1864 1860 1865 sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4]) 1861 denom = num.ones(4, num. Float)-timestep*sem1866 denom = num.ones(4, num.float)-timestep*sem 1862 1867 1863 1868 x = num.array([1., 2., 3., 4.]) … … 1901 1906 1902 1907 bed = domain.quantities['elevation'].vertex_values 1903 stage = num.zeros(bed.shape, num. Float)1908 stage = num.zeros(bed.shape, num.float) 1904 1909 1905 1910 h = 0.03 … … 1985 1990 1986 1991 bed = domain.quantities['elevation'].vertex_values 1987 stage = num.zeros(bed.shape, num. Float)1992 stage = num.zeros(bed.shape, num.float) 1988 1993 1989 1994 h = 0.03 … … 1999 2004 stage = domain.quantities['stage'] 2000 2005 A, V = stage.get_vertex_values(xy=False, smooth=False) 2001 Q = stage.vertex_values.flat 2006 Q = stage.vertex_values.flatten() 2002 2007 2003 2008 for k in range(8): … … 2508 2513 2509 2514 #------------------------------------------------------------- 2515 2510 2516 if __name__ == "__main__": 2511 2517 suite = unittest.makeSuite(Test_Quantity, 'test') 2512 #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')2513 2514 #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_subset_and_geo')2515 #print "restricted test"2516 #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')2517 2518 runner = unittest.TextTestRunner() 2518 2519 runner.run(suite) -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py
r6145 r7276 8 8 #from anuga.config import epsilon 9 9 10 import Numericas num10 import numpy as num 11 11 12 12 … … 42 42 43 43 def test_region_tags(self): 44 """ 45 get values based on triangle lists. 46 """ 47 from mesh_factory import rectangular 48 from shallow_water import Domain 49 50 #Create basic mesh 51 points, vertices, boundary = rectangular(1, 3) 52 53 #Create shallow water domain 54 domain = Domain(points, vertices, boundary) 55 domain.build_tagged_elements_dictionary({'bottom':[0,1], 56 'top':[4,5], 57 'all':[0,1,2,3,4,5]}) 58 44 """get values based on triangle lists.""" 45 46 from mesh_factory import rectangular 47 from shallow_water import Domain 48 49 #Create basic mesh 50 points, vertices, boundary = rectangular(1, 3) 51 52 #Create shallow water domain 53 domain = Domain(points, vertices, boundary) 54 domain.build_tagged_elements_dictionary({'bottom': [0,1], 55 'top': [4,5], 56 'all': [0,1,2,3,4,5]}) 59 57 60 58 #Set friction … … 65 63 b = Set_region('top', 'friction', 1.0) 66 64 domain.set_region([a, b]) 67 #print domain.quantities['friction'].get_values() 68 assert num.allclose(domain.quantities['friction'].get_values(),\ 69 [[ 0.09, 0.09, 0.09], 70 [ 0.09, 0.09, 0.09], 71 [ 0.07, 0.07, 0.07], 72 [ 0.07, 0.07, 0.07], 73 [ 1.0, 1.0, 1.0], 74 [ 1.0, 1.0, 1.0]]) 65 66 expected = [[ 0.09, 0.09, 0.09], 67 [ 0.09, 0.09, 0.09], 68 [ 0.07, 0.07, 0.07], 69 [ 0.07, 0.07, 0.07], 70 [ 1.0, 1.0, 1.0], 71 [ 1.0, 1.0, 1.0]] 72 msg = ("\ndomain.quantities['friction']=%s\nexpected value=%s" 73 % (str(domain.quantities['friction'].get_values()), 74 str(expected))) 75 assert num.allclose(domain.quantities['friction'].get_values(), 76 expected), msg 75 77 76 78 #c = Add_Value_To_region('all', 'friction', 10.0) … … 126 128 127 129 def test_unique_vertices(self): 128 """ 129 get values based on triangle lists. 130 """ 130 """get values based on triangle lists.""" 131 131 132 from mesh_factory import rectangular 132 133 from shallow_water import Domain … … 147 148 a = Set_region('bottom', 'friction', 0.09, location = 'unique vertices') 148 149 domain.set_region(a) 149 #print domain.quantities['friction'].get_values() 150 assert num.allclose(domain.quantities['friction'].get_values(),\ 150 assert num.allclose(domain.quantities['friction'].get_values(), 151 151 [[ 0.09, 0.09, 0.09], 152 152 [ 0.09, 0.09, 0.09], … … 189 189 190 190 def test_unique_vertices_average_loc_vert(self): 191 """ 192 get values based on triangle lists. 193 """194 from mesh_factory import rectangular195 from shallow_water import Domain 196 197 #Create basic mesh198 points, vertices, boundary = rectangular(1, 3) 199 #Create shallow water domain 200 domain = Domain(points, vertices, boundary) 201 domain.build_tagged_elements_dictionary({'bottom': [0,1],202 'top': [4,5],203 'not_bottom': [2,3,4,5]})191 """Get values based on triangle lists.""" 192 193 from mesh_factory import rectangular 194 from shallow_water import Domain 195 196 #Create basic mesh 197 points, vertices, boundary = rectangular(1, 3) 198 199 #Create shallow water domain 200 domain = Domain(points, vertices, boundary) 201 domain.build_tagged_elements_dictionary({'bottom': [0, 1], 202 'top': [4, 5], 203 'not_bottom': [2, 3, 4, 5]}) 204 204 205 205 #Set friction 206 206 domain.set_quantity('friction', add_x_y) 207 av_bottom = 2.0 /3.0207 av_bottom = 2.0 / 3.0 208 208 add = 60.0 209 209 calc_frict = av_bottom + add 210 210 domain.set_region(Add_value_to_region('bottom', 'friction', add, 211 initial_quantity='friction', 212 #location='unique vertices', 213 location='vertices', 214 average=True 215 )) 216 217 #print domain.quantities['friction'].get_values() 211 initial_quantity='friction', 212 location='vertices', 213 average=True)) 214 218 215 frict_points = domain.quantities['friction'].get_values() 219 assert num.allclose(frict_points[0],\ 220 [ calc_frict, calc_frict, calc_frict]) 221 assert num.allclose(frict_points[1],\ 222 [ calc_frict, calc_frict, calc_frict]) 216 217 expected = [calc_frict, calc_frict, calc_frict] 218 msg = ('frict_points[0]=%s\nexpected=%s' 219 % (str(frict_points[0]), str(expected))) 220 assert num.allclose(frict_points[0], expected), msg 221 222 msg = ('frict_points[1]=%s\nexpected=%s' 223 % (str(frict_points[1]), str(expected))) 224 assert num.allclose(frict_points[1], expected), msg 223 225 224 226 def test_unique_vertices_average_loc_unique_vert(self): … … 261 263 262 264 #------------------------------------------------------------- 265 263 266 if __name__ == "__main__": 264 suite = unittest.makeSuite(Test_Region, 'test')267 suite = unittest.makeSuite(Test_Region, 'test') 265 268 runner = unittest.TextTestRunner() 266 269 runner.run(suite) -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py
r6318 r7276 22 22 import string 23 23 24 import Numericas num24 import numpy as num 25 25 26 26 … … 183 183 184 184 last_time_index = len(time)-1 #Last last_time_index 185 d_stage = num.reshape(num.take(stage[last_time_index, :], [0,5,10,15]), (4,1)) 186 d_uh = num.reshape(num.take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1)) 187 d_vh = num.reshape(num.take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1)) 185 d_stage = num.reshape(num.take(stage[last_time_index, :], 186 [0,5,10,15], 187 axis=0), 188 (4,1)) 189 d_uh = num.reshape(num.take(xmomentum[last_time_index, :], 190 [0,5,10,15], 191 axis=0), 192 (4,1)) 193 d_vh = num.reshape(num.take(ymomentum[last_time_index, :], 194 [0,5,10,15], 195 axis=0), 196 (4,1)) 188 197 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 189 198 … … 195 204 196 205 #And the midpoints are found now 197 Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15] )198 Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15] )206 Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15], axis=0) 207 Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15], axis=0) 199 208 200 209 diag = num.concatenate( (Dx, Dy), axis=1) … … 219 228 220 229 timestep = 0 #First timestep 221 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15] ), (4,1))222 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15] ), (4,1))223 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15] ), (4,1))230 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1)) 231 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 232 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 224 233 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 225 234 … … 241 250 242 251 timestep = 33 243 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15] ), (4,1))244 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15] ), (4,1))245 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15] ), (4,1))252 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1)) 253 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 254 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 246 255 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 247 256 … … 262 271 263 272 timestep = 15 264 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15] ), (4,1))265 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15] ), (4,1))266 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15] ), (4,1))273 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1)) 274 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 275 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 267 276 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 268 277 … … 275 284 # 276 285 timestep = 16 277 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15] ), (4,1))278 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15] ), (4,1))279 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15] ), (4,1))286 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1)) 287 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 288 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 280 289 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 281 290 … … 377 386 x = fid.variables['x'][:] 378 387 y = fid.variables['y'][:] 388 # we 'cast' to 64 bit floats to pass this test 389 # SWW file quantities are stored as 32 bits 390 x = num.array(x, num.float) 391 y = num.array(y, num.float) 392 379 393 stage = fid.variables['stage'][:] 380 394 xmomentum = fid.variables['xmomentum'][:] … … 386 400 387 401 last_time_index = len(time)-1 #Last last_time_index 388 d_stage = num.reshape(num.take(stage[last_time_index, :], [0,5,10,15]), (4,1)) 389 d_uh = num.reshape(num.take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1)) 390 d_vh = num.reshape(num.take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1)) 402 d_stage = num.reshape(num.take(stage[last_time_index, :], 403 [0,5,10,15], 404 axis=0), 405 (4,1)) 406 d_uh = num.reshape(num.take(xmomentum[last_time_index, :], 407 [0,5,10,15], 408 axis=0), 409 (4,1)) 410 d_vh = num.reshape(num.take(ymomentum[last_time_index, :], 411 [0,5,10,15], 412 axis=0), 413 (4,1)) 391 414 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 392 415 … … 398 421 399 422 #And the midpoints are found now 400 Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15] )401 Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15] )402 403 diag = num.concatenate( 423 Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15], axis=0) 424 Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15], axis=0) 425 426 diag = num.concatenate((Dx, Dy), axis=1) 404 427 d_midpoints = (diag[1:] + diag[:-1])/2 405 428 … … 415 438 416 439 t = time[last_time_index] 417 q = f(t, point_id=0); assert num.allclose(r0, q) 418 q = f(t, point_id=1); assert num.allclose(r1, q) 419 q = f(t, point_id=2); assert num.allclose(r2, q) 440 441 q = f(t, point_id=0) 442 msg = '\nr0=%s\nq=%s' % (str(r0), str(q)) 443 assert num.allclose(r0, q), msg 444 445 q = f(t, point_id=1) 446 msg = '\nr1=%s\nq=%s' % (str(r1), str(q)) 447 assert num.allclose(r1, q), msg 448 449 q = f(t, point_id=2) 450 msg = '\nr2=%s\nq=%s' % (str(r2), str(q)) 451 assert num.allclose(r2, q), msg 420 452 421 453 … … 424 456 425 457 timestep = 0 #First timestep 426 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1)) 427 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) 428 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) 458 d_stage = num.reshape(num.take(stage[timestep, :], 459 [0,5,10,15], 460 axis=0), 461 (4,1)) 462 d_uh = num.reshape(num.take(xmomentum[timestep, :], 463 [0,5,10,15], 464 axis=0), 465 (4,1)) 466 d_vh = num.reshape(num.take(ymomentum[timestep, :], 467 [0,5,10,15], 468 axis=0), 469 (4,1)) 429 470 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 430 471 … … 446 487 447 488 timestep = 33 448 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1)) 449 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) 450 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) 489 d_stage = num.reshape(num.take(stage[timestep, :], 490 [0,5,10,15], 491 axis=0), 492 (4,1)) 493 d_uh = num.reshape(num.take(xmomentum[timestep, :], 494 [0,5,10,15], 495 axis=0), 496 (4,1)) 497 d_vh = num.reshape(num.take(ymomentum[timestep, :], 498 [0,5,10,15], 499 axis=0), 500 (4,1)) 451 501 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 452 502 … … 467 517 468 518 timestep = 15 469 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1)) 470 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) 471 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) 519 d_stage = num.reshape(num.take(stage[timestep, :], 520 [0,5,10,15], 521 axis=0), 522 (4,1)) 523 d_uh = num.reshape(num.take(xmomentum[timestep, :], 524 [0,5,10,15], 525 axis=0), 526 (4,1)) 527 d_vh = num.reshape(num.take(ymomentum[timestep, :], 528 [0,5,10,15], 529 axis=0), 530 (4,1)) 472 531 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 473 532 … … 480 539 # 481 540 timestep = 16 482 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1)) 483 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) 484 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) 541 d_stage = num.reshape(num.take(stage[timestep, :], 542 [0,5,10,15], 543 axis=0), 544 (4,1)) 545 d_uh = num.reshape(num.take(xmomentum[timestep, :], 546 [0,5,10,15], 547 axis=0), 548 (4,1)) 549 d_vh = num.reshape(num.take(ymomentum[timestep, :], 550 [0,5,10,15], 551 axis=0), 552 (4,1)) 485 553 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 486 554 … … 630 698 q1 = F(t+60, point_id=id) 631 699 632 if q0 == NAN:700 if num.alltrue(q0 == NAN): 633 701 actual = q0 634 702 else: … … 641 709 #print "actual", actual 642 710 #print 643 if q0 == NAN:644 self.failUnless( q == actual, 'Fail!')711 if num.alltrue(q0 == NAN): 712 self.failUnless(num.alltrue(q == actual), 'Fail!') 645 713 else: 646 714 assert num.allclose(q, actual) … … 1181 1249 #FIXME: Division is not expected to work for integers. 1182 1250 #This must be caught. 1183 foo = num.array([[1,2,3], [4,5,6]], num. Float)1184 1185 bar = num.array([[-1,0,5], [6,1,1]], num. Float)1251 foo = num.array([[1,2,3], [4,5,6]], num.float) 1252 1253 bar = num.array([[-1,0,5], [6,1,1]], num.float) 1186 1254 1187 1255 D = {'X': foo, 'Y': bar} … … 1201 1269 1202 1270 # make an error for zero on zero 1203 # this is really an error in Numeric, SciPy core can handle it1271 # this is really an error in numeric, SciPy core can handle it 1204 1272 # Z = apply_expression_to_dictionary('0/Y', D) 1205 1273 … … 1919 1987 if 314 < angle < 316: v=1 1920 1988 assert v==1 1921 1922 1923 1924 1925 1989 1926 1990 1927 1991 #------------------------------------------------------------- 1992 1928 1993 if __name__ == "__main__": 1929 suite = unittest.makeSuite(Test_Util,'test') 1930 # suite = unittest.makeSuite(Test_Util,'test_sww2csv_gauges') 1994 suite = unittest.makeSuite(Test_Util, 'test') 1931 1995 # runner = unittest.TextTestRunner(verbosity=2) 1932 1996 runner = unittest.TextTestRunner(verbosity=1) 1933 1997 runner.run(suite) 1934 1935 1936 1937 -
anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
r7012 r7276 27 27 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a 28 28 29 import Numericas num29 import numpy as num 30 30 31 31 … … 476 476 if boundary_polygon is not None: 477 477 #removes sts points that do not lie on boundary 478 quantities[name] = num.take(quantities[name], gauge_id, 1)478 quantities[name] = num.take(quantities[name], gauge_id, axis=1) 479 479 480 480 # Close sww, tms or sts netcdf file … … 554 554 expression e.g. by overloading 555 555 556 Due to a limitation with Numeric, this can not evaluate 0/0556 Due to a limitation with numeric, this can not evaluate 0/0 557 557 In general, the user can fix by adding 1e-30 to the numerator. 558 558 SciPy core can handle this situation. … … 789 789 - easting, northing, name , elevation? 790 790 - OR (this is not yet done) 791 - structure which can be converted to a Numeric array,791 - structure which can be converted to a numeric array, 792 792 such as a geospatial data object 793 793 … … 1256 1256 n0 = int(n0) 1257 1257 m = len(locations) 1258 model_time = num.zeros((n0, m, p), num. Float)1259 stages = num.zeros((n0, m, p), num. Float)1260 elevations = num.zeros((n0, m, p), num. Float)1261 momenta = num.zeros((n0, m, p), num. Float)1262 xmom = num.zeros((n0, m, p), num. Float)1263 ymom = num.zeros((n0, m, p), num. Float)1264 speed = num.zeros((n0, m, p), num. Float)1265 bearings = num.zeros((n0, m, p), num. Float)1266 due_east = 90.0*num.ones((n0, 1), num. Float)1267 due_west = 270.0*num.ones((n0, 1), num. Float)1268 depths = num.zeros((n0, m, p), num. Float)1269 eastings = num.zeros((n0, m, p), num. Float)1258 model_time = num.zeros((n0, m, p), num.float) 1259 stages = num.zeros((n0, m, p), num.float) 1260 elevations = num.zeros((n0, m, p), num.float) 1261 momenta = num.zeros((n0, m, p), num.float) 1262 xmom = num.zeros((n0, m, p), num.float) 1263 ymom = num.zeros((n0, m, p), num.float) 1264 speed = num.zeros((n0, m, p), num.float) 1265 bearings = num.zeros((n0, m, p), num.float) 1266 due_east = 90.0*num.ones((n0, 1), num.float) 1267 due_west = 270.0*num.ones((n0, 1), num.float) 1268 depths = num.zeros((n0, m, p), num.float) 1269 eastings = num.zeros((n0, m, p), num.float) 1270 1270 min_stages = [] 1271 1271 max_stages = [] … … 1279 1279 min_speeds = [] 1280 1280 max_depths = [] 1281 model_time_plot3d = num.zeros((n0, m), num. Float)1282 stages_plot3d = num.zeros((n0, m), num. Float)1283 eastings_plot3d = num.zeros((n0, m),num. Float)1281 model_time_plot3d = num.zeros((n0, m), num.float) 1282 stages_plot3d = num.zeros((n0, m), num.float) 1283 eastings_plot3d = num.zeros((n0, m),num.float) 1284 1284 if time_unit is 'mins': scale = 60.0 1285 1285 if time_unit is 'hours': scale = 3600.0 … … 1800 1800 # Remove the loners from verts 1801 1801 # Could've used X=compress(less(loners,N),loners) 1802 # verts=num.take(verts,X ) to Remove the loners from verts1802 # verts=num.take(verts,X,axis=0) to Remove the loners from verts 1803 1803 # but I think it would use more memory 1804 1804 new_i = lone_start # point at first loner - 'shuffle down' target … … 1834 1834 """ 1835 1835 1836 xc = num.zeros(triangles.shape[0], num. Float) # Space for centroid info1836 xc = num.zeros(triangles.shape[0], num.float) # Space for centroid info 1837 1837 1838 1838 for k in range(triangles.shape[0]): … … 2121 2121 #add tide to stage if provided 2122 2122 if quantity == 'stage': 2123 quantity_value[quantity] = num.array(quantity_value[quantity], num.Float) \2124 + directory_add_tide2123 quantity_value[quantity] = num.array(quantity_value[quantity], 2124 num.float) + directory_add_tide 2125 2125 2126 2126 #condition to find max and mins for all the plots … … 2370 2370 file.close() 2371 2371 2372 2373 2372 ## 2374 2373 # @brief ?? … … 2389 2388 2390 2389 Inputs: 2391 2392 2390 NOTE: if using csv2timeseries_graphs after creating csv file, 2393 2391 it is essential to export quantities 'depth' and 'elevation'. … … 2414 2412 myfile_2_point1.csv if <out_name> ='myfile_2_' 2415 2413 2416 2417 2414 They will all have a header 2418 2415 … … 2438 2435 import string 2439 2436 from anuga.shallow_water.data_manager import get_all_swwfiles 2440 2441 # quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']2442 #print "points",points2443 2437 2444 2438 assert type(gauge_file) == type(''), 'Gauge filename must be a string' … … 2457 2451 point_name = [] 2458 2452 2459 # read point info from file2453 # read point info from file 2460 2454 for i,row in enumerate(point_reader): 2461 # read header and determine the column numbers to read correcty.2455 # read header and determine the column numbers to read correctly. 2462 2456 if i==0: 2463 2457 for j,value in enumerate(row): … … 2471 2465 2472 2466 #convert to array for file_function 2473 points_array = num.array(points,num. Float)2467 points_array = num.array(points,num.float) 2474 2468 2475 2469 points_array = ensure_absolute(points_array) … … 2525 2519 for sww_file in sww_files: 2526 2520 sww_file = join(dir_name, sww_file+'.sww') 2527 #print 'sww file = ',sww_file2528 2521 callable_sww = file_function(sww_file, 2529 2522 quantities=core_quantities, … … 2579 2572 momentum = sqrt(point_quantities[2]**2 2580 2573 + point_quantities[3]**2) 2581 # vel = momentum/depth2582 2574 vel = momentum / (point_quantities[0] 2583 2575 - point_quantities[1]) 2584 # vel = momentum/(depth + 1.e-6/depth)2585 2576 else: 2586 2577 momentum = 0 … … 2593 2584 point_quantities[3])) 2594 2585 2595 #print 'point list before write (writer %s) = %s' % (str(point_name[point_i]), str(points_list))2596 2586 points_writer[point_i].writerow(points_list) 2597 2598 2587 2599 2588 ## -
anuga_core/source/anuga/advection/advection.py
r6146 r7276 33 33 from anuga.abstract_2d_finite_volumes.domain import * 34 34 35 import Numericas num35 import numpy as num 36 36 37 37 … … 252 252 stage_bdry = Stage.boundary_values 253 253 254 flux = num.zeros(1, num. Float) #Work array for summing up fluxes254 flux = num.zeros(1, num.float) #Work array for summing up fluxes 255 255 256 256 #Loop -
anuga_core/source/anuga/advection/advection_ext.c
r5162 r7276 10 10 11 11 #include "Python.h" 12 #include " Numeric/arrayobject.h"12 #include "numpy/arrayobject.h" 13 13 #include "math.h" 14 14 #include "stdio.h" -
anuga_core/source/anuga/advection/test_advection.py
r6146 r7276 8 8 from anuga.advection.advection import Domain, Transmissive_boundary, Dirichlet_boundary 9 9 10 import Numericas num10 import numpy as num 11 11 12 12 -
anuga_core/source/anuga/alpha_shape/alpha_shape.py
r6268 r7276 25 25 from anuga.geospatial_data.geospatial_data import Geospatial_data 26 26 27 import Numericas num27 import numpy as num 28 28 29 29 … … 88 88 raise PointError, "Three points on a straight line" 89 89 90 #Convert input to Numeric arrays91 self.points = num.array(points, num. Float)90 #Convert input to numeric arrays 91 self.points = num.array(points, num.float) 92 92 93 93 … … 289 289 (denom[k]< EPSILON and denom[k] > -EPSILON)] 290 290 291 if num.a lltrue(denom != 0.0):292 dx = num.divide_safe(y31*dist21 - y21*dist31,denom)293 dy = num.divide_safe(x21*dist31 - x31*dist21,denom) 294 else:295 raise AlphaError296 291 if num.any(denom == 0.0): 292 raise AlphaError 293 294 dx = num.divide(y31*dist21 - y21*dist31, denom) 295 dy = num.divide(x21*dist31 - x31*dist21, denom) 296 297 297 self.triradius = 0.5*num.sqrt(dx*dx + dy*dy) 298 298 #print "triangle radii", self.triradius -
anuga_core/source/anuga/alpha_shape/test_alpha_shape.py
r6147 r7276 5 5 import unittest 6 6 7 import Numericas num7 import numpy as num 8 8 9 9 try: … … 394 394 (46, 45)] 395 395 assert num.allclose(answer, result) 396 396 397 #------------------------------------------------------------- 398 397 399 if __name__ == "__main__": 398 #print "starting tests \n"399 400 suite = unittest.makeSuite(TestCase,'test') 400 401 runner = unittest.TextTestRunner(verbosity=1) 401 402 runner.run(suite) 402 #print "finished tests \n"403 404 405 406 407 408 409 -
anuga_core/source/anuga/caching/caching.py
r7197 r7276 51 51 unix = True 52 52 53 import Numericas num53 import numpy as num 54 54 55 55 … … 1394 1394 I.sort() 1395 1395 val = myhash(I, ids) 1396 elif type(T) == num.ArrayType:1396 elif isinstance(T, num.ndarray): 1397 1397 T = num.array(T) # Ensure array is contiguous 1398 1398 … … 1465 1465 identical = compare(a, b, ids) 1466 1466 1467 elif type(A) == num.ArrayType:1467 elif isinstance(A, num.ndarray): 1468 1468 # Use element by element comparison 1469 1469 identical = num.alltrue(A==B) … … 2439 2439 argstr = argstr + "'"+str(args)+"'" 2440 2440 else: 2441 # Truncate large Numeric arrays before using str()2442 if type(args) == num.ArrayType:2441 # Truncate large numeric arrays before using str() 2442 if isinstance(args, num.ndarray): 2443 2443 # if len(args.flat) > textwidth: 2444 2444 # Changed by Duncan and Nick 21/2/07 .flat has problems with -
anuga_core/source/anuga/caching/test_caching.py
r6406 r7276 7 7 from anuga.caching.dummy_classes_for_testing import Dummy, Dummy_memorytest 8 8 9 import Numericas num9 import numpy as num 10 10 11 11 … … 26 26 27 27 def f_numeric(A, B): 28 """Operation on Numeric arrays28 """Operation on numeric arrays 29 29 """ 30 30 … … 123 123 """test_caching_of_numeric_arrays 124 124 125 Test that Numeric arrays can be recognised by caching even if their id's are different125 Test that numeric arrays can be recognised by caching even if their id's are different 126 126 """ 127 127 … … 160 160 161 161 162 assert T1 == T2, 'Cached result does not match computed result'163 assert T2 == T3, 'Cached result does not match computed result'162 assert num.alltrue(T1 == T2), 'Cached result does not match computed result' 163 assert num.alltrue(T2 == T3), 'Cached result does not match computed result' 164 164 165 165 166 166 def test_hash_collision(self): 167 """test_hash_collision(self): 168 169 Test that hash collisons are dealt with correctly 170 """ 167 """Test that hash collisons are dealt with correctly""" 171 168 172 169 verbose = False … … 202 199 T2 = cache(f_numeric, (A1, A1), 203 200 compression=comp, verbose=verbose) 204 201 202 T1_ref = f_numeric(A0, A0) 203 T2_ref = f_numeric(A1, A1) 205 204 T1_ref = f_numeric(A0, A0) 206 205 T2_ref = f_numeric(A1, A1) … … 419 418 420 419 421 x = num.arange(10).astype(num. Float)420 x = num.arange(10).astype(num.float) 422 421 423 422 ref1 = f1(x) … … 903 902 #------------------------------------------------------------- 904 903 if __name__ == "__main__": 905 #suite = unittest.makeSuite(Test_Caching, 'test_caching_of_callable_objects')906 904 suite = unittest.makeSuite(Test_Caching, 'test') 907 905 runner = unittest.TextTestRunner() -
anuga_core/source/anuga/config.py
r7055 r7276 4 4 import os 5 5 import sys 6 6 7 7 8 ################################################################################ … … 162 163 # Too large (100) creates 'flopping' water 163 164 # Too small (0) creates 'creep' 164 165 165 166 maximum_froude_number = 100.0 # To be used in limiters. 166 167 … … 183 184 184 185 ################################################################################ 186 # NetCDF-specific type constants. Used when defining NetCDF file variables. 187 ################################################################################ 188 189 netcdf_char = 'c' 190 netcdf_byte = 'b' 191 netcdf_int = 'i' 192 netcdf_float = 'd' 193 netcdf_float64 = 'd' 194 netcdf_float32 = 'f' 195 196 ################################################################################