Changeset 6145
- Timestamp:
- Jan 13, 2009, 11:51:22 AM (16 years ago)
- Location:
- anuga_core/source/anuga/abstract_2d_finite_volumes
- Files:
-
- 21 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r6129 r6145 8 8 """ 9 9 10 from Numeric import allclose, argmax, zeros, Float11 10 from anuga.config import epsilon 12 11 from anuga.config import beta_euler, beta_rk2 … … 31 30 from anuga.abstract_2d_finite_volumes.util import get_textual_float 32 31 33 from Numeric import zeros, Float, Int, ones34 32 from quantity import Quantity 35 33 … … 37 35 from time import time as walltime 38 36 37 import Numeric as num 39 38 40 39 … … 156 155 for key in self.full_send_dict: 157 156 buffer_shape = self.full_send_dict[key][0].shape[0] 158 self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float)) 159 157 self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys), num.Float)) 160 158 161 159 for key in self.ghost_recv_dict: 162 160 buffer_shape = self.ghost_recv_dict[key][0].shape[0] 163 self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float)) 164 161 self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys), num.Float)) 165 162 166 163 # Setup cell full flag … … 169 166 N = len(self) #number_of_elements 170 167 self.number_of_elements = N 171 self.tri_full_flag = ones(N,Int)168 self.tri_full_flag = num.ones(N, num.Int) 172 169 for i in self.ghost_recv_dict.keys(): 173 170 for id in self.ghost_recv_dict[i][0]: … … 176 173 # Test the assumption that all full triangles are store before 177 174 # the ghost triangles. 178 if not allclose(self.tri_full_flag[:self.number_of_full_nodes],1):175 if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes],1): 179 176 print 'WARNING: Not all full triangles are store before ghost triangles' 180 177 … … 229 226 # calculation 230 227 N = len(self) # Number_of_triangles 231 self.already_computed_flux = zeros((N, 3),Int)228 self.already_computed_flux = num.zeros((N, 3), num.Int) 232 229 233 230 # Storage for maximal speeds computed for each triangle by 234 231 # compute_fluxes 235 232 # This is used for diagnostics only (reset at every yieldstep) 236 self.max_speed = zeros(N,Float)233 self.max_speed = num.zeros(N, num.Float) 237 234 238 235 if mesh_filename is not None: … … 265 262 """ 266 263 267 from Numeric import zeros, Float268 269 264 if not (vertex is None or edge is None): 270 265 msg = 'Values for both vertex and edge was specified.' … … 272 267 raise msg 273 268 274 q = zeros( len(self.conserved_quantities),Float)269 q = num.zeros(len(self.conserved_quantities), num.Float) 275 270 276 271 for i, name in enumerate(self.conserved_quantities): … … 791 786 # Find index of largest computed flux speed 792 787 if triangle_id is None: 793 k = self.k = argmax(self.max_speed)788 k = self.k = num.argmax(self.max_speed) 794 789 else: 795 790 errmsg = 'Triangle_id %d does not exist in mesh: %s' %(triangle_id, … … 1224 1219 self.number_of_steps = 0 1225 1220 self.number_of_first_order_steps = 0 1226 self.max_speed = zeros(N,Float)1221 self.max_speed = num.zeros(N, num.Float) 1227 1222 1228 1223 def evolve_one_euler_step(self, yieldstep, finaltime): … … 1588 1583 """ 1589 1584 1590 from Numeric import ones, sum, equal, Float1591 1592 1585 N = len(self) # Number_of_triangles 1593 1586 d = len(self.conserved_quantities) -
anuga_core/source/anuga/abstract_2d_finite_volumes/ermapper_grids.py
r5897 r6145 2 2 3 3 # from os import open, write, read 4 import Numeric 5 6 celltype_map = {'IEEE4ByteReal': Numeric.Float32, 'IEEE8ByteReal': Numeric.Float64}4 import Numeric as num 5 6 celltype_map = {'IEEE4ByteReal': num.Float32, 'IEEE8ByteReal': num.Float64} 7 7 8 8 … … 22 22 ofile: string - filename for output (note the output will consist of two files 23 23 ofile and ofile.ers. Either of these can be entered into this function 24 data: Numeric.array - 2D array containing the data to be output to the grid24 data: num.array - 2D array containing the data to be output to the grid 25 25 header: dictionary - contains spatial information about the grid, in particular: 26 26 header['datum'] datum for the data ('"GDA94"') … … 57 57 58 58 # Check that the data is a 2 dimensional array 59 data_size = Numeric.shape(data)59 data_size = num.shape(data) 60 60 assert len(data_size) == 2 61 61 … … 83 83 nrofcellsperlines = int(header['nrofcellsperline']) 84 84 data = read_ermapper_data(data_file) 85 data = Numeric.reshape(data,(nroflines,nrofcellsperlines))85 data = num.reshape(data,(nroflines,nrofcellsperlines)) 86 86 return data 87 87 … … 163 163 return header 164 164 165 def write_ermapper_data(grid, ofile, data_format = Numeric.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 = Numeric.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 Numeric.Float32)202 grid_as_float = Numeric.fromstring(input_string,data_format)201 # convert input string to required format (Note default format is num.Float32) 202 grid_as_float = num.fromstring(input_string,data_format) 203 203 return grid_as_float 204 204 -
anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
r5897 r6145 1 2 from Numeric import concatenate, reshape, take, allclose 3 from Numeric import array, zeros, Int, Float, sqrt, sum, arange 1 import Numeric as num 4 2 5 3 from anuga.coordinate_transforms.geo_reference import Geo_reference … … 90 88 if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' 91 89 92 self.triangles = array(triangles,Int)93 self.nodes = array(nodes,Float)90 self.triangles = num.array(triangles, num.Int) 91 self.nodes = num.array(nodes, num.Float) 94 92 95 93 … … 146 144 max(self.nodes[:,0]), max(self.nodes[:,1]) ] 147 145 148 self.xy_extent = array(xy_extent,Float)146 self.xy_extent = num.array(xy_extent, num.Float) 149 147 150 148 151 149 # Allocate space for geometric quantities 152 self.normals = zeros((N, 6),Float)153 self.areas = zeros(N,Float)154 self.edgelengths = zeros((N, 3),Float)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) 155 153 156 154 # Get x,y coordinates for all triangles and store … … 187 185 # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle 188 186 189 n0 = array([x2 - x1, y2 - y1])190 l0 = sqrt(sum(n0**2))191 192 n1 = array([x0 - x2, y0 - y2])193 l1 = sqrt(sum(n1**2))194 195 n2 = array([x1 - x0, y1 - y0])196 l2 = sqrt(sum(n2**2))187 n0 = num.array([x2 - x1, y2 - y1]) 188 l0 = num.sqrt(num.sum(n0**2)) 189 190 n1 = num.array([x0 - x2, y0 - y2]) 191 l1 = num.sqrt(num.sum(n1**2)) 192 193 n2 = num.array([x1 - x0, y1 - y0]) 194 l2 = num.sqrt(num.sum(n2**2)) 197 195 198 196 # Normalise … … 274 272 if absolute is True: 275 273 if not self.geo_reference.is_absolute(): 276 return V + array([self.geo_reference.get_xllcorner(),277 self.geo_reference.get_yllcorner()])274 return V + num.array([self.geo_reference.get_xllcorner(), 275 self.geo_reference.get_yllcorner()]) 278 276 else: 279 277 return V … … 316 314 i3 = 3*i 317 315 if absolute is True and not self.geo_reference.is_absolute(): 318 offset= array([self.geo_reference.get_xllcorner(),316 offset=num.array([self.geo_reference.get_xllcorner(), 319 317 self.geo_reference.get_yllcorner()]) 320 return array([V[i3,:]+offset,321 V[i3+1,:]+offset,322 V[i3+2,:]+offset])318 return num.array([V[i3,:]+offset, 319 V[i3+1,:]+offset, 320 V[i3+2,:]+offset]) 323 321 else: 324 return array([V[i3,:], V[i3+1,:], V[i3+2,:]])322 return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]]) 325 323 326 324 … … 349 347 350 348 M = self.number_of_triangles 351 vertex_coordinates = zeros((3*M, 2),Float)349 vertex_coordinates = num.zeros((3*M, 2), num.Float) 352 350 353 351 for i in range(M): … … 376 374 indices = range(M) 377 375 378 return take(self.triangles, indices)376 return num.take(self.triangles, indices) 379 377 380 378 … … 408 406 K = 3*M # Total number of unique vertices 409 407 410 #T = reshape(array(range(K)).astype(Int), (M,3)) 411 T = reshape(arange(K).astype(Int), (M,3)) # Faster 408 T = num.reshape(num.arange(K).astype(num.Int), (M,3)) 412 409 413 410 return T … … 439 436 if node is not None: 440 437 # Get index for this node 441 first = sum(self.number_of_triangles_per_node[:node])438 first = num.sum(self.number_of_triangles_per_node[:node]) 442 439 443 440 # Get number of triangles for this node … … 452 449 triangle_list.append( (volume_id, vertex_id) ) 453 450 454 triangle_list = array(triangle_list)451 triangle_list = num.array(triangle_list) 455 452 else: 456 453 # Get info for all nodes recursively. … … 529 526 530 527 # Count number of triangles per node 531 number_of_triangles_per_node = zeros(self.number_of_full_nodes)528 number_of_triangles_per_node = num.zeros(self.number_of_full_nodes) 532 529 for volume_id, triangle in enumerate(self.get_triangles()): 533 530 for vertex_id in triangle: … … 535 532 536 533 # Allocate space for inverted structure 537 number_of_entries = sum(number_of_triangles_per_node)538 vertex_value_indices = zeros(number_of_entries)534 number_of_entries = num.sum(number_of_triangles_per_node) 535 vertex_value_indices = num.zeros(number_of_entries) 539 536 540 537 # Register (triangle, vertex) indices for each node … … 598 595 """ 599 596 600 return sum(self.areas)601 602 603 597 return num.sum(self.areas) 598 599 600 -
anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
r6011 r6145 7 7 from anuga.fit_interpolate.interpolate import Modeltime_too_late 8 8 from anuga.fit_interpolate.interpolate import Modeltime_too_early 9 10 import Numeric as num 9 11 10 12 … … 67 69 raise Exception, msg 68 70 69 from Numeric import array, Float 70 self.conserved_quantities=array(conserved_quantities).astype(Float) 71 self.conserved_quantities=num.array(conserved_quantities).astype(num.Float) 71 72 72 73 def __repr__(self): … … 97 98 98 99 99 from Numeric import array, Float100 100 try: 101 q = array(q).astype(Float)101 q = num.array(q).astype(num.Float) 102 102 except: 103 103 msg = 'Return value from time boundary function could ' … … 164 164 165 165 import time 166 from Numeric import array, zeros, Float167 166 from anuga.config import time_format 168 167 from anuga.abstract_2d_finite_volumes.util import file_function … … 180 179 181 180 if verbose: print 'Find midpoint coordinates of entire boundary' 182 self.midpoint_coordinates = zeros( (len(domain.boundary), 2),Float)181 self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.Float) 183 182 boundary_keys = domain.boundary.keys() 184 183 … … 200 199 201 200 # Compute midpoints 202 if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])203 if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])204 if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])201 if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2]) 202 if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2]) 203 if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2]) 205 204 206 205 # Convert to absolute UTM coordinates -
anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py
r5897 r6145 2 2 mesh file formats 3 3 """ 4 5 import Numeric as num 4 6 5 7 … … 73 75 74 76 from anuga.config import epsilon 75 from Numeric import zeros, Float, Int76 77 77 78 delta1 = float(len1)/m … … 93 94 index = Index(n,m) 94 95 95 points = zeros( (Np,2),Float)96 points = num.zeros((Np, 2), num.Float) 96 97 97 98 for i in range(m+1): … … 105 106 106 107 107 elements = zeros( (Nt,3),Int)108 elements = num.zeros((Nt, 3), num.Int) 108 109 boundary = {} 109 110 nt = -1 … … 149 150 150 151 from anuga.config import epsilon 151 from Numeric import zeros, Float, Int152 152 153 153 delta1 = float(len1)/m … … 208 208 """ 209 209 210 from Numeric import array211 210 import math 212 211 … … 510 509 """ 511 510 512 from Numeric import array513 511 import math 514 512 … … 592 590 """ 593 591 594 from Numeric import array595 592 import math 596 593 … … 686 683 """ 687 684 688 from Numeric import array689 685 import math 690 686 -
anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
r5897 r6145 11 11 from anuga.caching import cache 12 12 from math import pi, sqrt 13 from Numeric import array, allclose 13 14 import Numeric as num 14 15 15 16 … … 81 82 82 83 83 from Numeric import array, zeros, Int, Float, maximum, sqrt, sum84 84 85 85 General_mesh.__init__(self, coordinates, triangles, … … 98 98 99 99 #Allocate space for geometric quantities 100 self.centroid_coordinates = zeros((N, 2),Float)101 102 self.radii = zeros(N,Float)103 104 self.neighbours = zeros((N, 3),Int)105 self.neighbour_edges = zeros((N, 3),Int)106 self.number_of_boundaries = zeros(N,Int)107 self.surrogate_neighbours = zeros((N, 3),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 … … 124 124 125 125 #Compute centroid 126 centroid = array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])126 centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3]) 127 127 self.centroid_coordinates[i] = centroid 128 128 … … 133 133 134 134 #Midpoints 135 m0 = array([(x1 + x2)/2, (y1 + y2)/2])136 m1 = array([(x0 + x2)/2, (y0 + y2)/2])137 m2 = array([(x1 + x0)/2, (y1 + y0)/2])135 m0 = num.array([(x1 + x2)/2, (y1 + y2)/2]) 136 m1 = num.array([(x0 + x2)/2, (y0 + y2)/2]) 137 m2 = num.array([(x1 + x0)/2, (y1 + y0)/2]) 138 138 139 139 #The radius is the distance from the centroid of 140 140 #a triangle to the midpoint of the side of the triangle 141 141 #closest to the centroid 142 d0 = sqrt(sum( (centroid-m0)**2 ))143 d1 = sqrt(sum( (centroid-m1)**2 ))144 d2 = sqrt(sum( (centroid-m2)**2 ))142 d0 = num.sqrt(num.sum( (centroid-m0)**2 )) 143 d1 = num.sqrt(num.sum( (centroid-m1)**2 )) 144 d2 = num.sqrt(num.sum( (centroid-m2)**2 )) 145 145 146 146 self.radii[i] = min(d0, d1, d2) … … 150 150 #of inscribed circle is computed 151 151 152 a = sqrt((x0-x1)**2+(y0-y1)**2)153 b = sqrt((x1-x2)**2+(y1-y2)**2)154 c = sqrt((x2-x0)**2+(y2-y0)**2)152 a = num.sqrt((x0-x1)**2+(y0-y1)**2) 153 b = num.sqrt((x1-x2)**2+(y1-y2)**2) 154 c = num.sqrt((x2-x0)**2+(y2-y0)**2) 155 155 156 156 self.radii[i]=2.0*self.areas[i]/(a+b+c) … … 212 212 x1 = V[i, 2]; y1 = V[i, 3] 213 213 x2 = V[i, 4]; y2 = V[i, 5] 214 a = sqrt((x0-x1)**2+(y0-y1)**2)215 b = sqrt((x1-x2)**2+(y1-y2)**2)216 c = sqrt((x2-x0)**2+(y2-y0)**2)214 a = num.sqrt((x0-x1)**2+(y0-y1)**2) 215 b = num.sqrt((x1-x2)**2+(y1-y2)**2) 216 c = num.sqrt((x2-x0)**2+(y2-y0)**2) 217 217 ratio = old_rad/self.radii[i] 218 218 max_ratio = ratio … … 224 224 x1 = V[i, 2]; y1 = V[i, 3] 225 225 x2 = V[i, 4]; y2 = V[i, 5] 226 a = sqrt((x0-x1)**2+(y0-y1)**2)227 b = sqrt((x1-x2)**2+(y1-y2)**2)228 c = sqrt((x2-x0)**2+(y2-y0)**2)226 a = num.sqrt((x0-x1)**2+(y0-y1)**2) 227 b = num.sqrt((x1-x2)**2+(y1-y2)**2) 228 c = num.sqrt((x2-x0)**2+(y2-y0)**2) 229 229 self.radii[i]=self.areas[i]/(2*(a+b+c))*safety_factor 230 230 ratio = old_rad/self.radii[i] … … 388 388 self.element_tag is defined 389 389 """ 390 from Numeric import array, Int391 390 392 391 if tagged_elements is None: … … 395 394 #Check that all keys in given boundary exist 396 395 for tag in tagged_elements.keys(): 397 tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)396 tagged_elements[tag] = num.array(tagged_elements[tag]).astype(num.Int) 398 397 399 398 msg = 'Not all elements exist. ' … … 463 462 """ 464 463 465 from Numeric import allclose, sqrt, array, minimum, maximum466 464 from anuga.utilities.numerical_tools import angle, ensure_numeric 467 465 … … 477 475 inverse_segments = {} 478 476 p0 = None 479 mindist = sqrt(sum((pmax-pmin)**2)) # Start value across entire mesh477 mindist = num.sqrt(num.sum((pmax-pmin)**2)) # Start value across entire mesh 480 478 for i, edge_id in self.boundary.keys(): 481 479 # Find vertex ids for boundary segment … … 491 489 # Note: Could be arbitrary, but nice to have 492 490 # a unique way of selecting 493 dist_A = sqrt(sum((A-pmin)**2))494 dist_B = sqrt(sum((B-pmin)**2))491 dist_A = num.sqrt(num.sum((A-pmin)**2)) 492 dist_B = num.sqrt(num.sum((B-pmin)**2)) 495 493 496 494 # Find lower leftmost point … … 593 591 # We have reached a point already visited. 594 592 595 if allclose(p1, polygon[0]):593 if num.allclose(p1, polygon[0]): 596 594 # If it is the initial point, the polygon is complete. 597 595 … … 629 627 from anuga.utilities.numerical_tools import anglediff 630 628 631 from Numeric import sort, allclose632 633 629 N = len(self) 634 630 … … 672 668 673 669 # Normalise 674 l_u = sqrt(u[0]*u[0] + u[1]*u[1])675 l_v = sqrt(v[0]*v[0] + v[1]*v[1])670 l_u = num.sqrt(u[0]*u[0] + u[1]*u[1]) 671 l_v = num.sqrt(v[0]*v[0] + v[1]*v[1]) 676 672 677 673 msg = 'Normal vector in triangle %d does not have unit length' %i 678 assert allclose(l_v, 1), msg674 assert num.allclose(l_v, 1), msg 679 675 680 676 x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product … … 730 726 731 727 V = self.vertex_value_indices[:] #Take a copy 732 V = sort(V)733 assert allclose(V, range(3*N))734 735 assert sum(self.number_of_triangles_per_node) ==\736 len(self.vertex_value_indices)728 V = num.sort(V) 729 assert num.allclose(V, range(3*N)) 730 731 assert num.sum(self.number_of_triangles_per_node) ==\ 732 len(self.vertex_value_indices) 737 733 738 734 # Check number of triangles per node … … 742 738 count[i] += 1 743 739 744 assert allclose(count, self.number_of_triangles_per_node)740 assert num.allclose(count, self.number_of_triangles_per_node) 745 741 746 742 … … 805 801 """ 806 802 807 from Numeric import arange808 803 from anuga.utilities.numerical_tools import histogram, create_bins 809 804 … … 1141 1136 1142 1137 # Distances from line origin to the two intersections 1143 z0 = array([x0 - xi0, y0 - eta0])1144 z1 = array([x1 - xi0, y1 - eta0])1145 d0 = sqrt(sum(z0**2))1146 d1 = sqrt(sum(z1**2))1138 z0 = num.array([x0 - xi0, y0 - eta0]) 1139 z1 = num.array([x1 - xi0, y1 - eta0]) 1140 d0 = num.sqrt(num.sum(z0**2)) 1141 d1 = num.sqrt(num.sum(z1**2)) 1147 1142 1148 1143 if d1 < d0: … … 1157 1152 # Normal direction: 1158 1153 # Right hand side relative to line direction 1159 vector = array([x1 - x0, y1 - y0]) # Segment vector1160 length = sqrt(sum(vector**2)) # Segment length1161 normal = array([vector[1], -vector[0]])/length1154 vector = num.array([x1 - x0, y1 - y0]) # Segment vector 1155 length = num.sqrt(num.sum(vector**2)) # Segment length 1156 normal = num.array([vector[1], -vector[0]])/length 1162 1157 1163 1158 … … 1239 1234 assert isinstance(segment, Triangle_intersection), msg 1240 1235 1241 midpoint = sum(array(segment.segment))/21236 midpoint = num.sum(num.array(segment.segment))/2 1242 1237 midpoints.append(midpoint) 1243 1238 -
anuga_core/source/anuga/abstract_2d_finite_volumes/pmesh2domain.py
r5897 r6145 9 9 import sys 10 10 11 import Numeric as num 11 12 12 13 … … 161 162 """ 162 163 163 from Numeric import transpose164 164 from load_mesh.loadASCII import import_mesh_file 165 165 … … 172 172 volumes = mesh_dict['triangles'] 173 173 vertex_quantity_dict = {} 174 point_atts = transpose(mesh_dict['vertex_attributes'])174 point_atts = num.transpose(mesh_dict['vertex_attributes']) 175 175 point_titles = mesh_dict['vertex_attribute_titles'] 176 176 geo_reference = mesh_dict['geo_reference'] -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
r5990 r6145 15 15 """ 16 16 17 from Numeric import array, zeros, Float, less, concatenate, NewAxis,\18 argmax, argmin, allclose, take, reshape, alltrue19 20 17 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar 21 18 from anuga.utilities.polygon import inside_polygon … … 26 23 from anuga.config import epsilon 27 24 25 import Numeric as num 26 27 28 28 class Quantity: 29 29 … … 38 38 if vertex_values is None: 39 39 N = len(domain) # number_of_elements 40 self.vertex_values = zeros((N, 3),Float)40 self.vertex_values = num.zeros((N, 3), num.Float) 41 41 else: 42 self.vertex_values = array(vertex_values).astype(Float)42 self.vertex_values = num.array(vertex_values).astype(num.Float) 43 43 44 44 N, V = self.vertex_values.shape … … 57 57 58 58 # Allocate space for other quantities 59 self.centroid_values = zeros(N,Float)60 self.edge_values = zeros((N, 3),Float)59 self.centroid_values = num.zeros(N, num.Float) 60 self.edge_values = num.zeros((N, 3), num.Float) 61 61 62 62 # Allocate space for Gradient 63 self.x_gradient = zeros(N,Float)64 self.y_gradient = zeros(N,Float)63 self.x_gradient = num.zeros(N, num.Float) 64 self.y_gradient = num.zeros(N, num.Float) 65 65 66 66 # Allocate space for Limiter Phi 67 self.phi = zeros(N,Float)67 self.phi = num.zeros(N, num.Float) 68 68 69 69 # Intialise centroid and edge_values … … 72 72 # Allocate space for boundary values 73 73 L = len(domain.boundary) 74 self.boundary_values = zeros(L,Float)74 self.boundary_values = num.zeros(L, num.Float) 75 75 76 76 # Allocate space for updates of conserved quantities by … … 78 78 79 79 # Allocate space for update fields 80 self.explicit_update = zeros(N,Float )81 self.semi_implicit_update = zeros(N,Float )82 self.centroid_backup_values = zeros(N,Float)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) 83 83 84 84 self.set_beta(1.0) … … 382 382 from anuga.geospatial_data.geospatial_data import Geospatial_data 383 383 from types import FloatType, IntType, LongType, ListType, NoneType 384 from Numeric import ArrayType385 384 386 385 # Treat special case: Polygon situation … … 450 449 451 450 msg = 'Indices must be a list or None' 452 assert type(indices) in [ListType, NoneType, ArrayType], msg451 assert type(indices) in [ListType, NoneType, num.ArrayType], msg 453 452 454 453 … … 460 459 self.set_values_from_constant(numeric, 461 460 location, indices, verbose) 462 elif type(numeric) in [ ArrayType, ListType]:461 elif type(numeric) in [num.ArrayType, ListType]: 463 462 self.set_values_from_array(numeric, 464 463 location, indices, verbose) … … 612 611 """ 613 612 614 from Numeric import array, Float, Int, allclose 615 616 values = array(values).astype(Float) 613 values = num.array(values).astype(num.Float) 617 614 618 615 if indices is not None: 619 indices = array(indices).astype(Int)616 indices = num.array(indices).astype(num.Int) 620 617 msg = 'Number of values must match number of indices:' 621 618 msg += ' You specified %d values and %d indices'\ … … 642 639 643 640 elif location == 'unique vertices': 644 assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\641 assert len(values.shape) == 1 or num.allclose(values.shape[1:], 1),\ 645 642 'Values array must be 1d' 646 643 … … 678 675 A = q.vertex_values 679 676 680 from Numeric import allclose681 677 msg = 'Quantities are defined on different meshes. '+\ 682 678 'This might be a case for implementing interpolation '+\ 683 679 'between different meshes.' 684 assert allclose(A.shape, self.vertex_values.shape), msg680 assert num.allclose(A.shape, self.vertex_values.shape), msg 685 681 686 682 self.set_values(A, location='vertices', … … 718 714 indices = range(len(self)) 719 715 720 V = take(self.domain.get_centroid_coordinates(), indices)716 V = num.take(self.domain.get_centroid_coordinates(), indices) 721 717 self.set_values(f(V[:,0], V[:,1]), 722 718 location=location, … … 738 734 # more robust and simple. 739 735 740 #values = reshape(values, (M,3))736 #values = num.reshape(values, (M,3)) 741 737 #self.set_values(values, 742 738 # location='vertices', … … 782 778 783 779 784 points = ensure_numeric(points, Float)785 values = ensure_numeric(values, Float)780 points = ensure_numeric(points, num.Float) 781 values = ensure_numeric(values, num.Float) 786 782 787 783 if location != 'vertices': … … 915 911 # Always return absolute indices 916 912 if mode is None or mode == 'max': 917 i = argmax(V)913 i = num.argmax(V) 918 914 elif mode == 'min': 919 i = argmin(V)915 i = num.argmin(V) 920 916 921 917 … … 1120 1116 """ 1121 1117 1122 from Numeric import take1123 1118 1124 1119 # FIXME (Ole): I reckon we should have the option of passing a … … 1145 1140 raise msg 1146 1141 1147 import types , Numeric1142 import types 1148 1143 assert type(indices) in [types.ListType, types.NoneType, 1149 Numeric.ArrayType],\1144 num.ArrayType],\ 1150 1145 'Indices must be a list or None' 1151 1146 … … 1153 1148 if (indices == None): 1154 1149 indices = range(len(self)) 1155 return take(self.centroid_values,indices)1150 return num.take(self.centroid_values,indices) 1156 1151 elif location == 'edges': 1157 1152 if (indices == None): 1158 1153 indices = range(len(self)) 1159 return take(self.edge_values,indices)1154 return num.take(self.edge_values,indices) 1160 1155 elif location == 'unique vertices': 1161 1156 if (indices == None): … … 1180 1175 sum += self.vertex_values[triangle_id, vertex_id] 1181 1176 vert_values.append(sum/len(triangles)) 1182 return Numeric.array(vert_values)1177 return num.array(vert_values) 1183 1178 else: 1184 1179 if (indices is None): 1185 1180 indices = range(len(self)) 1186 return take(self.vertex_values, indices)1181 return num.take(self.vertex_values, indices) 1187 1182 1188 1183 … … 1198 1193 """ 1199 1194 1200 from Numeric import array, Float1201 1195 1202 1196 # Assert that A can be converted to a Numeric array of appropriate dim 1203 A = ensure_numeric(A, Float)1197 A = ensure_numeric(A, num.Float) 1204 1198 1205 1199 # print 'SHAPE A', A.shape … … 1275 1269 """ 1276 1270 1277 from Numeric import concatenate, zeros, Float, Int, array, reshape1278 1279 1271 1280 1272 if smooth is None: … … 1286 1278 1287 1279 if precision is None: 1288 precision = Float1280 precision = num.Float 1289 1281 1290 1282 … … 1295 1287 V = self.domain.get_triangles() 1296 1288 N = self.domain.number_of_full_nodes # Ignore ghost nodes if any 1297 A = zeros(N,Float)1289 A = num.zeros(N, num.Float) 1298 1290 points = self.domain.get_nodes() 1299 1291 -
anuga_core/source/anuga/abstract_2d_finite_volumes/region.py
r5897 r6145 8 8 # FIXME (DSG-DSG) add better comments 9 9 10 from Numeric import average 10 import Numeric as num 11 12 11 13 class Region: 12 14 """Base class for modifying quantities based on a region. … … 101 103 values = Q.get_values(indices=self.build_indices(elements, domain), 102 104 location=self.location) 103 av = average(values)105 av = num.average(values) 104 106 if self.location == "vertices": 105 av = average(av)107 av = num.average(av) 106 108 new_values = av + self.X 107 109 else: -
anuga_core/source/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py
r5897 r6145 18 18 19 19 from mesh_factory import rectangular 20 from Numeric import array 21 20 22 21 23 22 ###################### … … 78 77 domain.set_quantity('stage', Z) 79 78 80 from Numeric import allclose81 82 79 #Evolve 83 80 for t in domain.evolve(yieldstep = 0.1, finaltime = 30): -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py
r6129 r6145 6 6 from domain import * 7 7 from anuga.config import epsilon 8 from Numeric import allclose, array, ones, Float 8 9 import Numeric as num 9 10 10 11 … … 95 96 #Centroids 96 97 q = domain.get_conserved_quantities(0) 97 assert allclose(q, [2., 2., 0.])98 assert num.allclose(q, [2., 2., 0.]) 98 99 99 100 q = domain.get_conserved_quantities(1) 100 assert allclose(q, [5., 5., 0.])101 assert num.allclose(q, [5., 5., 0.]) 101 102 102 103 q = domain.get_conserved_quantities(2) 103 assert allclose(q, [3., 3., 0.])104 assert num.allclose(q, [3., 3., 0.]) 104 105 105 106 q = domain.get_conserved_quantities(3) 106 assert allclose(q, [0., 0., 0.])107 assert num.allclose(q, [0., 0., 0.]) 107 108 108 109 109 110 #Edges 110 111 q = domain.get_conserved_quantities(0, edge=0) 111 assert allclose(q, [2.5, 2.5, 0.])112 assert num.allclose(q, [2.5, 2.5, 0.]) 112 113 q = domain.get_conserved_quantities(0, edge=1) 113 assert allclose(q, [2., 2., 0.])114 assert num.allclose(q, [2., 2., 0.]) 114 115 q = domain.get_conserved_quantities(0, edge=2) 115 assert allclose(q, [1.5, 1.5, 0.])116 assert num.allclose(q, [1.5, 1.5, 0.]) 116 117 117 118 for i in range(3): 118 119 q = domain.get_conserved_quantities(1, edge=i) 119 assert allclose(q, [5, 5, 0.])120 assert num.allclose(q, [5, 5, 0.]) 120 121 121 122 122 123 q = domain.get_conserved_quantities(2, edge=0) 123 assert allclose(q, [4.5, 4.5, 0.])124 assert num.allclose(q, [4.5, 4.5, 0.]) 124 125 q = domain.get_conserved_quantities(2, edge=1) 125 assert allclose(q, [4.5, 4.5, 0.])126 assert num.allclose(q, [4.5, 4.5, 0.]) 126 127 q = domain.get_conserved_quantities(2, edge=2) 127 assert allclose(q, [0., 0., 0.])128 assert num.allclose(q, [0., 0., 0.]) 128 129 129 130 130 131 q = domain.get_conserved_quantities(3, edge=0) 131 assert allclose(q, [3., 3., 0.])132 assert num.allclose(q, [3., 3., 0.]) 132 133 q = domain.get_conserved_quantities(3, edge=1) 133 assert allclose(q, [-1.5, -1.5, 0.])134 assert num.allclose(q, [-1.5, -1.5, 0.]) 134 135 q = domain.get_conserved_quantities(3, edge=2) 135 assert allclose(q, [-1.5, -1.5, 0.])136 assert num.allclose(q, [-1.5, -1.5, 0.]) 136 137 137 138 … … 179 180 Q = domain.create_quantity_from_expression(expression) 180 181 181 assert allclose(Q.vertex_values, [[2,3,4], [6,6,6],182 [1,1,10], [-5, 4, 4]])182 assert num.allclose(Q.vertex_values, [[2,3,4], [6,6,6], 183 [1,1,10], [-5, 4, 4]]) 183 184 184 185 expression = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5' … … 188 189 Y = domain.quantities['ymomentum'].vertex_values 189 190 190 assert allclose(Q.vertex_values, (X**2 + Y**2)**0.5)191 assert num.allclose(Q.vertex_values, (X**2 + Y**2)**0.5) 191 192 192 193 … … 314 315 Q = domain.quantities['depth'] 315 316 316 assert allclose(Q.vertex_values, [[2,3,4], [6,6,6],317 [1,1,10], [-5, 4, 4]])317 assert num.allclose(Q.vertex_values, [[2,3,4], [6,6,6], 318 [1,1,10], [-5, 4, 4]]) 318 319 319 320 … … 344 345 345 346 346 A = array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], 'f')347 B = array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], 'f')347 A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], 'f') 348 B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], 'f') 348 349 349 350 #print A … … 359 360 domain.set_quantity('elevation', A) 360 361 domain.add_quantity('elevation', B) 361 assert allclose(elevation.vertex_values, A+B)362 assert num.allclose(elevation.vertex_values, A+B) 362 363 363 364 domain.add_quantity('elevation', 4) 364 assert allclose(elevation.vertex_values, A+B+4)365 assert num.allclose(elevation.vertex_values, A+B+4) 365 366 366 367 … … 370 371 domain.set_quantity('depth', 1.0) 371 372 domain.add_quantity('depth', expression = 'stage - elevation') 372 assert allclose(depth.vertex_values, stage.vertex_values-elevation.vertex_values+1)373 assert num.allclose(depth.vertex_values, stage.vertex_values-elevation.vertex_values+1) 373 374 374 375 … … 376 377 reference = 2*stage.vertex_values - depth.vertex_values 377 378 domain.add_quantity('stage', expression = 'stage - depth') 378 assert allclose(stage.vertex_values, reference)379 assert num.allclose(stage.vertex_values, reference) 379 380 380 381 … … 388 389 389 390 domain.add_quantity('depth', f) 390 assert allclose(stage.vertex_values, depth.vertex_values)391 assert num.allclose(stage.vertex_values, depth.vertex_values) 391 392 392 393 … … 458 459 domain.check_integrity() 459 460 460 assert allclose(domain.neighbours, [[-1,-2,-3]])461 assert num.allclose(domain.neighbours, [[-1,-2,-3]]) 461 462 462 463 … … 547 548 [0,0,9], [-6, 3, 3]]) 548 549 549 assert allclose( domain.quantities['stage'].centroid_values,550 [2,5,3,0] )550 assert num.allclose( domain.quantities['stage'].centroid_values, 551 [2,5,3,0] ) 551 552 552 553 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], … … 560 561 561 562 #First order extrapolation 562 assert allclose( domain.quantities['stage'].vertex_values,563 [[ 2., 2., 2.],564 [ 5., 5., 5.],565 [ 3., 3., 3.],566 [ 0., 0., 0.]])563 assert num.allclose( domain.quantities['stage'].vertex_values, 564 [[ 2., 2., 2.], 565 [ 5., 5., 5.], 566 [ 3., 3., 3.], 567 [ 0., 0., 0.]]) 567 568 568 569 … … 603 604 604 605 for name in domain.conserved_quantities: 605 domain.quantities[name].explicit_update = array([4.,3.,2.,1.])606 domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.])606 domain.quantities[name].explicit_update = num.array([4.,3.,2.,1.]) 607 domain.quantities[name].semi_implicit_update = num.array([1.,1.,1.,1.]) 607 608 608 609 … … 611 612 domain.update_conserved_quantities() 612 613 613 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])614 denom = ones(4,Float)-domain.timestep*sem614 sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4]) 615 denom = num.ones(4, num.Float)-domain.timestep*sem 615 616 616 617 # x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) 617 618 # x /= denom 618 619 619 x = array([1., 2., 3., 4.])620 x = num.array([1., 2., 3., 4.]) 620 621 x /= denom 621 x += domain.timestep* array( [4,3,2,1] )622 x += domain.timestep*num.array( [4,3,2,1] ) 622 623 623 624 for name in domain.conserved_quantities: 624 assert allclose(domain.quantities[name].centroid_values, x)625 assert num.allclose(domain.quantities[name].centroid_values, x) 625 626 626 627 … … 654 655 [0,0,9], [-6, 3, 3]]) 655 656 656 assert allclose( domain.quantities['stage'].centroid_values,657 [2,5,3,0] )657 assert num.allclose( domain.quantities['stage'].centroid_values, 658 [2,5,3,0] ) 658 659 659 660 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], … … 667 668 668 669 #First order extrapolation 669 assert allclose( domain.quantities['stage'].vertex_values,670 [[ 2., 2., 2.],671 [ 5., 5., 5.],672 [ 3., 3., 3.],673 [ 0., 0., 0.]])670 assert num.allclose( domain.quantities['stage'].vertex_values, 671 [[ 2., 2., 2.], 672 [ 5., 5., 5.], 673 [ 3., 3., 3.], 674 [ 0., 0., 0.]]) 674 675 675 676 domain.build_tagged_elements_dictionary({'mound':[0,1]}) … … 686 687 from mesh_factory import rectangular 687 688 from shallow_water import Domain 688 from Numeric import zeros, Float689 689 690 690 #Create basic mesh … … 703 703 from mesh_factory import rectangular 704 704 from shallow_water import Domain 705 from Numeric import zeros, Float706 705 707 706 #Create basic mesh … … 721 720 domain.set_region([set_bottom_friction, set_top_friction]) 722 721 #print domain.quantities['friction'].get_values() 723 assert allclose(domain.quantities['friction'].get_values(),\724 [[ 0.09, 0.09, 0.09],725 [ 0.09, 0.09, 0.09],726 [ 0.07, 0.07, 0.07],727 [ 0.07, 0.07, 0.07],728 [ 1.0, 1.0, 1.0],729 [ 1.0, 1.0, 1.0]])722 assert num.allclose(domain.quantities['friction'].get_values(),\ 723 [[ 0.09, 0.09, 0.09], 724 [ 0.09, 0.09, 0.09], 725 [ 0.07, 0.07, 0.07], 726 [ 0.07, 0.07, 0.07], 727 [ 1.0, 1.0, 1.0], 728 [ 1.0, 1.0, 1.0]]) 730 729 731 730 domain.set_region([set_all_friction]) 732 731 #print domain.quantities['friction'].get_values() 733 assert allclose(domain.quantities['friction'].get_values(),734 [[ 10.09, 10.09, 10.09],735 [ 10.09, 10.09, 10.09],736 [ 10.07, 10.07, 10.07],737 [ 10.07, 10.07, 10.07],738 [ 11.0, 11.0, 11.0],739 [ 11.0, 11.0, 11.0]])732 assert num.allclose(domain.quantities['friction'].get_values(), 733 [[ 10.09, 10.09, 10.09], 734 [ 10.09, 10.09, 10.09], 735 [ 10.07, 10.07, 10.07], 736 [ 10.07, 10.07, 10.07], 737 [ 11.0, 11.0, 11.0], 738 [ 11.0, 11.0, 11.0]]) 740 739 741 740 … … 746 745 from mesh_factory import rectangular 747 746 from shallow_water import Domain 748 from Numeric import zeros, Float749 747 750 748 #Create basic mesh … … 766 764 767 765 #print domain.quantities['friction'].get_values() 768 assert allclose(domain.quantities['friction'].get_values(),\769 [[ 0.09, 0.09, 0.09],770 [ 0.09, 0.09, 0.09],771 [ 0.07, 0.07, 0.07],772 [ 0.07, 0.07, 0.07],773 [ 1.0, 1.0, 1.0],774 [ 1.0, 1.0, 1.0]])766 assert num.allclose(domain.quantities['friction'].get_values(), 767 [[ 0.09, 0.09, 0.09], 768 [ 0.09, 0.09, 0.09], 769 [ 0.07, 0.07, 0.07], 770 [ 0.07, 0.07, 0.07], 771 [ 1.0, 1.0, 1.0], 772 [ 1.0, 1.0, 1.0]]) 775 773 776 774 domain.set_region([set_bottom_friction, set_top_friction]) 777 775 #print domain.quantities['friction'].get_values() 778 assert allclose(domain.quantities['friction'].get_values(),\779 [[ 0.09, 0.09, 0.09],780 [ 0.09, 0.09, 0.09],781 [ 0.07, 0.07, 0.07],782 [ 0.07, 0.07, 0.07],783 [ 1.0, 1.0, 1.0],784 [ 1.0, 1.0, 1.0]])776 assert num.allclose(domain.quantities['friction'].get_values(), 777 [[ 0.09, 0.09, 0.09], 778 [ 0.09, 0.09, 0.09], 779 [ 0.07, 0.07, 0.07], 780 [ 0.07, 0.07, 0.07], 781 [ 1.0, 1.0, 1.0], 782 [ 1.0, 1.0, 1.0]]) 785 783 786 784 domain.set_region([set_all_friction]) 787 785 #print domain.quantities['friction'].get_values() 788 assert allclose(domain.quantities['friction'].get_values(),789 [[ 10.09, 10.09, 10.09],790 [ 10.09, 10.09, 10.09],791 [ 10.07, 10.07, 10.07],792 [ 10.07, 10.07, 10.07],793 [ 11.0, 11.0, 11.0],794 [ 11.0, 11.0, 11.0]])786 assert num.allclose(domain.quantities['friction'].get_values(), 787 [[ 10.09, 10.09, 10.09], 788 [ 10.09, 10.09, 10.09], 789 [ 10.07, 10.07, 10.07], 790 [ 10.07, 10.07, 10.07], 791 [ 11.0, 11.0, 11.0], 792 [ 11.0, 11.0, 11.0]]) 795 793 796 794 #------------------------------------------------------------- -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_ermapper.py
r5897 r6145 5 5 6 6 import ermapper_grids 7 import Numeric8 7 from os import remove 8 9 import Numeric as num 10 9 11 10 12 class Test_ERMapper(unittest.TestCase): … … 19 21 data_filename = 'test_write_ermapper_grid' 20 22 21 original_grid = Numeric.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]])23 original_grid = num.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]]) 22 24 23 25 # Check that the function works when passing the filename without … … 25 27 ermapper_grids.write_ermapper_grid(data_filename, original_grid) 26 28 new_grid = ermapper_grids.read_ermapper_grid(data_filename) 27 assert Numeric.allclose(original_grid, new_grid)29 assert num.allclose(original_grid, new_grid) 28 30 29 31 # Check that the function works when passing the filename with … … 31 33 ermapper_grids.write_ermapper_grid(header_filename, original_grid) 32 34 new_grid = ermapper_grids.read_ermapper_grid(header_filename) 33 assert Numeric.allclose(original_grid, new_grid)35 assert num.allclose(original_grid, new_grid) 34 36 35 37 # Clean up created files … … 40 42 # Setup test data 41 43 filename = 'test_write_ermapper_grid' 42 original_grid = Numeric.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])44 original_grid = num.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0]) 43 45 44 46 # Write test data 45 ermapper_grids.write_ermapper_data(original_grid, filename, Numeric.Float64)47 ermapper_grids.write_ermapper_data(original_grid, filename, num.Float64) 46 48 47 49 # Read in the test data 48 new_grid = ermapper_grids.read_ermapper_data(filename, Numeric.Float64)50 new_grid = ermapper_grids.read_ermapper_data(filename, num.Float64) 49 51 50 52 # Check that the test data that has been read in matches the original data 51 assert Numeric.allclose(original_grid, new_grid)53 assert num.allclose(original_grid, new_grid) 52 54 53 55 # Clean up created files … … 57 59 # Setup test data 58 60 filename = 'test_write_ermapper_grid' 59 original_grid = Numeric.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])61 original_grid = num.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0]) 60 62 61 63 # Write test data … … 66 68 67 69 # Check that the test data that has been read in matches the original data 68 assert Numeric.allclose(original_grid, new_grid)70 assert num.allclose(original_grid, new_grid) 69 71 70 72 # Clean up created files … … 75 77 76 78 # setup test data 77 original_grid = Numeric.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]])79 original_grid = num.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]]) 78 80 # Write test data 79 81 ermapper_grids.write_ermapper_data(original_grid, data_filename) -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py
r5897 r6145 7 7 8 8 from anuga.config import epsilon 9 from Numeric import allclose, array, ones, Float10 9 from general_mesh import General_mesh 11 10 from anuga.coordinate_transforms.geo_reference import Geo_reference 12 11 12 import Numeric as num 13 13 14 14 … … 23 23 def test_get_vertex_coordinates(self): 24 24 from mesh_factory import rectangular 25 from Numeric import zeros, Float26 25 27 26 #Create basic mesh … … 30 29 31 30 32 assert allclose(domain.get_nodes(), nodes)31 assert num.allclose(domain.get_nodes(), nodes) 33 32 34 33 … … 41 40 for j in range(3): 42 41 k = triangles[i,j] #Index of vertex j in triangle i 43 assert allclose(V[3*i+j,:], nodes[k])42 assert num.allclose(V[3*i+j,:], nodes[k]) 44 43 45 44 def test_get_vertex_coordinates_with_geo_ref(self): … … 55 54 f = [4.0, 0.0] 56 55 57 nodes = array([a, b, c, d, e, f])56 nodes = num.array([a, b, c, d, e, f]) 58 57 59 58 nodes_absolute = geo.get_absolute(nodes) 60 59 61 60 #bac, bce, ecf, dbe, daf, dae 62 triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])61 triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 63 62 64 63 domain = General_mesh(nodes, triangles, 65 64 geo_reference = geo) 66 65 verts = domain.get_vertex_coordinates(triangle_id=0) 67 self.assert_( allclose(array([b,a,c]), verts))66 self.assert_(num.allclose(num.array([b,a,c]), verts)) 68 67 verts = domain.get_vertex_coordinates(triangle_id=0) 69 self.assert_( allclose(array([b,a,c]), verts))68 self.assert_(num.allclose(num.array([b,a,c]), verts)) 70 69 verts = domain.get_vertex_coordinates(triangle_id=0, 71 70 absolute=True) 72 self.assert_( allclose(array([nodes_absolute[1],71 self.assert_(num.allclose(num.array([nodes_absolute[1], 73 72 nodes_absolute[0], 74 73 nodes_absolute[2]]), verts)) 75 74 verts = domain.get_vertex_coordinates(triangle_id=0, 76 75 absolute=True) 77 self.assert_( allclose(array([nodes_absolute[1],76 self.assert_(num.allclose(num.array([nodes_absolute[1], 78 77 nodes_absolute[0], 79 78 nodes_absolute[2]]), verts)) … … 86 85 """ 87 86 from mesh_factory import rectangular 88 from Numeric import zeros, Float89 87 90 88 #Create basic mesh … … 93 91 94 92 95 assert allclose(domain.get_nodes(), nodes)93 assert num.allclose(domain.get_nodes(), nodes) 96 94 97 95 … … 104 102 for j in range(3): 105 103 k = triangles[i,j] #Index of vertex j in triangle i 106 assert allclose(V[j,:], nodes[k])104 assert num.allclose(V[j,:], nodes[k]) 107 105 108 106 … … 114 112 """ 115 113 from mesh_factory import rectangular 116 from Numeric import zeros, Float117 114 118 115 #Create basic mesh … … 121 118 122 119 value = [7] 123 assert allclose(domain.get_triangles(), triangles)124 assert allclose(domain.get_triangles([0,4]),120 assert num.allclose(domain.get_triangles(), triangles) 121 assert num.allclose(domain.get_triangles([0,4]), 125 122 [triangles[0], triangles[4]]) 126 123 … … 130 127 """ 131 128 from mesh_factory import rectangular 132 from Numeric import zeros, Float, array 133 134 a = [0.0, 0.0] 135 b = [0.0, 2.0] 136 c = [2.0, 0.0] 137 d = [0.0, 4.0] 138 e = [2.0, 2.0] 139 f = [4.0, 0.0] 140 141 nodes = array([a, b, c, d, e, f]) 129 130 a = [0.0, 0.0] 131 b = [0.0, 2.0] 132 c = [2.0, 0.0] 133 d = [0.0, 4.0] 134 e = [2.0, 2.0] 135 f = [4.0, 0.0] 136 137 nodes = num.array([a, b, c, d, e, f]) 142 138 #bac, bce, ecf, dbe, daf, dae 143 triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])139 triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 144 140 145 141 domain1 = General_mesh(nodes, triangles) … … 162 158 #print count 163 159 # 164 assert allclose(count, domain.number_of_triangles_per_node)160 assert num.allclose(count, domain.number_of_triangles_per_node) 165 161 166 162 # Check indices … … 196 192 f = [4.0, 0.0] 197 193 198 nodes = array([a, b, c, d, e, f])194 nodes = num.array([a, b, c, d, e, f]) 199 195 #bac, bce, ecf, dbe, daf, dae 200 triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])196 triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 201 197 202 198 domain = General_mesh(nodes, triangles) … … 204 200 # One node 205 201 L = domain.get_triangles_and_vertices_per_node(node=2) 206 assert allclose(L[0], [0, 2])207 assert allclose(L[1], [1, 1])208 assert allclose(L[2], [2, 1])202 assert num.allclose(L[0], [0, 2]) 203 assert num.allclose(L[1], [1, 1]) 204 assert num.allclose(L[2], [2, 1]) 209 205 210 206 # All nodes … … 213 209 for i, Lref in enumerate(ALL): 214 210 L = domain.get_triangles_and_vertices_per_node(node=i) 215 assert allclose(L, Lref)211 assert num.allclose(L, Lref) 216 212 217 213 … … 224 220 from mesh_factory import rectangular 225 221 from shallow_water import Domain 226 from Numeric import zeros, Float227 222 228 223 #Create basic mesh … … 239 234 from mesh_factory import rectangular 240 235 from shallow_water import Domain 241 from Numeric import zeros, Float242 236 243 237 #Create basic mesh … … 273 267 f = [4.0, 0.0] 274 268 275 nodes = array([a, b, c, d, e, f])269 nodes = num.array([a, b, c, d, e, f]) 276 270 277 271 nodes_absolute = geo.get_absolute(nodes) 278 272 279 273 #bac, bce, ecf, dbe, daf, dae 280 triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])274 triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 281 275 282 276 domain = General_mesh(nodes, triangles, … … 310 304 f = [4.0, 0.0] 311 305 312 nodes = array([a, b, c, d, e, f])306 nodes = num.array([a, b, c, d, e, f]) 313 307 314 308 nodes_absolute = geo.get_absolute(nodes) 315 309 316 310 # max index is 5, use 5, expect success 317 triangles = array([[1,5,2], [1,2,4], [4,2,5], [3,1,4]])311 triangles = num.array([[1,5,2], [1,2,4], [4,2,5], [3,1,4]]) 318 312 General_mesh(nodes, triangles, geo_reference=geo) 319 313 320 314 # max index is 5, use 6, expect assert failure 321 triangles = array([[1,6,2], [1,2,4], [4,2,5], [3,1,4]])315 triangles = num.array([[1,6,2], [1,2,4], [4,2,5], [3,1,4]]) 322 316 self.failUnlessRaises(AssertionError, General_mesh, 323 317 nodes, triangles, geo_reference=geo) 324 318 325 319 # max index is 5, use 10, expect assert failure 326 triangles = array([[1,10,2], [1,2,4], [4,2,5], [3,1,4]])320 triangles = num.array([[1,10,2], [1,2,4], [4,2,5], [3,1,4]]) 327 321 self.failUnlessRaises(AssertionError, General_mesh, 328 322 nodes, triangles, geo_reference=geo) -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py
r5897 r6145 6 6 from generic_boundary_conditions import * 7 7 from anuga.config import epsilon 8 from Numeric import allclose, array 8 9 import Numeric as num 9 10 10 11 … … 44 45 45 46 q = Bd.evaluate() 46 assert allclose(q, x)47 assert num.allclose(q, x) 47 48 48 49 … … 98 99 q = T.evaluate(0, 2) #Vol=0, edge=2 99 100 100 assert allclose(q, [1.5, 2.5])101 assert num.allclose(q, [1.5, 2.5]) 101 102 102 103 … … 175 176 176 177 #Check that midpoint coordinates at boundary are correctly computed 177 assert allclose( F.midpoint_coordinates,178 [[1.0, 0.0], [0.0, 1.0], [3.0, 0.0],179 [3.0, 1.0], [1.0, 3.0], [0.0, 3.0]])178 assert num.allclose( F.midpoint_coordinates, 179 [[1.0, 0.0], [0.0, 1.0], [3.0, 0.0], 180 [3.0, 1.0], [1.0, 3.0], [0.0, 3.0]]) 180 181 181 182 #assert allclose(F.midpoint_coordinates[(3,2)], [0.0, 3.0]) … … 193 194 domain.time = 5*30/2 #A quarter way through first step 194 195 q = F.evaluate() 195 assert allclose(q, [1.0/4, sin(2*pi/10)/4])196 assert num.allclose(q, [1.0/4, sin(2*pi/10)/4]) 196 197 197 198 198 199 domain.time = 2.5*5*60 #Half way between steps 2 and 3 199 200 q = F.evaluate() 200 assert allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2])201 assert num.allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2]) 201 202 202 203 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_ghost.py
r5897 r6145 6 6 from domain import * 7 7 from anuga.config import epsilon 8 from Numeric import allclose, array, ones, Float9 10 11 8 12 9 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py
r6003 r6145 13 13 from mesh_factory import rectangular 14 14 from anuga.config import epsilon 15 from Numeric import allclose, array, Int16 15 17 16 from anuga.coordinate_transforms.geo_reference import Geo_reference … … 19 18 from anuga.utilities.numerical_tools import ensure_numeric 20 19 20 import Numeric as num 21 22 21 23 def distance(x, y): 22 return sqrt( sum( ( array(x)-array(y))**2 ))24 return sqrt( sum( (num.array(x)-num.array(y))**2 )) 23 25 24 26 class Test_Mesh(unittest.TestCase): … … 63 65 #Normals 64 66 normals = mesh.get_normals() 65 assert allclose(normals[0, 0:2], [3.0/5, 4.0/5])66 assert allclose(normals[0, 2:4], [-1.0, 0.0])67 assert allclose(normals[0, 4:6], [0.0, -1.0])68 69 assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])70 assert allclose(mesh.get_normal(0,1), [-1.0, 0.0])71 assert allclose(mesh.get_normal(0,2), [0.0, -1.0])67 assert num.allclose(normals[0, 0:2], [3.0/5, 4.0/5]) 68 assert num.allclose(normals[0, 2:4], [-1.0, 0.0]) 69 assert num.allclose(normals[0, 4:6], [0.0, -1.0]) 70 71 assert num.allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5]) 72 assert num.allclose(mesh.get_normal(0,1), [-1.0, 0.0]) 73 assert num.allclose(mesh.get_normal(0,2), [0.0, -1.0]) 72 74 73 75 #Edge lengths 74 assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])76 assert num.allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0]) 75 77 76 78 … … 81 83 82 84 V = mesh.get_vertex_coordinates() 83 assert allclose(V, [ [0.0, 0.0],85 assert num.allclose(V, [ [0.0, 0.0], 84 86 [4.0, 0.0], 85 87 [0.0, 3.0] ]) 86 88 87 89 V0 = mesh.get_vertex_coordinate(0, 0) 88 assert allclose(V0, [0.0, 0.0])90 assert num.allclose(V0, [0.0, 0.0]) 89 91 90 92 V1 = mesh.get_vertex_coordinate(0, 1) 91 assert allclose(V1, [4.0, 0.0])93 assert num.allclose(V1, [4.0, 0.0]) 92 94 93 95 V2 = mesh.get_vertex_coordinate(0, 2) 94 assert allclose(V2, [0.0, 3.0])96 assert num.allclose(V2, [0.0, 3.0]) 95 97 96 98 … … 237 239 238 240 mesh = Mesh(points, vertices,use_inscribed_circle=False) 239 assert allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work'241 assert num.allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work' 240 242 241 243 mesh = Mesh(points, vertices,use_inscribed_circle=True) 242 assert allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work'244 assert num.allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work' 243 245 244 246 def test_inscribed_circle_rightangle_triangle(self): … … 252 254 253 255 mesh = Mesh(points, vertices,use_inscribed_circle=False) 254 assert allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work'256 assert num.allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work' 255 257 256 258 mesh = Mesh(points, vertices,use_inscribed_circle=True) 257 assert allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work'259 assert num.allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work' 258 260 259 261 … … 269 271 assert mesh.areas[0] == 2.0 270 272 271 assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])273 assert num.allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3]) 272 274 273 275 … … 303 305 assert mesh.edgelengths[1,2] == sqrt(8.0) 304 306 305 assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])306 assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])307 assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])308 assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])307 assert num.allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3]) 308 assert num.allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3]) 309 assert num.allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3]) 310 assert num.allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3]) 309 311 310 312 def test_mesh_and_neighbours(self): … … 698 700 """ 699 701 from mesh_factory import rectangular 700 from Numeric import zeros, Float701 702 702 703 #Create basic mesh … … 717 718 from mesh_factory import rectangular 718 719 #from mesh import Mesh 719 from Numeric import zeros, Float720 720 721 721 #Create basic mesh … … 727 727 728 728 assert len(P) == 8 729 assert allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],730 [1.0, 0.5], [1.0, 1.0], [0.5, 1.0],731 [0.0, 1.0], [0.0, 0.5]])729 assert num.allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0], 730 [1.0, 0.5], [1.0, 1.0], [0.5, 1.0], 731 [0.0, 1.0], [0.0, 0.5]]) 732 732 for p in points: 733 733 #print p, P … … 736 736 737 737 def test_boundary_polygon_II(self): 738 from Numeric import zeros, Float739 740 738 741 739 #Points … … 764 762 765 763 assert len(P) == 8 766 assert allclose(P, [a, d, f, i, h, e, c, b])764 assert num.allclose(P, [a, d, f, i, h, e, c, b]) 767 765 768 766 for p in points: … … 774 772 """Same as II but vertices ordered differently 775 773 """ 776 777 from Numeric import zeros, Float778 774 779 775 … … 804 800 805 801 assert len(P) == 8 806 assert allclose(P, [a, d, f, i, h, e, c, b])802 assert num.allclose(P, [a, d, f, i, h, e, c, b]) 807 803 808 804 for p in points: … … 815 811 is partitioned using pymetis. 816 812 """ 817 818 from Numeric import zeros, Float819 813 820 814 … … 847 841 848 842 # Note that point e appears twice! 849 assert allclose(P, [a, d, f, e, g, h, e, c, b])843 assert num.allclose(P, [a, d, f, e, g, h, e, c, b]) 850 844 851 845 for p in points: … … 864 858 """ 865 859 866 from Numeric import zeros, Float867 860 from mesh_factory import rectangular 868 861 … … 907 900 908 901 """ 909 from Numeric import zeros, Float910 911 902 912 903 #Points … … 955 946 956 947 assert len(P) == 8 957 assert allclose(P, [a, d, f, i, h, e, c, b])958 assert allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)])948 assert num.allclose(P, [a, d, f, i, h, e, c, b]) 949 assert num.allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)]) 959 950 960 951 … … 1101 1092 [ 35406.3359375 , 79332.9140625 ]] 1102 1093 1103 scaled_points = ensure_numeric(points, Int)/1000 # Simplify for ease of interpretation1094 scaled_points = ensure_numeric(points, num.Int)/1000 # Simplify for ease of interpretation 1104 1095 1105 1096 triangles = [[ 0, 1, 2], … … 1142 1133 assert is_inside_polygon(p, P) 1143 1134 1144 assert allclose(P, Pref)1135 assert num.allclose(P, Pref) 1145 1136 1146 1137 def test_lone_vertices(self): … … 1199 1190 boundary_polygon = mesh.get_boundary_polygon() 1200 1191 1201 assert allclose(absolute_points, boundary_polygon)1192 assert num.allclose(absolute_points, boundary_polygon) 1202 1193 1203 1194 def test_get_triangle_containing_point(self): … … 1246 1237 1247 1238 neighbours = mesh.get_triangle_neighbours(0) 1248 assert allclose(neighbours, [-1,1,-1])1239 assert num.allclose(neighbours, [-1,1,-1]) 1249 1240 neighbours = mesh.get_triangle_neighbours(-10) 1250 1241 assert neighbours == [] … … 1295 1286 for x in L: 1296 1287 if x.triangle_id % 2 == 0: 1297 assert allclose(x.length, ceiling-y_line)1288 assert num.allclose(x.length, ceiling-y_line) 1298 1289 else: 1299 assert allclose(x.length, y_line-floor)1290 assert num.allclose(x.length, y_line-floor) 1300 1291 1301 1292 1302 assert allclose(x.normal, [0,-1])1303 1304 assert allclose(x.segment[1][0], x.segment[0][0] + x.length)1305 assert allclose(x.segment[0][1], y_line)1306 assert allclose(x.segment[1][1], y_line)1293 assert num.allclose(x.normal, [0,-1]) 1294 1295 assert num.allclose(x.segment[1][0], x.segment[0][0] + x.length) 1296 assert num.allclose(x.segment[0][1], y_line) 1297 assert num.allclose(x.segment[1][1], y_line) 1307 1298 1308 1299 assert x.triangle_id in intersected_triangles … … 1311 1302 1312 1303 msg = 'Segments do not add up' 1313 assert allclose(total_length, 2), msg1304 assert num.allclose(total_length, 2), msg 1314 1305 1315 1306 … … 1346 1337 total_length = 0 1347 1338 for x in L: 1348 assert allclose(x.length, 1.0)1349 assert allclose(x.normal, [0,-1])1350 1351 assert allclose(x.segment[1][0], x.segment[0][0] + x.length)1352 assert allclose(x.segment[0][1], y_line)1353 assert allclose(x.segment[1][1], y_line)1339 assert num.allclose(x.length, 1.0) 1340 assert num.allclose(x.normal, [0,-1]) 1341 1342 assert num.allclose(x.segment[1][0], x.segment[0][0] + x.length) 1343 assert num.allclose(x.segment[0][1], y_line) 1344 assert num.allclose(x.segment[1][1], y_line) 1354 1345 1355 1346 … … 1360 1351 1361 1352 msg = 'Segments do not add up' 1362 assert allclose(total_length, 2), msg1353 assert num.allclose(total_length, 2), msg 1363 1354 1364 1355 … … 1400 1391 for x in L: 1401 1392 if x.triangle_id == 1: 1402 assert allclose(x.length, 1)1403 assert allclose(x.normal, [0, -1])1393 assert num.allclose(x.length, 1) 1394 assert num.allclose(x.normal, [0, -1]) 1404 1395 1405 1396 if x.triangle_id == 5: 1406 assert allclose(x.length, 0.5)1407 assert allclose(x.normal, [0, -1])1397 assert num.allclose(x.length, 0.5) 1398 assert num.allclose(x.normal, [0, -1]) 1408 1399 1409 1400 … … 1413 1404 1414 1405 msg = 'Segments do not add up' 1415 assert allclose(total_length, 1.5), msg1406 assert num.allclose(total_length, 1.5), msg 1416 1407 1417 1408 … … 1447 1438 total_length = 0 1448 1439 for i, x in enumerate(L): 1449 assert allclose(x.length, s2)1450 assert allclose(x.normal, [-s2, -s2])1451 assert allclose(sum(x.normal**2), 1)1440 assert num.allclose(x.length, s2) 1441 assert num.allclose(x.normal, [-s2, -s2]) 1442 assert num.allclose(sum(x.normal**2), 1) 1452 1443 1453 1444 assert x.triangle_id in intersected_triangles … … 1456 1447 1457 1448 msg = 'Segments do not add up' 1458 assert allclose(total_length, 4*s2), msg1449 assert num.allclose(total_length, 4*s2), msg 1459 1450 1460 1451 … … 1471 1462 total_length = 0 1472 1463 for i, x in enumerate(L): 1473 assert allclose(x.length, s2)1474 assert allclose(x.normal, [s2, s2])1475 assert allclose(sum(x.normal**2), 1)1464 assert num.allclose(x.length, s2) 1465 assert num.allclose(x.normal, [s2, s2]) 1466 assert num.allclose(sum(x.normal**2), 1) 1476 1467 1477 1468 assert x.triangle_id in intersected_triangles … … 1480 1471 1481 1472 msg = 'Segments do not add up' 1482 assert allclose(total_length, 4*s2), msg1473 assert num.allclose(total_length, 4*s2), msg 1483 1474 1484 1475 … … 1496 1487 total_length = 0 1497 1488 for i, x in enumerate(L): 1498 assert allclose(x.length, 2*s2)1499 assert allclose(x.normal, [-s2, s2])1500 assert allclose(sum(x.normal**2), 1)1489 assert num.allclose(x.length, 2*s2) 1490 assert num.allclose(x.normal, [-s2, s2]) 1491 assert num.allclose(sum(x.normal**2), 1) 1501 1492 1502 1493 assert x.triangle_id in intersected_triangles … … 1505 1496 1506 1497 msg = 'Segments do not add up' 1507 assert allclose(total_length, 4*s2), msg1498 assert num.allclose(total_length, 4*s2), msg 1508 1499 1509 1500 … … 1520 1511 total_length = 0 1521 1512 for i, x in enumerate(L): 1522 assert allclose(x.length, 2*s2)1523 assert allclose(x.normal, [s2, -s2])1524 assert allclose(sum(x.normal**2), 1)1513 assert num.allclose(x.length, 2*s2) 1514 assert num.allclose(x.normal, [s2, -s2]) 1515 assert num.allclose(sum(x.normal**2), 1) 1525 1516 1526 1517 assert x.triangle_id in intersected_triangles … … 1529 1520 1530 1521 msg = 'Segments do not add up' 1531 assert allclose(total_length, 4*s2), msg1522 assert num.allclose(total_length, 4*s2), msg 1532 1523 1533 1524 … … 1545 1536 total_length = 0 1546 1537 for i, x in enumerate(L): 1547 assert allclose(x.length, s2)1548 assert allclose(x.normal, [-s2, -s2])1549 assert allclose(sum(x.normal**2), 1)1538 assert num.allclose(x.length, s2) 1539 assert num.allclose(x.normal, [-s2, -s2]) 1540 assert num.allclose(sum(x.normal**2), 1) 1550 1541 1551 1542 assert x.triangle_id in intersected_triangles … … 1554 1545 1555 1546 msg = 'Segments do not add up' 1556 assert allclose(total_length, 2*s2), msg1547 assert num.allclose(total_length, 2*s2), msg 1557 1548 1558 1549 … … 1567 1558 total_length = 0 1568 1559 for i, x in enumerate(L): 1569 assert allclose(x.normal, [-s2, -s2])1570 assert allclose(sum(x.normal**2), 1)1560 assert num.allclose(x.normal, [-s2, -s2]) 1561 assert num.allclose(sum(x.normal**2), 1) 1571 1562 1572 1563 msg = 'Triangle %d' %x.triangle_id + ' is not in %s' %(intersected_triangles) … … 1603 1594 assert len(L) == 1 1604 1595 assert L[0].triangle_id == 3 1605 assert allclose(L[0].length, 0.5)1606 assert allclose(L[0].normal, [-1,0])1596 assert num.allclose(L[0].length, 0.5) 1597 assert num.allclose(L[0].normal, [-1,0]) 1607 1598 1608 1599 … … 1615 1606 assert len(L) == 1 1616 1607 assert L[0].triangle_id == 3 1617 assert allclose(L[0].length, 0.4)1618 assert allclose(L[0].normal, [-1,0])1608 assert num.allclose(L[0].length, 0.4) 1609 assert num.allclose(L[0].normal, [-1,0]) 1619 1610 1620 1611 intersected_triangles = [3] … … 1628 1619 assert len(L) == 1 1629 1620 assert L[0].triangle_id == 3 1630 assert allclose(L[0].length, 0.4)1631 assert allclose(L[0].normal, [1,0])1621 assert num.allclose(L[0].length, 0.4) 1622 assert num.allclose(L[0].normal, [1,0]) 1632 1623 1633 1624 … … 1663 1654 for x in L: 1664 1655 if x.triangle_id == 3: 1665 assert allclose(x.length, 0.5)1666 assert allclose(x.normal, [-1,0])1656 assert num.allclose(x.length, 0.5) 1657 assert num.allclose(x.normal, [-1,0]) 1667 1658 1668 1659 if x.triangle_id == 2: 1669 assert allclose(x.length, s2)1670 assert allclose(x.normal, [-s2,-s2])1660 assert num.allclose(x.length, s2) 1661 assert num.allclose(x.normal, [-s2,-s2]) 1671 1662 1672 1663 … … 1701 1692 for x in L: 1702 1693 if x.triangle_id == 3: 1703 assert allclose(x.length, 0.5)1704 assert allclose(x.normal, [-1,0])1694 assert num.allclose(x.length, 0.5) 1695 assert num.allclose(x.normal, [-1,0]) 1705 1696 1706 1697 if x.triangle_id == 2: 1707 1698 msg = str(x.length) 1708 assert allclose(x.length, s2), msg1709 assert allclose(x.normal, [-s2,-s2])1699 assert num.allclose(x.length, s2), msg 1700 assert num.allclose(x.normal, [-s2,-s2]) 1710 1701 1711 1702 if x.triangle_id == 5: 1712 segvec = array([line[2][0]-1,1713 line[2][1]-1])1703 segvec = num.array([line[2][0]-1, 1704 line[2][1]-1]) 1714 1705 msg = str(x.length) 1715 assert allclose(x.length, sqrt(sum(segvec**2))), msg1716 assert allclose(x.normal, [-s2,-s2])1706 assert num.allclose(x.length, sqrt(sum(segvec**2))), msg 1707 assert num.allclose(x.normal, [-s2,-s2]) 1717 1708 1718 1709 … … 1752 1743 for x in L: 1753 1744 if x.triangle_id == 3: 1754 assert allclose(x.length, 0.5)1755 assert allclose(x.normal, [-1,0])1745 assert num.allclose(x.length, 0.5) 1746 assert num.allclose(x.normal, [-1,0]) 1756 1747 1757 1748 if x.triangle_id == 2: 1758 1749 msg = str(x.length) 1759 assert allclose(x.length, s2), msg1760 assert allclose(x.normal, [-s2,-s2])1750 assert num.allclose(x.length, s2), msg 1751 assert num.allclose(x.normal, [-s2,-s2]) 1761 1752 1762 1753 if x.triangle_id == 5: 1763 1754 if x.segment == ((1.0, 1.0), (1.25, 0.75)): 1764 segvec = array([line[2][0]-1,1765 line[2][1]-1])1755 segvec = num.array([line[2][0]-1, 1756 line[2][1]-1]) 1766 1757 msg = str(x.length) 1767 assert allclose(x.length, sqrt(sum(segvec**2))), msg1768 assert allclose(x.normal, [-s2,-s2])1758 assert num.allclose(x.length, sqrt(sum(segvec**2))), msg 1759 assert num.allclose(x.normal, [-s2,-s2]) 1769 1760 elif x.segment == ((1.25, 0.75), (1.5, 1.0)): 1770 segvec = array([1.5-line[2][0],1771 1.0-line[2][1]])1761 segvec = num.array([1.5-line[2][0], 1762 1.0-line[2][1]]) 1772 1763 1773 assert allclose(x.length, sqrt(sum(segvec**2))), msg1774 assert allclose(x.normal, [s2,-s2])1764 assert num.allclose(x.length, sqrt(sum(segvec**2))), msg 1765 assert num.allclose(x.normal, [s2,-s2]) 1775 1766 else: 1776 1767 msg = 'Unknow segment: %s' %x.segment … … 1780 1771 1781 1772 if x.triangle_id == 6: 1782 assert allclose(x.normal, [s2,-s2])1783 assert allclose(x.segment, ((1.5, 1.0), (2, 1.5)))1773 assert num.allclose(x.normal, [s2,-s2]) 1774 assert num.allclose(x.segment, ((1.5, 1.0), (2, 1.5))) 1784 1775 1785 1776 … … 1791 1782 #xi1 = line[1][0] 1792 1783 #eta1 = line[1][1] 1793 #linevector = array([xi1-xi0, eta1-eta0])1784 #linevector = num.array([xi1-xi0, eta1-eta0]) 1794 1785 #linelength = sqrt(sum(linevector**2)) 1795 1786 # … … 1847 1838 ref_length = line[1][1] - line[0][1] 1848 1839 #print ref_length, total_length 1849 assert allclose(total_length, ref_length)1840 assert num.allclose(total_length, ref_length) 1850 1841 1851 1842 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py
r5897 r6145 3 3 4 4 import unittest 5 6 from Numeric import allclose, array7 8 5 9 6 #from anuga.pyvolution.pmesh2domain import * … … 16 13 from anuga.coordinate_transforms.geo_reference import Geo_reference 17 14 from anuga.pmesh.mesh import importMeshFromFile 15 16 import Numeric as num 17 18 18 19 19 class Test_pmesh2domain(unittest.TestCase): … … 66 66 67 67 tags = {} 68 b1 = Dirichlet_boundary(conserved_quantities = array([0.0]))69 b2 = Dirichlet_boundary(conserved_quantities = array([1.0]))70 b3 = Dirichlet_boundary(conserved_quantities = array([2.0]))68 b1 = Dirichlet_boundary(conserved_quantities = num.array([0.0])) 69 b2 = Dirichlet_boundary(conserved_quantities = num.array([1.0])) 70 b3 = Dirichlet_boundary(conserved_quantities = num.array([2.0])) 71 71 tags["1"] = b1 72 72 tags["2"] = b2 … … 80 80 answer = [[0., 8., 0.], 81 81 [0., 10., 8.]] 82 assert allclose(domain.quantities['elevation'].vertex_values,83 answer)82 assert num.allclose(domain.quantities['elevation'].vertex_values, 83 answer) 84 84 85 85 #print domain.quantities['stage'].vertex_values 86 86 answer = [[0., 12., 10.], 87 87 [0., 10., 12.]] 88 assert allclose(domain.quantities['stage'].vertex_values,89 answer)88 assert num.allclose(domain.quantities['stage'].vertex_values, 89 answer) 90 90 91 91 #print domain.quantities['friction'].vertex_values 92 92 answer = [[0.01, 0.04, 0.03], 93 93 [0.01, 0.02, 0.04]] 94 assert allclose(domain.quantities['friction'].vertex_values,95 answer)96 97 #print domain.quantities['friction'].vertex_values 98 assert allclose(domain.tagged_elements['dsg'][0],0)99 assert allclose(domain.tagged_elements['ole nielsen'][0],1)94 assert num.allclose(domain.quantities['friction'].vertex_values, 95 answer) 96 97 #print domain.quantities['friction'].vertex_values 98 assert num.allclose(domain.tagged_elements['dsg'][0],0) 99 assert num.allclose(domain.tagged_elements['ole nielsen'][0],1) 100 100 101 101 self.failUnless( domain.boundary[(1, 0)] == '1', … … 159 159 160 160 tags = {} 161 b1 = Dirichlet_boundary(conserved_quantities = array([0.0]))162 b2 = Dirichlet_boundary(conserved_quantities = array([1.0]))163 b3 = Dirichlet_boundary(conserved_quantities = array([2.0]))161 b1 = Dirichlet_boundary(conserved_quantities = num.array([0.0])) 162 b2 = Dirichlet_boundary(conserved_quantities = num.array([1.0])) 163 b3 = Dirichlet_boundary(conserved_quantities = num.array([2.0])) 164 164 tags["1"] = b1 165 165 tags["2"] = b2 … … 174 174 answer = [[0., 8., 0.], 175 175 [0., 10., 8.]] 176 assert allclose(domain.quantities['elevation'].vertex_values,177 answer)176 assert num.allclose(domain.quantities['elevation'].vertex_values, 177 answer) 178 178 179 179 #print domain.quantities['stage'].vertex_values 180 180 answer = [[0., 12., 10.], 181 181 [0., 10., 12.]] 182 assert allclose(domain.quantities['stage'].vertex_values,183 answer)182 assert num.allclose(domain.quantities['stage'].vertex_values, 183 answer) 184 184 185 185 #print domain.quantities['friction'].vertex_values 186 186 answer = [[0.01, 0.04, 0.03], 187 187 [0.01, 0.02, 0.04]] 188 assert allclose(domain.quantities['friction'].vertex_values,189 answer)190 191 #print domain.quantities['friction'].vertex_values 192 assert allclose(domain.tagged_elements['dsg'][0],0)193 assert allclose(domain.tagged_elements['ole nielsen'][0],1)188 assert num.allclose(domain.quantities['friction'].vertex_values, 189 answer) 190 191 #print domain.quantities['friction'].vertex_values 192 assert num.allclose(domain.tagged_elements['dsg'][0],0) 193 assert num.allclose(domain.tagged_elements['ole nielsen'][0],1) 194 194 195 195 self.failUnless( domain.boundary[(1, 0)] == '1', … … 228 228 #domain.set_tag_dict(tag_dict) 229 229 #Boundary tests 230 b1 = Dirichlet_boundary(conserved_quantities = array([0.0]))231 b2 = Dirichlet_boundary(conserved_quantities = array([1.0]))232 b3 = Dirichlet_boundary(conserved_quantities = array([1.0]))230 b1 = Dirichlet_boundary(conserved_quantities = num.array([0.0])) 231 b2 = Dirichlet_boundary(conserved_quantities = num.array([1.0])) 232 b3 = Dirichlet_boundary(conserved_quantities = num.array([1.0])) 233 233 #test adding a boundary 234 234 tags = {} -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
r6131 r6145 7 7 from quantity import * 8 8 from anuga.config import epsilon 9 from Numeric import allclose, array, ones, Float10 9 11 10 from anuga.fit_interpolate.fit import fit_to_mesh … … 16 15 from anuga.utilities.polygon import * 17 16 17 import Numeric as num 18 19 18 20 #Aux for fit_interpolate.fit example 19 21 def linear_function(point): 20 point = array(point)22 point = num.array(point) 21 23 return point[:,0]+point[:,1] 22 24 … … 63 65 64 66 quantity = Quantity(self.mesh1, [[1,2,3]]) 65 assert allclose(quantity.vertex_values, [[1.,2.,3.]])67 assert num.allclose(quantity.vertex_values, [[1.,2.,3.]]) 66 68 67 69 try: … … 84 86 85 87 quantity = Quantity(self.mesh1) 86 assert allclose(quantity.vertex_values, [[0.,0.,0.]])87 88 89 quantity = Quantity(self.mesh4) 90 assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],91 [0.,0.,0.], [0.,0.,0.]])88 assert num.allclose(quantity.vertex_values, [[0.,0.,0.]]) 89 90 91 quantity = Quantity(self.mesh4) 92 assert num.allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.], 93 [0.,0.,0.], [0.,0.,0.]]) 92 94 93 95 94 96 def test_interpolation(self): 95 97 quantity = Quantity(self.mesh1, [[1,2,3]]) 96 assert allclose(quantity.centroid_values, [2.0]) #Centroid97 98 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])98 assert num.allclose(quantity.centroid_values, [2.0]) #Centroid 99 100 assert num.allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 99 101 100 102 … … 102 104 quantity = Quantity(self.mesh4, 103 105 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 104 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid106 assert num.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 105 107 106 108 … … 108 110 109 111 #print quantity.vertex_values 110 assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],111 [3.+2./3, 6.+2./3, 4.+2./3],112 assert num.allclose(quantity.vertex_values, [[3.5, -1.0, 3.5], 113 [3.+2./3, 6.+2./3, 4.+2./3], 112 114 [4.6, 3.4, 1.], 113 115 [-5.0, 1.0, 4.0]]) 114 116 115 117 #print quantity.edge_values 116 assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25],117 [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],118 [2.2, 2.8, 4.0],119 [2.5, -0.5, -2.0]])118 assert num.allclose(quantity.edge_values, [[1.25, 3.5, 1.25], 119 [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6], 120 [2.2, 2.8, 4.0], 121 [2.5, -0.5, -2.0]]) 120 122 121 123 … … 128 130 quantity = Quantity(self.mesh4, 129 131 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 130 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids132 assert num.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids 131 133 132 134 v = quantity.get_maximum_value() … … 148 150 149 151 v = quantity.get_values(interpolation_points = [[x,y]]) 150 assert allclose(v, 5)152 assert num.allclose(v, 5) 151 153 152 154 153 155 x,y = quantity.get_minimum_location() 154 156 v = quantity.get_values(interpolation_points = [[x,y]]) 155 assert allclose(v, 0)157 assert num.allclose(v, 0) 156 158 157 159 … … 192 194 193 195 v = quantity.get_values(interpolation_points = [[x,y]]) 194 assert allclose(v, 6)196 assert num.allclose(v, 6) 195 197 196 198 x,y = quantity.get_minimum_location() 197 199 v = quantity.get_values(interpolation_points = [[x,y]]) 198 assert allclose(v, 2)200 assert num.allclose(v, 2) 199 201 200 202 #Multiple locations for maximum - 201 203 #Test that the algorithm picks the first occurrence 202 204 v = quantity.get_maximum_value(indices=[0,1,2]) 203 assert allclose(v, 4)205 assert num.allclose(v, 4) 204 206 205 207 i = quantity.get_maximum_index(indices=[0,1,2]) … … 212 214 213 215 v = quantity.get_values(interpolation_points = [[x,y]]) 214 assert allclose(v, 4)216 assert num.allclose(v, 4) 215 217 216 218 # More test of indices...... 217 219 v = quantity.get_maximum_value(indices=[2,3]) 218 assert allclose(v, 6)220 assert num.allclose(v, 6) 219 221 220 222 i = quantity.get_maximum_index(indices=[2,3]) … … 227 229 228 230 v = quantity.get_values(interpolation_points = [[x,y]]) 229 assert allclose(v, 6)231 assert num.allclose(v, 6) 230 232 231 233 … … 244 246 quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], 245 247 location = 'vertices') 246 assert allclose(quantity.vertex_values,247 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])248 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid249 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],250 [5., 5., 5.],251 [4.5, 4.5, 0.],252 [3.0, -1.5, -1.5]])248 assert num.allclose(quantity.vertex_values, 249 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 250 assert num.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 251 assert num.allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 252 [5., 5., 5.], 253 [4.5, 4.5, 0.], 254 [3.0, -1.5, -1.5]]) 253 255 254 256 255 257 # Test default 256 258 quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 257 assert allclose(quantity.vertex_values,258 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])259 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid260 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],261 [5., 5., 5.],262 [4.5, 4.5, 0.],263 [3.0, -1.5, -1.5]])259 assert num.allclose(quantity.vertex_values, 260 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 261 assert num.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 262 assert num.allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 263 [5., 5., 5.], 264 [4.5, 4.5, 0.], 265 [3.0, -1.5, -1.5]]) 264 266 265 267 # Test centroids 266 268 quantity.set_values([1,2,3,4], location = 'centroids') 267 assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid269 assert num.allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid 268 270 269 271 # Test exceptions … … 288 290 289 291 quantity.set_values(1.0, location = 'vertices') 290 assert allclose(quantity.vertex_values,291 [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])292 293 assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid294 assert allclose(quantity.edge_values, [[1, 1, 1],295 [1, 1, 1],296 [1, 1, 1],297 [1, 1, 1]])292 assert num.allclose(quantity.vertex_values, 293 [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]]) 294 295 assert num.allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid 296 assert num.allclose(quantity.edge_values, [[1, 1, 1], 297 [1, 1, 1], 298 [1, 1, 1], 299 [1, 1, 1]]) 298 300 299 301 300 302 quantity.set_values(2.0, location = 'centroids') 301 assert allclose(quantity.centroid_values, [2, 2, 2, 2])303 assert num.allclose(quantity.centroid_values, [2, 2, 2, 2]) 302 304 303 305 … … 310 312 quantity.set_values(f, location = 'vertices') 311 313 #print "quantity.vertex_values",quantity.vertex_values 312 assert allclose(quantity.vertex_values,313 [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])314 assert allclose(quantity.centroid_values,315 [4.0/3, 8.0/3, 10.0/3, 10.0/3])316 assert allclose(quantity.edge_values,317 [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])314 assert num.allclose(quantity.vertex_values, 315 [[2,0,2], [2,2,4], [4,2,4], [4,2,4]]) 316 assert num.allclose(quantity.centroid_values, 317 [4.0/3, 8.0/3, 10.0/3, 10.0/3]) 318 assert num.allclose(quantity.edge_values, 319 [[1,2,1], [3,3,2], [3,4,3], [3,4,3]]) 318 320 319 321 320 322 quantity.set_values(f, location = 'centroids') 321 assert allclose(quantity.centroid_values,322 [4.0/3, 8.0/3, 10.0/3, 10.0/3])323 assert num.allclose(quantity.centroid_values, 324 [4.0/3, 8.0/3, 10.0/3, 10.0/3]) 323 325 324 326 … … 331 333 #print 'Q', quantity.get_integral() 332 334 333 assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)335 assert num.allclose(quantity.get_integral(), self.mesh4.get_area() * const) 334 336 335 337 # Try with a linear function … … 342 344 ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 343 345 344 assert allclose (quantity.get_integral(), ref_integral)346 assert num.allclose (quantity.get_integral(), ref_integral) 345 347 346 348 … … 350 352 quantity.set_vertex_values([0,1,2,3,4,5]) 351 353 352 assert allclose(quantity.vertex_values,353 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])354 assert allclose(quantity.centroid_values,355 [1., 7./3, 11./3, 8./3]) #Centroid356 assert allclose(quantity.edge_values, [[1., 1.5, 0.5],357 [3., 2.5, 1.5],358 [3.5, 4.5, 3.],359 [2.5, 3.5, 2]])354 assert num.allclose(quantity.vertex_values, 355 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 356 assert num.allclose(quantity.centroid_values, 357 [1., 7./3, 11./3, 8./3]) #Centroid 358 assert num.allclose(quantity.edge_values, [[1., 1.5, 0.5], 359 [3., 2.5, 1.5], 360 [3.5, 4.5, 3.], 361 [2.5, 3.5, 2]]) 360 362 361 363 … … 365 367 quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5]) 366 368 367 assert allclose(quantity.vertex_values,368 [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])369 assert num.allclose(quantity.vertex_values, 370 [[1,0,20], [1,20,4], [4,20,50], [30,1,4]]) 369 371 370 372 … … 376 378 377 379 378 assert allclose(quantity.vertex_values,379 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])380 assert num.allclose(quantity.vertex_values, 381 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 380 382 381 383 #Centroid 382 assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])383 384 assert allclose(quantity.edge_values, [[1., 1.5, 0.5],385 [3., 2.5, 1.5],386 [3.5, 4.5, 3.],387 [2.5, 3.5, 2]])384 assert num.allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) 385 386 assert num.allclose(quantity.edge_values, [[1., 1.5, 0.5], 387 [3., 2.5, 1.5], 388 [3.5, 4.5, 3.], 389 [2.5, 3.5, 2]]) 388 390 389 391 … … 399 401 400 402 quantity.set_values([0,2,3,5], indices=[0,2,3,5]) 401 assert allclose(quantity.vertex_values,402 [[0,0,2], [0,2,0], [0,2,5], [3,0,0]])403 assert num.allclose(quantity.vertex_values, 404 [[0,0,2], [0,2,0], [0,2,5], [3,0,0]]) 403 405 404 406 … … 408 410 409 411 # Indices refer to triangle numbers 410 assert allclose(quantity.vertex_values,411 [[3.14,3.14,3.14], [0,0,0],412 [3.14,3.14,3.14], [0,0,0]])412 assert num.allclose(quantity.vertex_values, 413 [[3.14,3.14,3.14], [0,0,0], 414 [3.14,3.14,3.14], [0,0,0]]) 413 415 414 416 … … 419 421 quantity.set_values(3.14, polygon=polygon) 420 422 421 assert allclose(quantity.vertex_values,422 [[0,0,0], [0,0,0], [0,0,0],423 [3.14,3.14,3.14]])423 assert num.allclose(quantity.vertex_values, 424 [[0,0,0], [0,0,0], [0,0,0], 425 [3.14,3.14,3.14]]) 424 426 425 427 … … 429 431 quantity.set_values(0.0) 430 432 quantity.set_values(3.14, location='centroids', polygon=polygon) 431 assert allclose(quantity.vertex_values,432 [[0,0,0],433 [3.14,3.14,3.14],434 [3.14,3.14,3.14],435 [0,0,0]])433 assert num.allclose(quantity.vertex_values, 434 [[0,0,0], 435 [3.14,3.14,3.14], 436 [3.14,3.14,3.14], 437 [0,0,0]]) 436 438 437 439 … … 441 443 #print 'Here 2' 442 444 quantity.set_values(3.14, polygon=polygon) 443 assert allclose(quantity.vertex_values,444 [[0,0,0],445 [3.14,3.14,3.14],446 [3.14,3.14,3.14],447 [0,0,0]])445 assert num.allclose(quantity.vertex_values, 446 [[0,0,0], 447 [3.14,3.14,3.14], 448 [3.14,3.14,3.14], 449 [0,0,0]]) 448 450 449 451 … … 478 480 479 481 # Indices refer to triangle numbers here - not vertices (why?) 480 assert allclose(quantity.vertex_values,481 [[3.14,3.14,3.14], [0,0,0],482 [3.14,3.14,3.14], [0,0,0]])482 assert num.allclose(quantity.vertex_values, 483 [[3.14,3.14,3.14], [0,0,0], 484 [3.14,3.14,3.14], [0,0,0]]) 483 485 484 486 485 487 486 488 # Now try with polygon (pick points where y>2) 487 polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]])489 polygon = num.array([[0,2.1], [4,2.1], [4,7], [0,7]]) 488 490 polygon += [G.xllcorner, G.yllcorner] 489 491 … … 491 493 quantity.set_values(3.14, polygon=polygon, location='centroids') 492 494 493 assert allclose(quantity.vertex_values,494 [[0,0,0], [0,0,0], [0,0,0],495 [3.14,3.14,3.14]])495 assert num.allclose(quantity.vertex_values, 496 [[0,0,0], [0,0,0], [0,0,0], 497 [3.14,3.14,3.14]]) 496 498 497 499 498 500 # Another polygon (pick triangle 1 and 2 (rightmost triangles) 499 polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])501 polygon = num.array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]) 500 502 polygon += [G.xllcorner, G.yllcorner] 501 503 … … 503 505 quantity.set_values(3.14, polygon=polygon) 504 506 505 assert allclose(quantity.vertex_values,506 [[0,0,0],507 [3.14,3.14,3.14],508 [3.14,3.14,3.14],509 [0,0,0]])507 assert num.allclose(quantity.vertex_values, 508 [[0,0,0], 509 [3.14,3.14,3.14], 510 [3.14,3.14,3.14], 511 [0,0,0]]) 510 512 511 513 … … 540 542 answer = linear_function(quantity.domain.get_vertex_coordinates()) 541 543 #print quantity.vertex_values, answer 542 assert allclose(quantity.vertex_values.flat, answer)544 assert num.allclose(quantity.vertex_values.flat, answer) 543 545 544 546 … … 553 555 #print vertex_attributes 554 556 quantity.set_values(vertex_attributes) 555 assert allclose(quantity.vertex_values.flat, answer)557 assert num.allclose(quantity.vertex_values.flat, answer) 556 558 557 559 … … 593 595 alpha = 0) 594 596 595 assert allclose( ref, [0,5,5] )597 assert num.allclose( ref, [0,5,5] ) 596 598 597 599 … … 610 612 # data_georef = data_georef, 611 613 # alpha = 0) 612 assert allclose(quantity.vertex_values.flat, ref)614 assert num.allclose(quantity.vertex_values.flat, ref) 613 615 614 616 … … 621 623 622 624 quantity.set_values(geospatial_data = geo, alpha = 0) 623 assert allclose(quantity.vertex_values.flat, ref)625 assert num.allclose(quantity.vertex_values.flat, ref) 624 626 625 627 … … 670 672 671 673 672 assert allclose(quantity.vertex_values.flat, answer)674 assert num.allclose(quantity.vertex_values.flat, answer) 673 675 674 676 675 677 #Check that values can be set from file using default attribute 676 678 quantity.set_values(filename = ptsfile, alpha = 0) 677 assert allclose(quantity.vertex_values.flat, answer)679 assert num.allclose(quantity.vertex_values.flat, answer) 678 680 679 681 #Cleanup … … 730 732 #print self.mesh4.nodes 731 733 #print inside_polygon(self.mesh4.nodes, polygon) 732 assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)734 assert num.allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 733 735 734 736 #print quantity.domain.get_vertex_coordinates() … … 752 754 753 755 # Check vertices in polygon have been set 754 assert allclose(take(quantity.vertex_values.flat, indices),755 756 assert num.allclose(take(quantity.vertex_values.flat, indices), 757 answer) 756 758 757 759 # Check vertices outside polygon are zero 758 760 indices = outside_polygon(quantity.domain.get_vertex_coordinates(), 759 761 polygon) 760 assert allclose(take(quantity.vertex_values.flat, indices),761 762 assert num.allclose(take(quantity.vertex_values.flat, indices), 763 0.0) 762 764 763 765 #Cleanup … … 815 817 verbose=False) 816 818 answer = linear_function(quantity.domain.get_vertex_coordinates()) 817 assert allclose(quantity.vertex_values.flat, answer)819 assert num.allclose(quantity.vertex_values.flat, answer) 818 820 819 821 … … 821 823 quantity.set_values(filename=ptsfile, 822 824 alpha=0) 823 assert allclose(quantity.vertex_values.flat, answer)825 assert num.allclose(quantity.vertex_values.flat, answer) 824 826 825 827 # Check cache … … 869 871 #print "answer",answer 870 872 871 assert allclose(quantity.vertex_values.flat, answer)873 assert num.allclose(quantity.vertex_values.flat, answer) 872 874 873 875 874 876 #Check that values can be set from file using default attribute 875 877 quantity.set_values(filename=txt_file, alpha=0) 876 assert allclose(quantity.vertex_values.flat, answer)878 assert num.allclose(quantity.vertex_values.flat, answer) 877 879 878 880 #Cleanup … … 913 915 #print "answer",answer 914 916 915 assert allclose(quantity.vertex_values.flat, answer)917 assert num.allclose(quantity.vertex_values.flat, answer) 916 918 917 919 918 920 #Check that values can be set from file using default attribute 919 921 quantity.set_values(filename=txt_file, alpha=0) 920 assert allclose(quantity.vertex_values.flat, answer)922 assert num.allclose(quantity.vertex_values.flat, answer) 921 923 922 924 #Cleanup … … 959 961 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 960 962 #print "answer",answer 961 assert allclose(quantity.vertex_values.flat, answer)963 assert num.allclose(quantity.vertex_values.flat, answer) 962 964 963 965 #Check that values can be set from file … … 967 969 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 968 970 #print "answer",answer 969 assert allclose(quantity.vertex_values.flat, answer)971 assert num.allclose(quantity.vertex_values.flat, answer) 970 972 971 973 972 974 #Check that values can be set from file using default attribute 973 975 quantity.set_values(filename=txt_file, alpha=0) 974 assert allclose(quantity.vertex_values.flat, answer)976 assert num.allclose(quantity.vertex_values.flat, answer) 975 977 976 978 #Cleanup … … 1041 1043 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 1042 1044 #print "answer",answer 1043 assert allclose(quantity.vertex_values.flat, answer)1045 assert num.allclose(quantity.vertex_values.flat, answer) 1044 1046 1045 1047 #Check that values can be set from file … … 1049 1051 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 1050 1052 #print "answer",answer 1051 assert allclose(quantity.vertex_values.flat, answer)1053 assert num.allclose(quantity.vertex_values.flat, answer) 1052 1054 1053 1055 1054 1056 #Check that values can be set from file using default attribute 1055 1057 quantity.set_values(filename=txt_file, alpha=0) 1056 assert allclose(quantity.vertex_values.flat, answer)1058 assert num.allclose(quantity.vertex_values.flat, answer) 1057 1059 1058 1060 #Cleanup … … 1128 1130 answer = linear_function(quantity.domain.get_vertex_coordinates()) 1129 1131 1130 assert allclose(quantity.vertex_values.flat, answer)1132 assert num.allclose(quantity.vertex_values.flat, answer) 1131 1133 1132 1134 1133 1135 #Check that values can be set from file using default attribute 1134 1136 quantity.set_values(filename=ptsfile, alpha=0) 1135 assert allclose(quantity.vertex_values.flat, answer)1137 assert num.allclose(quantity.vertex_values.flat, answer) 1136 1138 1137 1139 #Cleanup … … 1207 1209 1208 1210 1209 assert allclose(quantity.vertex_values.flat, answer)1211 assert num.allclose(quantity.vertex_values.flat, answer) 1210 1212 1211 1213 1212 1214 #Check that values can be set from file using default attribute 1213 1215 quantity.set_values(filename=ptsfile, alpha=0) 1214 assert allclose(quantity.vertex_values.flat, answer)1216 assert num.allclose(quantity.vertex_values.flat, answer) 1215 1217 1216 1218 #Cleanup … … 1226 1228 quantity1.set_vertex_values([0,1,2,3,4,5]) 1227 1229 1228 assert allclose(quantity1.vertex_values,1229 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])1230 assert num.allclose(quantity1.vertex_values, 1231 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1230 1232 1231 1233 1232 1234 quantity2 = Quantity(self.mesh4) 1233 1235 quantity2.set_values(quantity=quantity1) 1234 assert allclose(quantity2.vertex_values,1235 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])1236 assert num.allclose(quantity2.vertex_values, 1237 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1236 1238 1237 1239 quantity2.set_values(quantity = 2*quantity1) 1238 assert allclose(quantity2.vertex_values,1239 [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])1240 assert num.allclose(quantity2.vertex_values, 1241 [[2,0,4], [2,4,8], [8,4,10], [6,2,8]]) 1240 1242 1241 1243 quantity2.set_values(quantity = 2*quantity1 + 3) 1242 assert allclose(quantity2.vertex_values,1243 [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])1244 assert num.allclose(quantity2.vertex_values, 1245 [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 1244 1246 1245 1247 1246 1248 #Check detection of quantity as first orgument 1247 1249 quantity2.set_values(2*quantity1 + 3) 1248 assert allclose(quantity2.vertex_values,1249 [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])1250 assert num.allclose(quantity2.vertex_values, 1251 [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 1250 1252 1251 1253 … … 1262 1264 polygon = [[1.0, 1.0], [4.0, 1.0], 1263 1265 [4.0, 4.0], [1.0, 4.0]] 1264 assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)1266 assert num.allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 1265 1267 1266 1268 quantity1 = Quantity(self.mesh4) 1267 1269 quantity1.set_vertex_values([0,1,2,3,4,5]) 1268 1270 1269 assert allclose(quantity1.vertex_values,1270 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])1271 assert num.allclose(quantity1.vertex_values, 1272 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1271 1273 1272 1274 … … 1276 1278 1277 1279 msg = 'Only node #4(e) at (2,2) should have values applied ' 1278 assert allclose(quantity2.vertex_values,1279 [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg1280 #bac, bce, ecf, dbe1280 assert num.allclose(quantity2.vertex_values, 1281 [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg 1282 #bac, bce, ecf, dbe 1281 1283 1282 1284 … … 1287 1289 quantity1.set_vertex_values([0,1,2,3,4,5]) 1288 1290 1289 assert allclose(quantity1.vertex_values,1290 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])1291 assert num.allclose(quantity1.vertex_values, 1292 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1291 1293 1292 1294 … … 1304 1306 # Negation 1305 1307 Q = -quantity1 1306 assert allclose(Q.vertex_values, -quantity1.vertex_values)1307 assert allclose(Q.centroid_values, -quantity1.centroid_values)1308 assert allclose(Q.edge_values, -quantity1.edge_values)1308 assert num.allclose(Q.vertex_values, -quantity1.vertex_values) 1309 assert num.allclose(Q.centroid_values, -quantity1.centroid_values) 1310 assert num.allclose(Q.edge_values, -quantity1.edge_values) 1309 1311 1310 1312 # Addition 1311 1313 Q = quantity1 + 7 1312 assert allclose(Q.vertex_values, quantity1.vertex_values + 7)1313 assert allclose(Q.centroid_values, quantity1.centroid_values + 7)1314 assert allclose(Q.edge_values, quantity1.edge_values + 7)1314 assert num.allclose(Q.vertex_values, quantity1.vertex_values + 7) 1315 assert num.allclose(Q.centroid_values, quantity1.centroid_values + 7) 1316 assert num.allclose(Q.edge_values, quantity1.edge_values + 7) 1315 1317 1316 1318 Q = 7 + quantity1 1317 assert allclose(Q.vertex_values, quantity1.vertex_values + 7)1318 assert allclose(Q.centroid_values, quantity1.centroid_values + 7)1319 assert allclose(Q.edge_values, quantity1.edge_values + 7)1319 assert num.allclose(Q.vertex_values, quantity1.vertex_values + 7) 1320 assert num.allclose(Q.centroid_values, quantity1.centroid_values + 7) 1321 assert num.allclose(Q.edge_values, quantity1.edge_values + 7) 1320 1322 1321 1323 Q = quantity1 + quantity2 1322 assert allclose(Q.vertex_values,1323 quantity1.vertex_values + quantity2.vertex_values)1324 assert allclose(Q.centroid_values,1325 quantity1.centroid_values + quantity2.centroid_values)1326 assert allclose(Q.edge_values,1327 quantity1.edge_values + quantity2.edge_values)1324 assert num.allclose(Q.vertex_values, 1325 quantity1.vertex_values + quantity2.vertex_values) 1326 assert num.allclose(Q.centroid_values, 1327 quantity1.centroid_values + quantity2.centroid_values) 1328 assert num.allclose(Q.edge_values, 1329 quantity1.edge_values + quantity2.edge_values) 1328 1330 1329 1331 1330 1332 Q = quantity1 + quantity2 - 3 1331 assert allclose(Q.vertex_values,1332 quantity1.vertex_values + quantity2.vertex_values - 3)1333 assert num.allclose(Q.vertex_values, 1334 quantity1.vertex_values + quantity2.vertex_values - 3) 1333 1335 1334 1336 Q = quantity1 - quantity2 1335 assert allclose(Q.vertex_values,1336 quantity1.vertex_values - quantity2.vertex_values)1337 assert num.allclose(Q.vertex_values, 1338 quantity1.vertex_values - quantity2.vertex_values) 1337 1339 1338 1340 #Scaling 1339 1341 Q = quantity1*3 1340 assert allclose(Q.vertex_values, quantity1.vertex_values*3)1341 assert allclose(Q.centroid_values, quantity1.centroid_values*3)1342 assert allclose(Q.edge_values, quantity1.edge_values*3)1342 assert num.allclose(Q.vertex_values, quantity1.vertex_values*3) 1343 assert num.allclose(Q.centroid_values, quantity1.centroid_values*3) 1344 assert num.allclose(Q.edge_values, quantity1.edge_values*3) 1343 1345 Q = 3*quantity1 1344 assert allclose(Q.vertex_values, quantity1.vertex_values*3)1346 assert num.allclose(Q.vertex_values, quantity1.vertex_values*3) 1345 1347 1346 1348 #Multiplication … … 1351 1353 #print quantity2.centroid_values 1352 1354 1353 assert allclose(Q.vertex_values,1354 quantity1.vertex_values * quantity2.vertex_values)1355 assert num.allclose(Q.vertex_values, 1356 quantity1.vertex_values * quantity2.vertex_values) 1355 1357 1356 1358 #Linear combinations 1357 1359 Q = 4*quantity1 + 2 1358 assert allclose(Q.vertex_values,1359 4*quantity1.vertex_values + 2)1360 assert num.allclose(Q.vertex_values, 1361 4*quantity1.vertex_values + 2) 1360 1362 1361 1363 Q = quantity1*quantity2 + 2 1362 assert allclose(Q.vertex_values,1363 quantity1.vertex_values * quantity2.vertex_values + 2)1364 assert num.allclose(Q.vertex_values, 1365 quantity1.vertex_values * quantity2.vertex_values + 2) 1364 1366 1365 1367 Q = quantity1*quantity2 + quantity3 1366 assert allclose(Q.vertex_values,1367 quantity1.vertex_values * quantity2.vertex_values +1368 assert num.allclose(Q.vertex_values, 1369 quantity1.vertex_values * quantity2.vertex_values + 1368 1370 quantity3.vertex_values) 1369 1371 Q = quantity1*quantity2 + 3*quantity3 1370 assert allclose(Q.vertex_values,1371 quantity1.vertex_values * quantity2.vertex_values +1372 3*quantity3.vertex_values)1372 assert num.allclose(Q.vertex_values, 1373 quantity1.vertex_values * quantity2.vertex_values + 1374 3*quantity3.vertex_values) 1373 1375 Q = quantity1*quantity2 + 3*quantity3 + 5.0 1374 assert allclose(Q.vertex_values,1375 quantity1.vertex_values * quantity2.vertex_values +1376 3*quantity3.vertex_values + 5)1376 assert num.allclose(Q.vertex_values, 1377 quantity1.vertex_values * quantity2.vertex_values + 1378 3*quantity3.vertex_values + 5) 1377 1379 1378 1380 Q = quantity1*quantity2 - quantity3 1379 assert allclose(Q.vertex_values,1380 quantity1.vertex_values * quantity2.vertex_values -1381 quantity3.vertex_values)1381 assert num.allclose(Q.vertex_values, 1382 quantity1.vertex_values * quantity2.vertex_values - 1383 quantity3.vertex_values) 1382 1384 Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0 1383 assert allclose(Q.vertex_values,1384 1.5*quantity1.vertex_values * quantity2.vertex_values -1385 3*quantity3.vertex_values + 5)1385 assert num.allclose(Q.vertex_values, 1386 1.5*quantity1.vertex_values * quantity2.vertex_values - 1387 3*quantity3.vertex_values + 5) 1386 1388 1387 1389 #Try combining quantities and arrays and scalars 1388 1390 Q = 1.5*quantity1*quantity2.vertex_values -\ 1389 1391 3*quantity3.vertex_values + 5.0 1390 assert allclose(Q.vertex_values,1391 1.5*quantity1.vertex_values * quantity2.vertex_values -1392 3*quantity3.vertex_values + 5)1392 assert num.allclose(Q.vertex_values, 1393 1.5*quantity1.vertex_values * quantity2.vertex_values - 1394 3*quantity3.vertex_values + 5) 1393 1395 1394 1396 1395 1397 #Powers 1396 1398 Q = quantity1**2 1397 assert allclose(Q.vertex_values, quantity1.vertex_values**2)1399 assert num.allclose(Q.vertex_values, quantity1.vertex_values**2) 1398 1400 1399 1401 Q = quantity1**2 +quantity2**2 1400 assert allclose(Q.vertex_values,1401 quantity1.vertex_values**2 + \1402 quantity2.vertex_values**2)1402 assert num.allclose(Q.vertex_values, 1403 quantity1.vertex_values**2 + \ 1404 quantity2.vertex_values**2) 1403 1405 1404 1406 Q = (quantity1**2 +quantity2**2)**0.5 1405 assert allclose(Q.vertex_values,1406 (quantity1.vertex_values**2 + \1407 quantity2.vertex_values**2)**0.5)1407 assert num.allclose(Q.vertex_values, 1408 (quantity1.vertex_values**2 + \ 1409 quantity2.vertex_values**2)**0.5) 1408 1410 1409 1411 … … 1431 1433 #The central triangle (1) 1432 1434 #(using standard gradient based on neigbours controid values) 1433 assert allclose(a[1], 2.0)1434 assert allclose(b[1], 0.0)1435 assert num.allclose(a[1], 2.0) 1436 assert num.allclose(b[1], 0.0) 1435 1437 1436 1438 … … 1438 1440 #q0 = q1 + a*(x0-x1) + b*(y0-y1) <=> 1439 1441 #2 = 4 + a*(-2/3) + b*(-2/3) 1440 assert allclose(a[0] + b[0], 3)1442 assert num.allclose(a[0] + b[0], 3) 1441 1443 #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 1442 assert allclose(a[0] - b[0], 0)1444 assert num.allclose(a[0] - b[0], 0) 1443 1445 1444 1446 … … 1446 1448 #q2 = q1 + a*(x2-x1) + b*(y2-y1) <=> 1447 1449 #6 = 4 + a*(4/3) + b*(-2/3) 1448 assert allclose(2*a[2] - b[2], 3)1450 assert num.allclose(2*a[2] - b[2], 3) 1449 1451 #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0) 1450 assert allclose(a[2] + 2*b[2], 0)1452 assert num.allclose(a[2] + 2*b[2], 0) 1451 1453 1452 1454 … … 1454 1456 #q3 = q1 + a*(x3-x1) + b*(y3-y1) <=> 1455 1457 #2 = 4 + a*(-2/3) + b*(4/3) 1456 assert allclose(a[3] - 2*b[3], 3)1458 assert num.allclose(a[3] - 2*b[3], 3) 1457 1459 #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 1458 assert allclose(2*a[3] + b[3], 0)1460 assert num.allclose(2*a[3] + b[3], 0) 1459 1461 1460 1462 … … 1464 1466 1465 1467 #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc) 1466 assert allclose(quantity.vertex_values[0,:], [3., 0., 3.])1467 assert allclose(quantity.vertex_values[1,:], [4./3, 16./3, 16./3])1468 assert num.allclose(quantity.vertex_values[0,:], [3., 0., 3.]) 1469 assert num.allclose(quantity.vertex_values[1,:], [4./3, 16./3, 16./3]) 1468 1470 1469 1471 1470 1472 #a = 1.2, b=-0.6 1471 1473 #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3) 1472 assert allclose(quantity.vertex_values[2,2], 8)1474 assert num.allclose(quantity.vertex_values[2,2], 8) 1473 1475 1474 1476 def test_get_gradients(self): … … 1489 1491 #The central triangle (1) 1490 1492 #(using standard gradient based on neigbours controid values) 1491 assert allclose(a[1], 2.0)1492 assert allclose(b[1], 0.0)1493 assert num.allclose(a[1], 2.0) 1494 assert num.allclose(b[1], 0.0) 1493 1495 1494 1496 … … 1496 1498 #q0 = q1 + a*(x0-x1) + b*(y0-y1) <=> 1497 1499 #2 = 4 + a*(-2/3) + b*(-2/3) 1498 assert allclose(a[0] + b[0], 3)1500 assert num.allclose(a[0] + b[0], 3) 1499 1501 #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 1500 assert allclose(a[0] - b[0], 0)1502 assert num.allclose(a[0] - b[0], 0) 1501 1503 1502 1504 … … 1504 1506 #q2 = q1 + a*(x2-x1) + b*(y2-y1) <=> 1505 1507 #6 = 4 + a*(4/3) + b*(-2/3) 1506 assert allclose(2*a[2] - b[2], 3)1508 assert num.allclose(2*a[2] - b[2], 3) 1507 1509 #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0) 1508 assert allclose(a[2] + 2*b[2], 0)1510 assert num.allclose(a[2] + 2*b[2], 0) 1509 1511 1510 1512 … … 1512 1514 #q3 = q1 + a*(x3-x1) + b*(y3-y1) <=> 1513 1515 #2 = 4 + a*(-2/3) + b*(4/3) 1514 assert allclose(a[3] - 2*b[3], 3)1516 assert num.allclose(a[3] - 2*b[3], 3) 1515 1517 #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 1516 assert allclose(2*a[3] + b[3], 0)1518 assert num.allclose(2*a[3] + b[3], 0) 1517 1519 1518 1520 … … 1532 1534 #print a, b 1533 1535 1534 assert allclose(a[1], 3.0)1535 assert allclose(b[1], 1.0)1536 assert num.allclose(a[1], 3.0) 1537 assert num.allclose(b[1], 1.0) 1536 1538 1537 1539 #Work out the others … … 1540 1542 1541 1543 #print quantity.vertex_values 1542 assert allclose(quantity.vertex_values[1,0], 2.0)1543 assert allclose(quantity.vertex_values[1,1], 6.0)1544 assert allclose(quantity.vertex_values[1,2], 8.0)1544 assert num.allclose(quantity.vertex_values[1,0], 2.0) 1545 assert num.allclose(quantity.vertex_values[1,1], 6.0) 1546 assert num.allclose(quantity.vertex_values[1,2], 8.0) 1545 1547 1546 1548 … … 1550 1552 1551 1553 #Set up for a gradient of (3,1), f(x) = 3x+y 1552 c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])1553 d_values = array([1.0, 2.0, 3.0, 4.0])1554 c_values = num.array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3]) 1555 d_values = num.array([1.0, 2.0, 3.0, 4.0]) 1554 1556 quantity.set_values(c_values, location = 'centroids') 1555 1557 … … 1558 1560 1559 1561 #print quantity.vertex_values 1560 assert allclose(quantity.centroid_values, quantity.centroid_backup_values)1562 assert num.allclose(quantity.centroid_values, quantity.centroid_backup_values) 1561 1563 1562 1564 … … 1574 1576 #Test centroids 1575 1577 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1576 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid1578 assert num.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1577 1579 1578 1580 #Extrapolate … … 1581 1583 #Check that gradient is zero 1582 1584 a,b = quantity.get_gradients() 1583 assert allclose(a, [0,0,0,0])1584 assert allclose(b, [0,0,0,0])1585 assert num.allclose(a, [0,0,0,0]) 1586 assert num.allclose(b, [0,0,0,0]) 1585 1587 1586 1588 #Check vertices but not edge values 1587 assert allclose(quantity.vertex_values,1588 [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])1589 assert num.allclose(quantity.vertex_values, 1590 [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 1589 1591 1590 1592 … … 1614 1616 1615 1617 #Assert that quantities are conserved 1616 from Numeric import sum1617 1618 for k in range(quantity.centroid_values.shape[0]): 1618 assert allclose (quantity.centroid_values[k],1619 sum(quantity.vertex_values[k,:])/3)1619 assert num.allclose (quantity.centroid_values[k], 1620 num.sum(quantity.vertex_values[k,:])/3) 1620 1621 1621 1622 … … 1646 1647 1647 1648 #Assert that quantities are conserved 1648 from Numeric import sum1649 1649 for k in range(quantity.centroid_values.shape[0]): 1650 assert allclose (quantity.centroid_values[k],1651 sum(quantity.vertex_values[k,:])/3)1650 assert num.allclose (quantity.centroid_values[k], 1651 num.sum(quantity.vertex_values[k,:])/3) 1652 1652 1653 1653 … … 1676 1676 1677 1677 #Assert that quantities are conserved 1678 from Numeric import sum1679 1678 for k in range(quantity.centroid_values.shape[0]): 1680 assert allclose (quantity.centroid_values[k],1681 sum(quantity.vertex_values[k,:])/3)1679 assert num.allclose (quantity.centroid_values[k], 1680 num.sum(quantity.vertex_values[k,:])/3) 1682 1681 1683 1682 … … 1705 1704 1706 1705 #Assert that quantities are conserved 1707 from Numeric import sum1708 1706 for k in range(quantity.centroid_values.shape[0]): 1709 assert allclose (quantity.centroid_values[k],1710 sum(quantity.vertex_values[k,:])/3)1707 assert num.allclose (quantity.centroid_values[k], 1708 num.sum(quantity.vertex_values[k,:])/3) 1711 1709 1712 1710 def test_limiter2(self): … … 1718 1716 #Test centroids 1719 1717 quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1720 assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid1718 assert num.allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1721 1719 1722 1720 … … 1724 1722 quantity.extrapolate_second_order() 1725 1723 1726 assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])1724 assert num.allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1727 1725 1728 1726 #Limit … … 1731 1729 # limited value for beta_w = 0.9 1732 1730 1733 assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])1731 assert num.allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9]) 1734 1732 # limited values for beta_w = 0.5 1735 1733 #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5]) … … 1737 1735 1738 1736 #Assert that quantities are conserved 1739 from Numeric import sum1740 1737 for k in range(quantity.centroid_values.shape[0]): 1741 assert allclose (quantity.centroid_values[k],1742 sum(quantity.vertex_values[k,:])/3)1738 assert num.allclose (quantity.centroid_values[k], 1739 num.sum(quantity.vertex_values[k,:])/3) 1743 1740 1744 1741 … … 1751 1748 #Test centroids 1752 1749 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1753 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid1750 assert num.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1754 1751 1755 1752 … … 1760 1757 #quantity.interpolate_from_vertices_to_edges() 1761 1758 1762 assert allclose(quantity.vertex_values,1763 [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])1764 assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],1765 [3,3,3], [4, 4, 4]])1759 assert num.allclose(quantity.vertex_values, 1760 [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 1761 assert num.allclose(quantity.edge_values, [[1,1,1], [2,2,2], 1762 [3,3,3], [4, 4, 4]]) 1766 1763 1767 1764 … … 1769 1766 quantity = Quantity(self.mesh4) 1770 1767 1771 quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)1768 quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],num.Float) 1772 1769 1773 1770 quantity.interpolate_from_vertices_to_edges() 1774 1771 1775 assert allclose(quantity.edge_values, [[1., 1.5, 0.5],1776 [3., 2.5, 1.5],1777 [3.5, 4.5, 3.],1778 [2.5, 3.5, 2]])1772 assert num.allclose(quantity.edge_values, [[1., 1.5, 0.5], 1773 [3., 2.5, 1.5], 1774 [3.5, 4.5, 3.], 1775 [2.5, 3.5, 2]]) 1779 1776 1780 1777 … … 1782 1779 quantity = Quantity(self.mesh4) 1783 1780 1784 quantity.edge_values = array([[1., 1.5, 0.5],1785 [3., 2.5, 1.5],1786 [3.5, 4.5, 3.],1787 [2.5, 3.5, 2]],Float)1781 quantity.edge_values = num.array([[1., 1.5, 0.5], 1782 [3., 2.5, 1.5], 1783 [3.5, 4.5, 3.], 1784 [2.5, 3.5, 2]],num.Float) 1788 1785 1789 1786 quantity.interpolate_from_edges_to_vertices() 1790 1787 1791 assert allclose(quantity.vertex_values,1792 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])1788 assert num.allclose(quantity.vertex_values, 1789 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1793 1790 1794 1791 … … 1799 1796 #Test centroids 1800 1797 quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1801 assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid1798 assert num.allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1802 1799 1803 1800 … … 1805 1802 quantity.extrapolate_second_order() 1806 1803 1807 assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])1804 assert num.allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1808 1805 1809 1806 … … 1813 1810 #Test centroids 1814 1811 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1815 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid1812 assert num.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1816 1813 1817 1814 #Set explicit_update 1818 quantity.explicit_update = array( [1.,1.,1.,1.] )1815 quantity.explicit_update = num.array( [1.,1.,1.,1.] ) 1819 1816 1820 1817 #Update with given timestep 1821 1818 quantity.update(0.1) 1822 1819 1823 x = array([1, 2, 3, 4]) +array( [.1,.1,.1,.1] )1824 assert allclose( quantity.centroid_values, x)1820 x = num.array([1, 2, 3, 4]) + num.array( [.1,.1,.1,.1] ) 1821 assert num.allclose( quantity.centroid_values, x) 1825 1822 1826 1823 def test_update_semi_implicit(self): … … 1829 1826 #Test centroids 1830 1827 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1831 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid1828 assert num.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1832 1829 1833 1830 #Set semi implicit update 1834 quantity.semi_implicit_update = array([1.,1.,1.,1.])1831 quantity.semi_implicit_update = num.array([1.,1.,1.,1.]) 1835 1832 1836 1833 #Update with given timestep … … 1838 1835 quantity.update(timestep) 1839 1836 1840 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])1841 denom = ones(4,Float)-timestep*sem1842 1843 x = array([1, 2, 3, 4])/denom1844 assert allclose( quantity.centroid_values, x)1837 sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4]) 1838 denom = num.ones(4, num.Float)-timestep*sem 1839 1840 x = num.array([1, 2, 3, 4])/denom 1841 assert num.allclose( quantity.centroid_values, x) 1845 1842 1846 1843 … … 1850 1847 #Test centroids 1851 1848 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1852 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid1849 assert num.allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1853 1850 1854 1851 #Set explicit_update 1855 quantity.explicit_update = array( [4.,3.,2.,1.] )1852 quantity.explicit_update = num.array( [4.,3.,2.,1.] ) 1856 1853 1857 1854 #Set semi implicit update 1858 quantity.semi_implicit_update = array( [1.,1.,1.,1.] )1855 quantity.semi_implicit_update = num.array( [1.,1.,1.,1.] ) 1859 1856 1860 1857 #Update with given timestep … … 1862 1859 quantity.update(0.1) 1863 1860 1864 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])1865 denom = ones(4,Float)-timestep*sem1866 1867 x = array([1., 2., 3., 4.])1861 sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4]) 1862 denom = num.ones(4, num.Float)-timestep*sem 1863 1864 x = num.array([1., 2., 3., 4.]) 1868 1865 x /= denom 1869 x += timestep* array( [4.0, 3.0, 2.0, 1.0] )1870 1871 assert allclose( quantity.centroid_values, x)1866 x += timestep*num.array( [4.0, 3.0, 2.0, 1.0] ) 1867 1868 assert num.allclose( quantity.centroid_values, x) 1872 1869 1873 1870 … … 1879 1876 from mesh_factory import rectangular 1880 1877 from shallow_water import Domain, Transmissive_boundary 1881 from Numeric import zeros, Float1882 1878 from anuga.utilities.numerical_tools import mean 1883 1879 … … 1906 1902 1907 1903 bed = domain.quantities['elevation'].vertex_values 1908 stage = zeros(bed.shape,Float)1904 stage = num.zeros(bed.shape, num.Float) 1909 1905 1910 1906 h = 0.03 … … 1929 1925 1930 1926 #First four points 1931 assert allclose(A[0], (Q[0,2] + Q[1,1])/2)1932 assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)1933 assert allclose(A[2], Q[3,0])1934 assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)1927 assert num.allclose(A[0], (Q[0,2] + Q[1,1])/2) 1928 assert num.allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3) 1929 assert num.allclose(A[2], Q[3,0]) 1930 assert num.allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3) 1935 1931 1936 1932 #Center point 1937 assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\1938 Q[5,0] + Q[6,2] + Q[7,1])/6)1933 assert num.allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\ 1934 Q[5,0] + Q[6,2] + Q[7,1])/6) 1939 1935 1940 1936 1941 1937 #Check V 1942 assert allclose(V[0,:], [3,4,0])1943 assert allclose(V[1,:], [1,0,4])1944 assert allclose(V[2,:], [4,5,1])1945 assert allclose(V[3,:], [2,1,5])1946 assert allclose(V[4,:], [6,7,3])1947 assert allclose(V[5,:], [4,3,7])1948 assert allclose(V[6,:], [7,8,4])1949 assert allclose(V[7,:], [5,4,8])1938 assert num.allclose(V[0,:], [3,4,0]) 1939 assert num.allclose(V[1,:], [1,0,4]) 1940 assert num.allclose(V[2,:], [4,5,1]) 1941 assert num.allclose(V[3,:], [2,1,5]) 1942 assert num.allclose(V[4,:], [6,7,3]) 1943 assert num.allclose(V[5,:], [4,3,7]) 1944 assert num.allclose(V[6,:], [7,8,4]) 1945 assert num.allclose(V[7,:], [5,4,8]) 1950 1946 1951 1947 #Get smoothed stage with XY 1952 1948 X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True) 1953 1949 1954 assert allclose(A, A1)1955 assert allclose(V, V1)1950 assert num.allclose(A, A1) 1951 assert num.allclose(V, V1) 1956 1952 1957 1953 #Check XY 1958 assert allclose(X[4], 0.5)1959 assert allclose(Y[4], 0.5)1960 1961 assert allclose(X[7], 1.0)1962 assert allclose(Y[7], 0.5)1954 assert num.allclose(X[4], 0.5) 1955 assert num.allclose(Y[4], 0.5) 1956 1957 assert num.allclose(X[7], 1.0) 1958 assert num.allclose(Y[7], 0.5) 1963 1959 1964 1960 … … 1969 1965 from mesh_factory import rectangular 1970 1966 from shallow_water import Domain, Transmissive_boundary 1971 from Numeric import zeros, Float1972 1967 from anuga.utilities.numerical_tools import mean 1973 1968 … … 1991 1986 1992 1987 bed = domain.quantities['elevation'].vertex_values 1993 stage = zeros(bed.shape,Float)1988 stage = num.zeros(bed.shape, num.Float) 1994 1989 1995 1990 h = 0.03 … … 2008 2003 2009 2004 for k in range(8): 2010 assert allclose(A[k], Q[k])2005 assert num.allclose(A[k], Q[k]) 2011 2006 2012 2007 … … 2021 2016 2022 2017 2023 assert allclose(A, A1)2024 assert allclose(V, V1)2018 assert num.allclose(A, A1) 2019 assert num.allclose(V, V1) 2025 2020 2026 2021 #Check XY 2027 assert allclose(X[1], 0.5)2028 assert allclose(Y[1], 0.5)2029 assert allclose(X[4], 0.0)2030 assert allclose(Y[4], 0.0)2031 assert allclose(X[12], 1.0)2032 assert allclose(Y[12], 0.0)2022 assert num.allclose(X[1], 0.5) 2023 assert num.allclose(Y[1], 0.5) 2024 assert num.allclose(X[4], 0.0) 2025 assert num.allclose(Y[4], 0.0) 2026 assert num.allclose(X[12], 1.0) 2027 assert num.allclose(Y[12], 0.0) 2033 2028 2034 2029 … … 2038 2033 from mesh_factory import rectangular 2039 2034 from shallow_water import Domain 2040 from Numeric import zeros, Float2041 2035 2042 2036 #Create basic mesh … … 2054 2048 #print "quantity.centroid_values",quantity.centroid_values 2055 2049 2056 assert allclose(quantity.centroid_values, [1,7])2050 assert num.allclose(quantity.centroid_values, [1,7]) 2057 2051 2058 2052 quantity.set_array_values([15,20,25], indices = indices) 2059 assert allclose(quantity.centroid_values, [1,20])2053 assert num.allclose(quantity.centroid_values, [1,20]) 2060 2054 2061 2055 quantity.set_array_values([15,20,25], indices = indices) 2062 assert allclose(quantity.centroid_values, [1,20])2056 assert num.allclose(quantity.centroid_values, [1,20]) 2063 2057 2064 2058 def test_setting_some_vertex_values(self): … … 2068 2062 from mesh_factory import rectangular 2069 2063 from shallow_water import Domain 2070 from Numeric import zeros, Float2071 2064 2072 2065 #Create basic mesh … … 2087 2080 indices = indices) 2088 2081 #print "quantity.centroid_values",quantity.centroid_values 2089 assert allclose(quantity.centroid_values, [1,7,3,4,5,6])2082 assert num.allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2090 2083 2091 2084 value = [7] … … 2095 2088 indices = indices) 2096 2089 #print "quantity.centroid_values",quantity.centroid_values 2097 assert allclose(quantity.centroid_values, [1,7,3,4,5,6])2090 assert num.allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2098 2091 2099 2092 value = [[15,20,25]] 2100 2093 quantity.set_values(value, indices = indices) 2101 2094 #print "1 quantity.vertex_values",quantity.vertex_values 2102 assert allclose(quantity.vertex_values[1], value[0])2095 assert num.allclose(quantity.vertex_values[1], value[0]) 2103 2096 2104 2097 … … 2107 2100 quantity.set_values(values, indices = [0,1,5], location = 'centroids') 2108 2101 #print "2 quantity.vertex_values",quantity.vertex_values 2109 assert allclose(quantity.vertex_values[0], [10,10,10])2110 assert allclose(quantity.vertex_values[5], [50,50,50])2102 assert num.allclose(quantity.vertex_values[0], [10,10,10]) 2103 assert num.allclose(quantity.vertex_values[5], [50,50,50]) 2111 2104 #quantity.interpolate() 2112 2105 #print "quantity.centroid_values",quantity.centroid_values 2113 assert allclose(quantity.centroid_values, [10,100,3,4,5,50])2106 assert num.allclose(quantity.centroid_values, [10,100,3,4,5,50]) 2114 2107 2115 2108 … … 2121 2114 quantity.set_values(values, indices = [0,1,5]) 2122 2115 #print "quantity.vertex_values",quantity.vertex_values 2123 assert allclose(quantity.vertex_values[0], [1,50,10])2124 assert allclose(quantity.vertex_values[5], [6,6,6])2125 assert allclose(quantity.vertex_values[1], [100,10,50])2116 assert num.allclose(quantity.vertex_values[0], [1,50,10]) 2117 assert num.allclose(quantity.vertex_values[5], [6,6,6]) 2118 assert num.allclose(quantity.vertex_values[1], [100,10,50]) 2126 2119 2127 2120 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], … … 2130 2123 quantity.set_values(values, indices = [3,3,5]) 2131 2124 quantity.interpolate() 2132 assert allclose(quantity.centroid_values, [1,2,3,400,5,999])2125 assert num.allclose(quantity.centroid_values, [1,2,3,400,5,999]) 2133 2126 2134 2127 values = [[1,1,1],[2,2,2],[3,3,3], … … 2142 2135 quantity.set_values(values) 2143 2136 #print "1 quantity.vertex_values",quantity.vertex_values 2144 assert allclose(quantity.vertex_values,[[ 4., 5., 0.],2145 [ 1., 0., 5.],2146 [ 5., 6., 1.],2147 [ 2., 1., 6.],2148 [ 6., 7., 2.],2149 [ 3., 2., 7.]])2137 assert num.allclose(quantity.vertex_values,[[ 4., 5., 0.], 2138 [ 1., 0., 5.], 2139 [ 5., 6., 1.], 2140 [ 2., 1., 6.], 2141 [ 6., 7., 2.], 2142 [ 3., 2., 7.]]) 2150 2143 2151 2144 def test_setting_unique_vertex_values(self): … … 2155 2148 from mesh_factory import rectangular 2156 2149 from shallow_water import Domain 2157 from Numeric import zeros, Float2158 2150 2159 2151 #Create basic mesh … … 2171 2163 indices = indices) 2172 2164 #print "quantity.centroid_values",quantity.centroid_values 2173 assert allclose(quantity.vertex_values[0], [0,7,0])2174 assert allclose(quantity.vertex_values[1], [7,1,7])2175 assert allclose(quantity.vertex_values[2], [7,2,7])2165 assert num.allclose(quantity.vertex_values[0], [0,7,0]) 2166 assert num.allclose(quantity.vertex_values[1], [7,1,7]) 2167 assert num.allclose(quantity.vertex_values[2], [7,2,7]) 2176 2168 2177 2169 … … 2182 2174 from mesh_factory import rectangular 2183 2175 from shallow_water import Domain 2184 from Numeric import zeros, Float2185 2176 2186 2177 #Create basic mesh … … 2205 2196 2206 2197 answer = [0.5,2,4,5,0,1,3,4.5] 2207 assert allclose(answer,2208 quantity.get_values(location = 'unique vertices'))2198 assert num.allclose(answer, 2199 quantity.get_values(location = 'unique vertices')) 2209 2200 2210 2201 indices = [0,5,3] 2211 2202 answer = [0.5,1,5] 2212 assert allclose(answer,2213 quantity.get_values(indices=indices, \2214 location = 'unique vertices'))2203 assert num.allclose(answer, 2204 quantity.get_values(indices=indices, 2205 location = 'unique vertices')) 2215 2206 #print "quantity.centroid_values",quantity.centroid_values 2216 2207 #print "quantity.get_values(location = 'centroids') ",\ … … 2241 2232 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2242 2233 2243 assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])2244 assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])2245 2246 2247 assert allclose(quantity.get_values(location='vertices'), [[4,0,2],2248 [4,2,6],2249 [6,2,4],2250 [8,4,6]])2251 2252 assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],2253 [8,4,6]])2254 2255 2256 assert allclose(quantity.get_values(location='edges'), [[1,3,2],2257 [4,5,3],2258 [3,5,4],2259 [5,7,6]])2260 assert allclose(quantity.get_values(location='edges', indices=[1,3]),2261 [[4,5,3],2262 [5,7,6]])2234 assert num.allclose(quantity.get_values(location='centroids'), [2,4,4,6]) 2235 assert num.allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6]) 2236 2237 2238 assert num.allclose(quantity.get_values(location='vertices'), [[4,0,2], 2239 [4,2,6], 2240 [6,2,4], 2241 [8,4,6]]) 2242 2243 assert num.allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6], 2244 [8,4,6]]) 2245 2246 2247 assert num.allclose(quantity.get_values(location='edges'), [[1,3,2], 2248 [4,5,3], 2249 [3,5,4], 2250 [5,7,6]]) 2251 assert num.allclose(quantity.get_values(location='edges', indices=[1,3]), 2252 [[4,5,3], 2253 [5,7,6]]) 2263 2254 2264 2255 # Check averaging over vertices … … 2280 2271 from mesh_factory import rectangular 2281 2272 from shallow_water import Domain 2282 from Numeric import zeros, Float2283 2273 2284 2274 #Create basic mesh … … 2298 2288 2299 2289 #print quantity.get_values(points=interpolation_points) 2300 assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))2290 assert num.allclose(answer, quantity.get_values(interpolation_points=interpolation_points)) 2301 2291 2302 2292 … … 2311 2301 #print answer 2312 2302 #print quantity.get_values(interpolation_points=interpolation_points) 2313 assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points,2314 verbose=False))2303 assert num.allclose(answer, quantity.get_values(interpolation_points=interpolation_points, 2304 verbose=False)) 2315 2305 2316 2306 … … 2345 2335 x, y = 2.0/3, 8.0/3 2346 2336 v = quantity.get_values(interpolation_points = [[x,y]]) 2347 assert allclose(v, 6)2337 assert num.allclose(v, 6) 2348 2338 2349 2339 # Then another to test that algorithm won't blindly … … 2351 2341 x, y = 4.0/3, 4.0/3 2352 2342 v = quantity.get_values(interpolation_points = [[x,y]]) 2353 assert allclose(v, 4)2343 assert num.allclose(v, 4) 2354 2344 2355 2345 … … 2380 2370 x, y = 2.0/3, 8.0/3 2381 2371 v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2382 assert allclose(v, 6)2372 assert num.allclose(v, 6) 2383 2373 2384 2374 … … 2387 2377 x, y = 4.0/3, 4.0/3 2388 2378 v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2389 assert allclose(v, 4)2379 assert num.allclose(v, 4) 2390 2380 2391 2381 # Try two points … … 2393 2383 [4.0/3 + xllcorner, 4.0/3 + yllcorner]] 2394 2384 v = quantity.get_values(interpolation_points=pts) 2395 assert allclose(v, [6, 4])2385 assert num.allclose(v, [6, 4]) 2396 2386 2397 2387 # Test it using the geospatial data format with absolute input points and default georef 2398 2388 pts = Geospatial_data(data_points=pts) 2399 2389 v = quantity.get_values(interpolation_points=pts) 2400 assert allclose(v, [6, 4])2390 assert num.allclose(v, [6, 4]) 2401 2391 2402 2392 … … 2405 2395 geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2406 2396 v = quantity.get_values(interpolation_points=pts) 2407 assert allclose(v, [6, 4])2397 assert num.allclose(v, [6, 4]) 2408 2398 2409 2399 … … 2416 2406 from mesh_factory import rectangular 2417 2407 from shallow_water import Domain 2418 from Numeric import zeros, Float2419 2408 2420 2409 #Create basic mesh … … 2438 2427 #print "quantity.get_values(location = 'centroids') ",\ 2439 2428 # quantity.get_values(location = 'centroids') 2440 assert allclose(quantity.centroid_values,2441 quantity.get_values(location = 'centroids'))2429 assert num.allclose(quantity.centroid_values, 2430 quantity.get_values(location = 'centroids')) 2442 2431 2443 2432 … … 2445 2434 quantity.set_values(value, indices = indices) 2446 2435 #print "1 quantity.vertex_values",quantity.vertex_values 2447 assert allclose(quantity.vertex_values, quantity.get_values())2448 2449 assert allclose(quantity.edge_values,2450 quantity.get_values(location = 'edges'))2436 assert num.allclose(quantity.vertex_values, quantity.get_values()) 2437 2438 assert num.allclose(quantity.edge_values, 2439 quantity.get_values(location = 'edges')) 2451 2440 2452 2441 # get a subset of elements 2453 2442 subset = quantity.get_values(location='centroids', indices=[0,5]) 2454 2443 answer = [quantity.centroid_values[0],quantity.centroid_values[5]] 2455 assert allclose(subset, answer)2444 assert num.allclose(subset, answer) 2456 2445 2457 2446 … … 2460 2449 #print "subset",subset 2461 2450 #print "answer",answer 2462 assert allclose(subset, answer)2451 assert num.allclose(subset, answer) 2463 2452 2464 2453 subset = quantity.get_values( indices=[1,5]) … … 2466 2455 #print "subset",subset 2467 2456 #print "answer",answer 2468 assert allclose(subset, answer)2457 assert num.allclose(subset, answer) 2469 2458 2470 2459 def test_smooth_vertex_values(self): … … 2474 2463 from mesh_factory import rectangular 2475 2464 from shallow_water import Domain 2476 from Numeric import zeros, Float2477 2465 2478 2466 #Create basic mesh … … 2512 2500 [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]] 2513 2501 2514 assert allclose(answer_vertex_values,2515 quantity.vertex_values)2502 assert num.allclose(answer_vertex_values, 2503 quantity.vertex_values) 2516 2504 #print "quantity.centroid_values",quantity.centroid_values 2517 2505 #print "quantity.get_values(location = 'centroids') ",\ -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py
r5897 r6145 7 7 from region import * 8 8 #from anuga.config import epsilon 9 from Numeric import allclose, average #, array, ones, Float 9 10 import Numeric as num 11 12 10 13 """ 11 14 This is what the mesh in these tests look like; … … 44 47 from mesh_factory import rectangular 45 48 from shallow_water import Domain 46 from Numeric import zeros, Float47 49 48 50 #Create basic mesh … … 64 66 domain.set_region([a, b]) 65 67 #print domain.quantities['friction'].get_values() 66 assert allclose(domain.quantities['friction'].get_values(),\67 [[ 0.09, 0.09, 0.09],68 [ 0.09, 0.09, 0.09],69 [ 0.07, 0.07, 0.07],70 [ 0.07, 0.07, 0.07],71 [ 1.0, 1.0, 1.0],72 [ 1.0, 1.0, 1.0]])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]]) 73 75 74 76 #c = Add_Value_To_region('all', 'friction', 10.0) 75 77 domain.set_region(Add_value_to_region('all', 'friction', 10.0)) 76 78 #print domain.quantities['friction'].get_values() 77 assert allclose(domain.quantities['friction'].get_values(),78 [[ 10.09, 10.09, 10.09],79 [ 10.09, 10.09, 10.09],80 [ 10.07, 10.07, 10.07],81 [ 10.07, 10.07, 10.07],82 [ 11.0, 11.0, 11.0],83 [ 11.0, 11.0, 11.0]])79 assert num.allclose(domain.quantities['friction'].get_values(), 80 [[ 10.09, 10.09, 10.09], 81 [ 10.09, 10.09, 10.09], 82 [ 10.07, 10.07, 10.07], 83 [ 10.07, 10.07, 10.07], 84 [ 11.0, 11.0, 11.0], 85 [ 11.0, 11.0, 11.0]]) 84 86 85 87 # trying a function 86 88 domain.set_region(Set_region('top', 'friction', add_x_y)) 87 89 #print domain.quantities['friction'].get_values() 88 assert allclose(domain.quantities['friction'].get_values(),89 [[ 10.09, 10.09, 10.09],90 [ 10.09, 10.09, 10.09],91 [ 10.07, 10.07, 10.07],92 [ 10.07, 10.07, 10.07],93 [ 5./3, 2.0, 2./3],94 [ 1.0, 2./3, 2.0]])90 assert num.allclose(domain.quantities['friction'].get_values(), 91 [[ 10.09, 10.09, 10.09], 92 [ 10.09, 10.09, 10.09], 93 [ 10.07, 10.07, 10.07], 94 [ 10.07, 10.07, 10.07], 95 [ 5./3, 2.0, 2./3], 96 [ 1.0, 2./3, 2.0]]) 95 97 96 98 domain.set_quantity('elevation', 10.0) … … 98 100 domain.set_region(Add_value_to_region('top', 'stage', 1.0,initial_quantity='elevation')) 99 101 #print domain.quantities['stage'].get_values() 100 assert allclose(domain.quantities['stage'].get_values(),101 [[ 10., 10., 10.],102 [ 10., 10., 10.],103 [ 10., 10., 10.],104 [ 10., 10., 10.],105 [ 11.0, 11.0, 11.0],106 [ 11.0, 11.0, 11.0]])102 assert num.allclose(domain.quantities['stage'].get_values(), 103 [[ 10., 10., 10.], 104 [ 10., 10., 10.], 105 [ 10., 10., 10.], 106 [ 10., 10., 10.], 107 [ 11.0, 11.0, 11.0], 108 [ 11.0, 11.0, 11.0]]) 107 109 108 110 … … 115 117 domain.set_region(Add_quantities('top', 'elevation','stage')) 116 118 #print domain.quantities['stage'].get_values() 117 assert allclose(domain.quantities['elevation'].get_values(),118 [[ 10., 10., 10.],119 [ 10., 10., 10.],120 [ 10., 10., 10.],121 [ 10., 10., 10.],122 [ 33., 33.0, 33.],123 [ 33.0, 33., 33.]])119 assert num.allclose(domain.quantities['elevation'].get_values(), 120 [[ 10., 10., 10.], 121 [ 10., 10., 10.], 122 [ 10., 10., 10.], 123 [ 10., 10., 10.], 124 [ 33., 33.0, 33.], 125 [ 33.0, 33., 33.]]) 124 126 125 127 def test_unique_vertices(self): … … 129 131 from mesh_factory import rectangular 130 132 from shallow_water import Domain 131 from Numeric import zeros, Float132 133 133 134 #Create basic mesh … … 147 148 domain.set_region(a) 148 149 #print domain.quantities['friction'].get_values() 149 assert allclose(domain.quantities['friction'].get_values(),\150 [[ 0.09, 0.09, 0.09],151 [ 0.09, 0.09, 0.09],152 [ 0.09, 0.07, 0.09],153 [ 0.07, 0.09, 0.07],154 [ 0.07, 0.07, 0.07],155 [ 0.07, 0.07, 0.07]])150 assert num.allclose(domain.quantities['friction'].get_values(),\ 151 [[ 0.09, 0.09, 0.09], 152 [ 0.09, 0.09, 0.09], 153 [ 0.09, 0.07, 0.09], 154 [ 0.07, 0.09, 0.07], 155 [ 0.07, 0.07, 0.07], 156 [ 0.07, 0.07, 0.07]]) 156 157 157 158 … … 162 163 from mesh_factory import rectangular 163 164 from shallow_water import Domain 164 from Numeric import zeros, Float165 165 166 166 #Create basic mesh … … 180 180 181 181 #print domain.quantities['friction'].get_values() 182 assert allclose(domain.quantities['friction'].get_values(),\183 [[ 1.07, 1.07, 1.07],184 [ 1.07, 1.07, 1.07],185 [ 1.07, 0.07, 1.07],186 [ 0.07, 1.07, 0.07],187 [ 0.07, 0.07, 0.07],188 [ 0.07, 0.07, 0.07]])182 assert num.allclose(domain.quantities['friction'].get_values(),\ 183 [[ 1.07, 1.07, 1.07], 184 [ 1.07, 1.07, 1.07], 185 [ 1.07, 0.07, 1.07], 186 [ 0.07, 1.07, 0.07], 187 [ 0.07, 0.07, 0.07], 188 [ 0.07, 0.07, 0.07]]) 189 189 190 190 def test_unique_vertices_average_loc_vert(self): … … 194 194 from mesh_factory import rectangular 195 195 from shallow_water import Domain 196 from Numeric import zeros, Float197 196 198 197 #Create basic mesh … … 218 217 #print domain.quantities['friction'].get_values() 219 218 frict_points = domain.quantities['friction'].get_values() 220 assert allclose(frict_points[0],\221 [ calc_frict, calc_frict, calc_frict])222 assert allclose(frict_points[1],\223 [ calc_frict, calc_frict, calc_frict])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]) 224 223 225 224 def test_unique_vertices_average_loc_unique_vert(self): … … 229 228 from mesh_factory import rectangular 230 229 from shallow_water import Domain 231 from Numeric import zeros, Float232 230 233 231 #Create basic mesh … … 252 250 #print domain.quantities['friction'].get_values() 253 251 frict_points = domain.quantities['friction'].get_values() 254 assert allclose(frict_points[0],\255 [ calc_frict, calc_frict, calc_frict])256 assert allclose(frict_points[1],\257 [ calc_frict, calc_frict, calc_frict])258 assert allclose(frict_points[2],\259 [ calc_frict, 1.0 + 2.0/3.0, calc_frict])260 assert allclose(frict_points[3],\261 [ 2.0/3.0,calc_frict, 1.0 + 2.0/3.0])252 assert num.allclose(frict_points[0],\ 253 [ calc_frict, calc_frict, calc_frict]) 254 assert num.allclose(frict_points[1],\ 255 [ calc_frict, calc_frict, calc_frict]) 256 assert num.allclose(frict_points[2],\ 257 [ calc_frict, 1.0 + 2.0/3.0, calc_frict]) 258 assert num.allclose(frict_points[3],\ 259 [ 2.0/3.0,calc_frict, 1.0 + 2.0/3.0]) 262 260 263 261 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py
r6086 r6145 3 3 4 4 import unittest 5 from Numeric import zeros, array, allclose, Float6 5 from math import sqrt, pi 7 6 import tempfile, os … … 22 21 import time 23 22 import string 23 24 import Numeric as num 25 24 26 25 27 def test_function(x, y): … … 91 93 92 94 #Exact linear intpolation 93 assert allclose(q[0], 2*t)95 assert num.allclose(q[0], 2*t) 94 96 if i%6 == 0: 95 assert allclose(q[1], t**2)96 assert allclose(q[2], sin(t*pi/600))97 assert num.allclose(q[1], t**2) 98 assert num.allclose(q[2], sin(t*pi/600)) 97 99 98 100 #Check non-exact … … 100 102 t = 90 #Halfway between 60 and 120 101 103 q = F(t) 102 assert allclose( (120**2 + 60**2)/2, q[1] )103 assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )104 assert num.allclose( (120**2 + 60**2)/2, q[1] ) 105 assert num.allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] ) 104 106 105 107 106 108 t = 100 #Two thirds of the way between between 60 and 120 107 109 q = F(t) 108 assert allclose( 2*120**2/3 + 60**2/3, q[1] )109 assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )110 assert num.allclose( 2*120**2/3 + 60**2/3, q[1] ) 111 assert num.allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] ) 110 112 111 113 os.remove(filename + '.txt') … … 125 127 from shallow_water import Domain, Dirichlet_boundary 126 128 from mesh_factory import rectangular 127 from Numeric import take, concatenate, reshape128 129 129 130 #Create basic mesh and shallow water domain … … 182 183 183 184 last_time_index = len(time)-1 #Last last_time_index 184 d_stage = reshape(take(stage[last_time_index, :], [0,5,10,15]), (4,1))185 d_uh = reshape(take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))186 d_vh = reshape(take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))187 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 188 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 188 189 189 190 #Reference interpolated values at midpoints on diagonal at … … 194 195 195 196 #And the midpoints are found now 196 Dx = take(reshape(x, (16,1)), [0,5,10,15])197 Dy = take(reshape(y, (16,1)), [0,5,10,15])198 199 diag = concatenate( (Dx, Dy), axis=1)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]) 199 200 diag = num.concatenate( (Dx, Dy), axis=1) 200 201 d_midpoints = (diag[1:] + diag[:-1])/2 201 202 … … 209 210 assert not T[-1] == T[-2], msg 210 211 t = time[last_time_index] 211 q = f(t, point_id=0); assert allclose(r0, q)212 q = f(t, point_id=1); assert allclose(r1, q)213 q = f(t, point_id=2); assert allclose(r2, q)212 q = f(t, point_id=0); assert num.allclose(r0, q) 213 q = f(t, point_id=1); assert num.allclose(r1, q) 214 q = f(t, point_id=2); assert num.allclose(r2, q) 214 215 215 216 … … 218 219 219 220 timestep = 0 #First timestep 220 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))221 d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))222 d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))223 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 224 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 224 225 225 226 #Reference interpolated values at midpoints on diagonal at … … 231 232 #Let us see if the file function can find the correct 232 233 #values 233 q = f(0, point_id=0); assert allclose(r0, q)234 q = f(0, point_id=1); assert allclose(r1, q)235 q = f(0, point_id=2); assert allclose(r2, q)234 q = f(0, point_id=0); assert num.allclose(r0, q) 235 q = f(0, point_id=1); assert num.allclose(r1, q) 236 q = f(0, point_id=2); assert num.allclose(r2, q) 236 237 237 238 … … 240 241 241 242 timestep = 33 242 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))243 d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))244 d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))245 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 246 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 246 247 247 248 #Reference interpolated values at midpoints on diagonal at … … 251 252 r2 = (D[2] + D[3])/2 252 253 253 q = f(timestep/10., point_id=0); assert allclose(r0, q)254 q = f(timestep/10., point_id=1); assert allclose(r1, q)255 q = f(timestep/10., point_id=2); assert allclose(r2, q)254 q = f(timestep/10., point_id=0); assert num.allclose(r0, q) 255 q = f(timestep/10., point_id=1); assert num.allclose(r1, q) 256 q = f(timestep/10., point_id=2); assert num.allclose(r2, q) 256 257 257 258 … … 261 262 262 263 timestep = 15 263 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))264 d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))265 d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))266 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 267 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 267 268 268 269 #Reference interpolated values at midpoints on diagonal at … … 274 275 # 275 276 timestep = 16 276 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))277 d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))278 d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))279 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 280 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 280 281 281 282 #Reference interpolated values at midpoints on diagonal at … … 290 291 r2 = (r2_0 + r2_1)/2 291 292 292 q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)293 q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)294 q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)293 q = f((timestep - 0.5)/10., point_id=0); assert num.allclose(r0, q) 294 q = f((timestep - 0.5)/10., point_id=1); assert num.allclose(r1, q) 295 q = f((timestep - 0.5)/10., point_id=2); assert num.allclose(r2, q) 295 296 296 297 ################## … … 304 305 305 306 #And the file function gives 306 q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)307 q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)308 q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)307 q = f((timestep - 1.0/3)/10., point_id=0); assert num.allclose(r0, q) 308 q = f((timestep - 1.0/3)/10., point_id=1); assert num.allclose(r1, q) 309 q = f((timestep - 1.0/3)/10., point_id=2); assert num.allclose(r2, q) 309 310 310 311 fid.close() … … 326 327 from shallow_water import Domain, Dirichlet_boundary 327 328 from mesh_factory import rectangular 328 from Numeric import take, concatenate, reshape329 329 330 330 … … 386 386 387 387 last_time_index = len(time)-1 #Last last_time_index 388 d_stage = reshape(take(stage[last_time_index, :], [0,5,10,15]), (4,1))389 d_uh = reshape(take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))390 d_vh = reshape(take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))391 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 391 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 392 392 393 393 #Reference interpolated values at midpoints on diagonal at … … 398 398 399 399 #And the midpoints are found now 400 Dx = take(reshape(x, (16,1)), [0,5,10,15])401 Dy = take(reshape(y, (16,1)), [0,5,10,15])402 403 diag = concatenate( (Dx, Dy), axis=1)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( (Dx, Dy), axis=1) 404 404 d_midpoints = (diag[1:] + diag[:-1])/2 405 405 … … 415 415 416 416 t = time[last_time_index] 417 q = f(t, point_id=0); assert allclose(r0, q)418 q = f(t, point_id=1); assert allclose(r1, q)419 q = f(t, point_id=2); assert allclose(r2, q)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) 420 420 421 421 … … 424 424 425 425 timestep = 0 #First timestep 426 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))427 d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))428 d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))429 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 429 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 430 430 431 431 #Reference interpolated values at midpoints on diagonal at … … 437 437 #Let us see if the file function can find the correct 438 438 #values 439 q = f(0, point_id=0); assert allclose(r0, q)440 q = f(0, point_id=1); assert allclose(r1, q)441 q = f(0, point_id=2); assert allclose(r2, q)439 q = f(0, point_id=0); assert num.allclose(r0, q) 440 q = f(0, point_id=1); assert num.allclose(r1, q) 441 q = f(0, point_id=2); assert num.allclose(r2, q) 442 442 443 443 … … 446 446 447 447 timestep = 33 448 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))449 d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))450 d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))451 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 451 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 452 452 453 453 #Reference interpolated values at midpoints on diagonal at … … 457 457 r2 = (D[2] + D[3])/2 458 458 459 q = f(timestep/10., point_id=0); assert allclose(r0, q)460 q = f(timestep/10., point_id=1); assert allclose(r1, q)461 q = f(timestep/10., point_id=2); assert allclose(r2, q)459 q = f(timestep/10., point_id=0); assert num.allclose(r0, q) 460 q = f(timestep/10., point_id=1); assert num.allclose(r1, q) 461 q = f(timestep/10., point_id=2); assert num.allclose(r2, q) 462 462 463 463 … … 467 467 468 468 timestep = 15 469 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))470 d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))471 d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))472 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 472 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 473 473 474 474 #Reference interpolated values at midpoints on diagonal at … … 480 480 # 481 481 timestep = 16 482 d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))483 d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))484 d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))485 D = concatenate( (d_stage, d_uh, d_vh), axis=1)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)) 485 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 486 486 487 487 #Reference interpolated values at midpoints on diagonal at … … 496 496 r2 = (r2_0 + r2_1)/2 497 497 498 q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)499 q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)500 q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)498 q = f((timestep - 0.5)/10., point_id=0); assert num.allclose(r0, q) 499 q = f((timestep - 0.5)/10., point_id=1); assert num.allclose(r1, q) 500 q = f((timestep - 0.5)/10., point_id=2); assert num.allclose(r2, q) 501 501 502 502 ################## … … 510 510 511 511 #And the file function gives 512 q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)513 q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)514 q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)512 q = f((timestep - 1.0/3)/10., point_id=0); assert num.allclose(r0, q) 513 q = f((timestep - 1.0/3)/10., point_id=1); assert num.allclose(r1, q) 514 q = f((timestep - 1.0/3)/10., point_id=2); assert num.allclose(r2, q) 515 515 516 516 fid.close() … … 542 542 import os, time 543 543 from anuga.config import time_format 544 from Numeric import sin, pi, exp545 544 from mesh_factory import rectangular 546 545 from shallow_water import Domain … … 584 583 domain.set_quantity('xmomentum', f2) 585 584 586 f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)585 f3 = lambda x,y: x**2 + y**2 * num.sin(t*num.pi/600) 587 586 domain.set_quantity('ymomentum', f3) 588 587 … … 604 603 605 604 #Check that FF updates fixes domain starttime 606 assert allclose(domain.starttime, start)605 assert num.allclose(domain.starttime, start) 607 606 608 607 #Check that domain.starttime isn't updated if later … … 611 610 quantities = domain.conserved_quantities, 612 611 interpolation_points = interpolation_points) 613 assert allclose(domain.starttime, start+1)612 assert num.allclose(domain.starttime, start+1) 614 613 domain.starttime = start 615 614 … … 645 644 self.failUnless( q == actual, 'Fail!') 646 645 else: 647 assert allclose(q, actual)646 assert num.allclose(q, actual) 648 647 649 648 … … 655 654 t = 90 #Halfway between 60 and 120 656 655 q = F(t, point_id=id) 657 assert allclose( (q120+q60)/2, q )656 assert num.allclose( (q120+q60)/2, q ) 658 657 659 658 t = 100 #Two thirds of the way between between 60 and 120 660 659 q = F(t, point_id=id) 661 assert allclose(q60/3 + 2*q120/3, q)660 assert num.allclose(q60/3 + 2*q120/3, q) 662 661 663 662 … … 670 669 quantities = domain.conserved_quantities, 671 670 interpolation_points = interpolation_points) 672 assert allclose(domain.starttime, start+delta)671 assert num.allclose(domain.starttime, start+delta) 673 672 674 673 … … 689 688 690 689 q = F(t-delta, point_id=id) 691 assert allclose(q, (k*q1 + (6-k)*q0)/6)690 assert num.allclose(q, (k*q1 + (6-k)*q0)/6) 692 691 693 692 … … 703 702 import os, time 704 703 from anuga.config import time_format 705 from Numeric import sin, pi, exp706 704 from mesh_factory import rectangular 707 705 from shallow_water import Domain … … 751 749 domain.set_quantity('xmomentum', f2) 752 750 753 f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)751 f3 = lambda x,y: x**2 + y**2 * num.sin(t*num.pi/600) 754 752 domain.set_quantity('ymomentum', f3) 755 753 … … 774 772 775 773 #Check that FF updates fixes domain starttime 776 assert allclose(domain.starttime, start)774 assert num.allclose(domain.starttime, start) 777 775 778 776 #Check that domain.starttime isn't updated if later … … 781 779 quantities = domain.conserved_quantities, 782 780 interpolation_points = interpolation_points) 783 assert allclose(domain.starttime, start+1)781 assert num.allclose(domain.starttime, start+1) 784 782 domain.starttime = start 785 783 … … 817 815 self.failUnless( q == actual, 'Fail!') 818 816 else: 819 assert allclose(q, actual)817 assert num.allclose(q, actual) 820 818 821 819 # now lets check points inside the mesh … … 863 861 self.failUnless( q == actual, 'Fail!') 864 862 else: 865 assert allclose(q, actual)863 assert num.allclose(q, actual) 866 864 867 865 … … 873 871 t = 90 #Halfway between 60 and 120 874 872 q = F(t, point_id=id) 875 assert allclose( (q120+q60)/2, q )873 assert num.allclose( (q120+q60)/2, q ) 876 874 877 875 t = 100 #Two thirds of the way between between 60 and 120 878 876 q = F(t, point_id=id) 879 assert allclose(q60/3 + 2*q120/3, q)877 assert num.allclose(q60/3 + 2*q120/3, q)