Changeset 5899
- Timestamp:
- Nov 6, 2008, 12:28:22 PM (17 years ago)
- Location:
- anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes
- Files:
-
- 22 edited
-
domain.py (modified) (10 diffs)
-
ermapper_grids.py (modified) (9 diffs)
-
general_mesh.py (modified) (10 diffs)
-
generic_boundary_conditions.py (modified) (6 diffs)
-
mesh_factory.py (modified) (9 diffs)
-
neighbour_mesh.py (modified) (9 diffs)
-
pmesh2domain.py (modified) (3 diffs)
-
quantity.py (modified) (22 diffs)
-
quantity_ext.c (modified) (1 diff)
-
region.py (modified) (2 diffs)
-
show_balanced_limiters.py (modified) (3 diffs)
-
test_domain.py (modified) (7 diffs)
-
test_ermapper.py (modified) (9 diffs)
-
test_general_mesh.py (modified) (11 diffs)
-
test_generic_boundary_conditions.py (modified) (2 diffs)
-
test_ghost.py (modified) (3 diffs)
-
test_neighbour_mesh.py (modified) (11 diffs)
-
test_pmesh2domain.py (modified) (2 diffs)
-
test_quantity.py (modified) (4 diffs)
-
test_region.py (modified) (8 diffs)
-
test_util.py (modified) (11 diffs)
-
util.py (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/domain.py
r5847 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 """Class Domain - 2D triangular domains for finite-volume computations of 2 4 conservation laws. … … 8 10 """ 9 11 10 from Numeric import allclose, argmax, zeros, Float 12 ##from numpy.oldnumeric import allclose, argmax, zeros, Float 13 from numpy import allclose, argmax, zeros, float 11 14 from anuga.config import epsilon 12 15 from anuga.config import beta_euler, beta_rk2 … … 105 108 106 109 if verbose: print 'Initialising Domain' 107 from Numeric import zeros, Float, Int, ones 110 ## from numpy.oldnumeric import zeros, Float, Int, ones 111 from numpy import zeros, float, int, ones 108 112 from quantity import Quantity 109 113 … … 155 159 for key in self.full_send_dict: 156 160 buffer_shape = self.full_send_dict[key][0].shape[0] 157 self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) , Float))161 self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,float)) 158 162 159 163 160 164 for key in self.ghost_recv_dict: 161 165 buffer_shape = self.ghost_recv_dict[key][0].shape[0] 162 self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) , Float))166 self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,float)) 163 167 164 168 … … 168 172 N = len(self) #number_of_elements 169 173 self.number_of_elements = N 170 self.tri_full_flag = ones(N, Int)174 self.tri_full_flag = ones(N, int) 171 175 for i in self.ghost_recv_dict.keys(): 172 176 for id in self.ghost_recv_dict[i][0]: … … 230 234 # calculation 231 235 N = len(self) # Number_of_triangles 232 self.already_computed_flux = zeros((N, 3), Int)236 self.already_computed_flux = zeros((N, 3), int) 233 237 234 238 # Storage for maximal speeds computed for each triangle by 235 239 # compute_fluxes 236 240 # This is used for diagnostics only (reset at every yieldstep) 237 self.max_speed = zeros(N, Float)241 self.max_speed = zeros(N, float) 238 242 239 243 if mesh_filename is not None: … … 266 270 """ 267 271 268 from Numeric import zeros, Float 272 ## from numpy.oldnumeric import zeros, Float 273 from numpy import zeros, float 269 274 270 275 if not (vertex is None or edge is None): … … 273 278 raise msg 274 279 275 q = zeros( len(self.conserved_quantities), Float)280 q = zeros( len(self.conserved_quantities), float) 276 281 277 282 for i, name in enumerate(self.conserved_quantities): … … 1202 1207 self.number_of_steps = 0 1203 1208 self.number_of_first_order_steps = 0 1204 self.max_speed = zeros(N, Float)1209 self.max_speed = zeros(N, float) 1205 1210 1206 1211 def evolve_one_euler_step(self, yieldstep, finaltime): … … 1568 1573 """ 1569 1574 1570 from Numeric import ones, sum, equal, Float 1575 ## from numpy.oldnumeric import ones, sum, equal, Float 1576 from numpy import ones, sum, equal, float 1571 1577 1572 1578 N = len(self) # Number_of_triangles -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/ermapper_grids.py
r2054 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 3 5 # from os import open, write, read 4 import Numeric 5 6 celltype_map = {'IEEE4ByteReal': Numeric.Float32, 'IEEE8ByteReal': Numeric.Float64} 6 ##import numpy.oldnumeric as Numeric 7 import numpy 8 9 celltype_map = {'IEEE4ByteReal': numpy.float32, 'IEEE8ByteReal': numpy.float64} 7 10 8 11 … … 11 14 write_ermapper_grid(ofile, data, header = {}): 12 15 13 Function to write a 2D Num ericarray to an ERMapper grid. There are a series of conventions adopted within16 Function to write a 2D NumPy array to an ERMapper grid. There are a series of conventions adopted within 14 17 this code, specifically: 15 18 1) The registration coordinate for the data is the SW (or lower-left) corner of the data 16 19 2) The registration coordinates refer to cell centres 17 3) The data is a 2D Num ericarray with the NW-most data in element (0,0) and the SE-most data in element (N,M)20 3) The data is a 2D NumPy array with the NW-most data in element (0,0) and the SE-most data in element (N,M) 18 21 where N is the last line and M is the last column 19 22 4) There has been no testng of the use of a rotated grid. Best to keep data in an NS orientation … … 22 25 ofile: string - filename for output (note the output will consist of two files 23 26 ofile and ofile.ers. Either of these can be entered into this function 24 data: Num eric.array - 2D array containing the data to be output to the grid27 data: NumPy.array - 2D array containing the data to be output to the grid 25 28 header: dictionary - contains spatial information about the grid, in particular: 26 29 header['datum'] datum for the data ('"GDA94"') … … 57 60 58 61 # Check that the data is a 2 dimensional array 59 data_size = Numeric.shape(data)62 data_size = numpy.shape(data) 60 63 assert len(data_size) == 2 61 64 … … 83 86 nrofcellsperlines = int(header['nrofcellsperline']) 84 87 data = read_ermapper_data(data_file) 85 data = Numeric.reshape(data,(nroflines,nrofcellsperlines))88 data = numpy.reshape(data,(nroflines,nrofcellsperlines)) 86 89 return data 87 90 … … 163 166 return header 164 167 165 def write_ermapper_data(grid, ofile, data_format = Numeric.Float32):168 def write_ermapper_data(grid, ofile, data_format = numpy.float32): 166 169 167 170 … … 181 184 182 185 183 # Convert the array to data_format (default format is Float32)186 # Convert the array to data_format (default format is float32) 184 187 grid_as_float = grid.astype(data_format) 185 188 … … 193 196 194 197 195 def read_ermapper_data(ifile, data_format = Numeric.Float32):198 def read_ermapper_data(ifile, data_format = numpy.float32): 196 199 # open input file in a binary format and read the input string 197 200 fid = open(ifile,'rb') … … 199 202 fid.close() 200 203 201 # convert input string to required format (Note default format is Numeric.Float32)202 grid_as_float = Numeric.fromstring(input_string,data_format)204 # convert input string to required format (Note default format is numpy.float32) 205 grid_as_float = numpy.fromstring(input_string,data_format) 203 206 return grid_as_float 204 207 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/general_mesh.py
r5842 r5899 1 2 from Numeric import concatenate, reshape, take, allclose 3 from Numeric import array, zeros, Int, Float, sqrt, sum, arange 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 3 4 ##from numpy.oldnumeric import concatenate, reshape, take, allclose 5 ##from numpy.oldnumeric import array, zeros, Int, Float, sqrt, sum, arange 6 from numpy import concatenate, reshape, take, allclose 7 from numpy import array, zeros, int, float, sqrt, sum, arange 4 8 5 9 from anuga.coordinate_transforms.geo_reference import Geo_reference … … 90 94 if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' 91 95 92 self.triangles = array(triangles, Int) 93 self.nodes = array(nodes, Float) 96 ## self.triangles = array(triangles, Int) 97 ## self.nodes = array(nodes, Float) 98 self.triangles = array(triangles, int) 99 self.nodes = array(nodes, float) 94 100 95 101 … … 138 144 139 145 msg = 'Vertex indices reference non-existing coordinate sets' 140 assert max(self.triangles. flat) < self.nodes.shape[0], msg146 assert max(self.triangles.ravel()) < self.nodes.shape[0], msg 141 147 142 148 … … 146 152 max(self.nodes[:,0]), max(self.nodes[:,1]) ] 147 153 148 self.xy_extent = array(xy_extent, Float) 154 ## self.xy_extent = array(xy_extent, Float) 155 self.xy_extent = array(xy_extent, float) 149 156 150 157 151 158 # 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) 159 ## self.normals = zeros((N, 6), Float) 160 ## self.areas = zeros(N, Float) 161 ## self.edgelengths = zeros((N, 3), Float) 162 self.normals = zeros((N, 6), float) 163 self.areas = zeros(N, float) 164 self.edgelengths = zeros((N, 3), float) 155 165 156 166 # Get x,y coordinates for all triangles and store … … 311 321 i = triangle_id 312 322 msg = 'triangle_id must be an integer' 323 print 'type(triangle_id)=%s. triangle_id=%s' % (type(triangle_id), str(triangle_id)) 313 324 assert int(i) == i, msg 314 325 assert 0 <= i < self.number_of_triangles … … 349 360 350 361 M = self.number_of_triangles 351 vertex_coordinates = zeros((3*M, 2), Float) 362 ## vertex_coordinates = zeros((3*M, 2), Float) 363 vertex_coordinates = zeros((3*M, 2), float) 352 364 353 365 for i in range(M): … … 376 388 indices = range(M) 377 389 378 return take(self.triangles, indices )390 return take(self.triangles, indices, axis=0) 379 391 380 392 … … 409 421 410 422 #T = reshape(array(range(K)).astype(Int), (M,3)) 411 T = reshape(arange(K).astype(Int), (M,3)) # Faster 423 ## T = reshape(arange(K).astype(Int), (M,3)) # Faster 424 T = reshape(arange(K).astype(int), (M,3)) # Faster 412 425 413 426 return T … … 442 455 443 456 # Get number of triangles for this node 444 count = self.number_of_triangles_per_node[node]457 count = int(self.number_of_triangles_per_node[node]) 445 458 446 459 for i in range(count): … … 582 595 Y = C[:,1:6:2].copy() 583 596 584 xmin = min(X. flat)585 xmax = max(X. flat)586 ymin = min(Y. flat)587 ymax = max(Y. flat)597 xmin = min(X.ravel()) 598 xmax = max(X.ravel()) 599 ymin = min(Y.ravel()) 600 ymax = max(Y.ravel()) 588 601 #print "C",C 589 602 return xmin, xmax, ymin, ymax -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
r5672 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 """boundary.py - Classes for implementing boundary conditions 2 4 """ … … 67 69 raise Exception, msg 68 70 69 from Numeric import array, Float 70 self.conserved_quantities=array(conserved_quantities).astype(Float) 71 # from numpy.oldnumeric import array, Float 72 from numpy import array, float 73 self.conserved_quantities=array(conserved_quantities).astype(float) 71 74 72 75 def __repr__(self): … … 97 100 98 101 99 from Numeric import array, Float 102 # from numpy.oldnumeric import array, Float 103 from numpy import array, float 100 104 try: 101 q = array(q).astype( Float)105 q = array(q).astype(float) 102 106 except: 103 107 msg = 'Return value from time boundary function could ' … … 136 140 def __init__(self, filename, domain): 137 141 import time 138 from Numeric import array 142 # from numpy.oldnumeric import array 143 from numpy import array 139 144 from anuga.config import time_format 140 145 from anuga.abstract_2d_finite_volumes.util import File_function … … 206 211 207 212 import time 208 from Numeric import array, zeros, Float 213 # from numpy.oldnumeric import array, zeros, Float 214 from numpy import array, zeros, float 209 215 from anuga.config import time_format 210 216 from anuga.abstract_2d_finite_volumes.util import file_function … … 222 228 223 229 if verbose: print 'Find midpoint coordinates of entire boundary' 224 self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)230 self.midpoint_coordinates = zeros( (len(domain.boundary), 2), float) 225 231 boundary_keys = domain.boundary.keys() 226 232 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/mesh_factory.py
r3678 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 """Library of standard meshes and facilities for reading various 2 4 mesh file formats … … 73 75 74 76 from anuga.config import epsilon 75 from Numeric import zeros, Float, Int 77 # from numpy.oldnumeric import zeros, Float, Int 78 from numpy import zeros, float, int 76 79 77 80 delta1 = float(len1)/m … … 93 96 index = Index(n,m) 94 97 95 points = zeros( (Np,2), Float)98 points = zeros( (Np,2), float) 96 99 97 100 for i in range(m+1): … … 105 108 106 109 107 elements = zeros( (Nt,3), Int)110 elements = zeros( (Nt,3), int) 108 111 boundary = {} 109 112 nt = -1 … … 149 152 150 153 from anuga.config import epsilon 151 from Numeric import zeros, Float, Int 154 # from numpy.oldnumeric import zeros, Float, Int 155 from numpy import zeros, float, int 152 156 153 157 delta1 = float(len1)/m … … 208 212 """ 209 213 210 from Numeric import array 214 # from numpy.oldnumeric import array 215 from numpy import array 211 216 import math 212 217 … … 510 515 """ 511 516 512 from Numeric import array 517 # from numpy.oldnumeric import array 518 from numpy import array 513 519 import math 514 520 … … 592 598 """ 593 599 594 from Numeric import array 600 # from numpy.oldnumeric import array 601 from numpy import array 595 602 import math 596 603 … … 686 693 """ 687 694 688 from Numeric import array 695 # from numpy.oldnumeric import array 696 from numpy import array 689 697 import math 690 698 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
r5866 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 """Classes implementing general 2D triangular mesh with neighbour structure. 2 4 … … 11 13 from anuga.caching import cache 12 14 from math import pi, sqrt 13 from Numeric import array, allclose 15 #from numpy.oldnumeric import array, allclose 16 from numpy import array, allclose 14 17 15 18 … … 81 84 82 85 83 from Numeric import array, zeros, Int, Float, maximum, sqrt, sum 86 # from numpy.oldnumeric import array, zeros, Int, Float, maximum, sqrt, sum 87 from numpy import array, zeros, int, float, maximum, sqrt, sum 84 88 85 89 General_mesh.__init__(self, coordinates, triangles, … … 98 102 99 103 #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)104 self.centroid_coordinates = zeros((N, 2), float) 105 106 self.radii = zeros(N, float) 107 108 self.neighbours = zeros((N, 3), int) 109 self.neighbour_edges = zeros((N, 3), int) 110 self.number_of_boundaries = zeros(N, int) 111 self.surrogate_neighbours = zeros((N, 3), int) 108 112 109 113 #Get x,y coordinates for all triangles and store … … 388 392 self.element_tag is defined 389 393 """ 390 from Numeric import array, Int 394 # from numpy.oldnumeric import array, Int 395 from numpy import array, int 391 396 392 397 if tagged_elements is None: … … 395 400 #Check that all keys in given boundary exist 396 401 for tag in tagged_elements.keys(): 397 tagged_elements[tag] = array(tagged_elements[tag]).astype( Int)402 tagged_elements[tag] = array(tagged_elements[tag]).astype(int) 398 403 399 404 msg = 'Not all elements exist. ' … … 463 468 """ 464 469 465 from Numeric import allclose, sqrt, array, minimum, maximum 470 # from numpy.oldnumeric import allclose, sqrt, array, minimum, maximum 471 from numpy import allclose, sqrt, array, minimum, maximum 466 472 from anuga.utilities.numerical_tools import angle, ensure_numeric 467 473 … … 629 635 from anuga.utilities.numerical_tools import anglediff 630 636 631 from Numeric import sort, allclose 637 # from numpy.oldnumeric import sort, allclose 638 from numpy import sort, allclose 632 639 633 640 N = len(self) … … 805 812 """ 806 813 807 from Numeric import arange 814 # from numpy.oldnumeric import arange 815 from numpy import arange 808 816 from anuga.utilities.numerical_tools import histogram, create_bins 809 817 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/pmesh2domain.py
r4902 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 """Class pmesh2domain - Converting .tsh files to doamains 2 4 … … 161 163 """ 162 164 163 from Numeric import transpose 165 # from numpy.oldnumeric import transpose 166 from numpy import transpose 164 167 from load_mesh.loadASCII import import_mesh_file 165 168 … … 176 179 geo_reference = mesh_dict['geo_reference'] 177 180 if point_atts != None: 178 for quantity, value_vector in map (None, point_titles, point_atts): 181 print 'type(point_atts)=%s' % type(point_atts) 182 for quantity, value_vector in map(None, point_titles, point_atts): 179 183 vertex_quantity_dict[quantity] = value_vector 180 184 tag_dict = pmesh_dict_to_tag_dict(mesh_dict) -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/quantity.py
r5866 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 """Class Quantity - Implements values at each triangular element 2 4 … … 15 17 """ 16 18 17 from Numeric import array, zeros, Float, less, concatenate, NewAxis,\ 18 argmax, argmin, allclose, take, reshape, alltrue 19 #from numpy.oldnumeric import array, zeros, Float, less, concatenate, NewAxis,\ 20 # argmax, argmin, allclose, take, reshape, alltrue, Int 21 from numpy import array, zeros, float, less, concatenate, \ 22 argmax, argmin, allclose, take, reshape, alltrue, int 19 23 20 24 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar … … 38 42 if vertex_values is None: 39 43 N = len(domain) # number_of_elements 40 self.vertex_values = zeros((N, 3), Float)44 self.vertex_values = zeros((N, 3), float) 41 45 else: 42 self.vertex_values = array(vertex_values).astype( Float)46 self.vertex_values = array(vertex_values).astype(float) 43 47 44 48 N, V = self.vertex_values.shape … … 57 61 58 62 # Allocate space for other quantities 59 self.centroid_values = zeros(N, Float)60 self.edge_values = zeros((N, 3), Float)63 self.centroid_values = zeros(N, float) 64 self.edge_values = zeros((N, 3), float) 61 65 62 66 # Allocate space for Gradient 63 self.x_gradient = zeros(N, Float)64 self.y_gradient = zeros(N, Float)67 self.x_gradient = zeros(N, float) 68 self.y_gradient = zeros(N, float) 65 69 66 70 # Allocate space for Limiter Phi 67 self.phi = zeros(N, Float)71 self.phi = zeros(N, float) 68 72 69 73 # Intialise centroid and edge_values … … 72 76 # Allocate space for boundary values 73 77 L = len(domain.boundary) 74 self.boundary_values = zeros(L, Float)78 self.boundary_values = zeros(L, float) 75 79 76 80 # Allocate space for updates of conserved quantities by … … 78 82 79 83 # Allocate space for update fields 80 self.explicit_update = zeros(N, Float )81 self.semi_implicit_update = zeros(N, Float )82 self.centroid_backup_values = zeros(N, Float)84 self.explicit_update = zeros(N, float ) 85 self.semi_implicit_update = zeros(N, float ) 86 self.centroid_backup_values = zeros(N, float) 83 87 84 88 self.set_beta(1.0) … … 382 386 from anuga.geospatial_data.geospatial_data import Geospatial_data 383 387 from types import FloatType, IntType, LongType, ListType, NoneType 384 from Numeric import ArrayType 388 ##NumPy from numpy.oldnumeric import ArrayType 389 from numpy import ndarray, float 385 390 386 391 # Treat special case: Polygon situation … … 448 453 449 454 msg = 'Indices must be a list or None' 450 assert type(indices) in [ListType, NoneType, ArrayType], msg 455 ##NumPy assert type(indices) in [ListType, NoneType, ArrayType], msg 456 assert type(indices) in [ListType, NoneType, ndarray], msg 451 457 452 458 … … 458 464 self.set_values_from_constant(numeric, 459 465 location, indices, verbose) 460 elif type(numeric) in [ArrayType, ListType]: 466 ##NumPy elif type(numeric) in [ArrayType, ListType]: 467 elif type(numeric) in [ndarray, ListType]: 461 468 self.set_values_from_array(numeric, 462 469 location, indices, verbose) … … 474 481 use_cache=use_cache) 475 482 else: 476 msg = 'Illegal type for argument numeric: %s' % str(numeric)477 raise msg483 msg = 'Illegal type for argument numeric: %s' % str(numeric) 484 raise TypeError, msg 478 485 479 486 elif quantity is not None: … … 610 617 """ 611 618 612 from Numeric import array, Float, Int, allclose 613 614 values = array(values).astype(Float) 619 # from numpy.oldnumeric import array, Float, Int, allclose 620 from numpy import array, float, int, allclose 621 622 values = array(values).astype(float) 615 623 616 624 if indices is not None: 617 indices = array(indices).astype( Int)625 indices = array(indices).astype(int) 618 626 msg = 'Number of values must match number of indices:' 619 627 msg += ' You specified %d values and %d indices'\ … … 643 651 'Values array must be 1d' 644 652 645 self.set_vertex_values(values. flat, indices=indices)653 self.set_vertex_values(values.ravel(), indices=indices) 646 654 647 655 else: … … 676 684 A = q.vertex_values 677 685 678 from Numeric import allclose 686 # from numpy.oldnumeric import allclose 687 from numpy import allclose 679 688 msg = 'Quantities are defined on different meshes. '+\ 680 689 'This might be a case for implementing interpolation '+\ … … 717 726 718 727 V = take(self.domain.get_centroid_coordinates(), indices) 728 print 'V=%s' % str(V) 719 729 self.set_values(f(V[:,0], V[:,1]), 720 730 location=location, … … 780 790 781 791 782 points = ensure_numeric(points, Float)783 values = ensure_numeric(values, Float)792 points = ensure_numeric(points, float) 793 values = ensure_numeric(values, float) 784 794 785 795 if location != 'vertices': … … 1118 1128 """ 1119 1129 1120 from Numeric import take 1130 # from numpy.oldnumeric import take 1131 from numpy import take 1121 1132 1122 1133 # FIXME (Ole): I reckon we should have the option of passing a … … 1143 1154 raise msg 1144 1155 1145 import types, Numeric 1146 assert type(indices) in [types.ListType, types.NoneType, 1147 Numeric.ArrayType],\ 1156 # import types, numpy.oldnumeric as Numeric 1157 import types 1158 from numpy import ndarray 1159 1160 assert type(indices) in [types.ListType, types.NoneType, ndarray], \ 1148 1161 'Indices must be a list or None' 1149 1162 … … 1196 1209 """ 1197 1210 1198 from Numeric import array, Float 1211 # from numpy.oldnumeric import array, Float 1212 from numpy import array, float 1199 1213 1200 1214 # Assert that A can be converted to a Numeric array of appropriate dim 1201 A = ensure_numeric(A, Float)1215 A = ensure_numeric(A, float) 1202 1216 1203 1217 # print 'SHAPE A', A.shape … … 1273 1287 """ 1274 1288 1275 from Numeric import concatenate, zeros, Float, Int, array, reshape 1289 # from numpy.oldnumeric import concatenate, zeros, Float, Int, array, reshape 1290 from numpy import concatenate, zeros, float, int, array, reshape 1276 1291 1277 1292 … … 1284 1299 1285 1300 if precision is None: 1286 precision = Float1301 precision = float 1287 1302 1288 1303 … … 1293 1308 V = self.domain.get_triangles() 1294 1309 N = self.domain.number_of_full_nodes # Ignore ghost nodes if any 1295 A = zeros(N, Float)1310 A = zeros(N, float) 1296 1311 points = self.domain.get_nodes() 1297 1312 … … 1343 1358 V = self.domain.get_disconnected_triangles() 1344 1359 points = self.domain.get_vertex_coordinates() 1345 A = self.vertex_values. flat.astype(precision)1360 A = self.vertex_values.ravel().astype(precision) 1346 1361 1347 1362 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/quantity_ext.c
r5306 r5899 11 11 12 12 #include "Python.h" 13 #include " Numeric/arrayobject.h"13 #include "numpy/arrayobject.h" 14 14 #include "math.h" 15 15 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/region.py
r5208 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 """region.py - Classes for implementing region conditions 2 4 … … 8 10 # FIXME (DSG-DSG) add better comments 9 11 10 from Numeric import average12 from numpy.oldnumeric import average 11 13 class Region: 12 14 """Base class for modifying quantities based on a region. -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py
r4204 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 """Example of shallow water wave equation. 2 4 … … 18 20 19 21 from mesh_factory import rectangular 20 from Numeric import array22 from numpy.oldnumeric import array 21 23 22 24 … … 78 80 domain.set_quantity('stage', Z) 79 81 80 from Numeric import allclose82 from numpy.oldnumeric import allclose 81 83 82 84 #Evolve -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_domain.py
r4932 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 … … 6 8 from domain import * 7 9 from anuga.config import epsilon 8 from Numeric import allclose, array, ones, Float 10 ##from numpy.oldnumeric import allclose, array, ones, Float, alltrue 11 from numpy import allclose, array, ones, float, alltrue 9 12 10 13 … … 64 67 65 68 66 assert domain.get_conserved_quantities(0, edge=1) == 0.69 assert alltrue(domain.get_conserved_quantities(0, edge=1) == 0.) 67 70 68 71 … … 536 539 537 540 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 538 denom = ones(4, Float)-domain.timestep*sem541 denom = ones(4, float)-domain.timestep*sem 539 542 540 543 # x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) … … 610 613 from mesh_factory import rectangular 611 614 from shallow_water import Domain 612 from Numeric import zeros, Float 615 ## from numpy.oldnumeric import zeros, Float 616 from numpy import zeros, float 613 617 614 618 #Create basic mesh … … 627 631 from mesh_factory import rectangular 628 632 from shallow_water import Domain 629 from Numeric import zeros, Float 633 ## from numpy.oldnumeric import zeros, Float 634 from numpy import zeros, float 630 635 631 636 #Create basic mesh … … 670 675 from mesh_factory import rectangular 671 676 from shallow_water import Domain 672 from Numeric import zeros, Float 677 ## from numpy.oldnumeric import zeros, Float 678 from numpy import zeros, float 673 679 674 680 #Create basic mesh -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_ermapper.py
r1813 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 … … 5 7 6 8 import ermapper_grids 7 import Numeric 9 ##import numpy.oldnumeric as Numeric 10 import numpy 8 11 from os import remove 9 12 … … 19 22 data_filename = 'test_write_ermapper_grid' 20 23 21 original_grid = Numeric.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]])24 original_grid = numpy.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]]) 22 25 23 26 # Check that the function works when passing the filename without … … 25 28 ermapper_grids.write_ermapper_grid(data_filename, original_grid) 26 29 new_grid = ermapper_grids.read_ermapper_grid(data_filename) 27 assert Numeric.allclose(original_grid, new_grid)30 assert numpy.allclose(original_grid, new_grid) 28 31 29 32 # Check that the function works when passing the filename with … … 31 34 ermapper_grids.write_ermapper_grid(header_filename, original_grid) 32 35 new_grid = ermapper_grids.read_ermapper_grid(header_filename) 33 assert Numeric.allclose(original_grid, new_grid)36 assert numpy.allclose(original_grid, new_grid) 34 37 35 38 # Clean up created files … … 40 43 # Setup test data 41 44 filename = 'test_write_ermapper_grid' 42 original_grid = Numeric.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])45 original_grid = numpy.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0]) 43 46 44 47 # Write test data 45 ermapper_grids.write_ermapper_data(original_grid, filename, Numeric.Float64)48 ermapper_grids.write_ermapper_data(original_grid, filename, numpy.float64) 46 49 47 50 # Read in the test data 48 new_grid = ermapper_grids.read_ermapper_data(filename, Numeric.Float64)51 new_grid = ermapper_grids.read_ermapper_data(filename, numpy.float64) 49 52 50 53 # Check that the test data that has been read in matches the original data 51 assert Numeric.allclose(original_grid, new_grid)54 assert numpy.allclose(original_grid, new_grid) 52 55 53 56 # Clean up created files … … 57 60 # Setup test data 58 61 filename = 'test_write_ermapper_grid' 59 original_grid = Numeric.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0])62 original_grid = numpy.array([0.0, 0.1, 1.0, 2.0, 3.0, 4.0]) 60 63 61 64 # Write test data … … 66 69 67 70 # Check that the test data that has been read in matches the original data 68 assert Numeric.allclose(original_grid, new_grid)71 assert numpy.allclose(original_grid, new_grid) 69 72 70 73 # Clean up created files … … 75 78 76 79 # setup test data 77 original_grid = Numeric.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]])80 original_grid = numpy.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]]) 78 81 # Write test data 79 82 ermapper_grids.write_ermapper_data(original_grid, data_filename) -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_general_mesh.py
r5843 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 … … 7 9 8 10 from anuga.config import epsilon 9 from Numeric import allclose, array, ones, Float 11 ##from numpy.oldnumeric import allclose, array, ones, Float, alltrue 12 from numpy import allclose, array, ones, float, alltrue 10 13 from general_mesh import General_mesh 11 14 from anuga.coordinate_transforms.geo_reference import Geo_reference … … 23 26 def test_get_vertex_coordinates(self): 24 27 from mesh_factory import rectangular 25 from Numeric import zeros, Float 28 ## from numpy.oldnumeric import zeros, Float 29 from numpy import zeros, float 26 30 27 31 #Create basic mesh … … 86 90 """ 87 91 from mesh_factory import rectangular 88 from Numeric import zeros, Float 92 ## from numpy.oldnumeric import zeros, Float 93 from numpy import zeros, float 89 94 90 95 #Create basic mesh … … 114 119 """ 115 120 from mesh_factory import rectangular 116 from Numeric import zeros, Float 121 ## from numpy.oldnumeric import zeros, Float 122 from numpy import zeros, float 117 123 118 124 #Create basic mesh … … 130 136 """ 131 137 from mesh_factory import rectangular 132 from Numeric import zeros, Float, array 138 ## from numpy.oldnumeric import zeros, Float, array 139 from numpy import zeros, float, array 133 140 134 141 a = [0.0, 0.0] … … 204 211 # One node 205 212 L = domain.get_triangles_and_vertices_per_node(node=2) 213 print 'L=%s' % str(L) 206 214 assert allclose(L[0], [0, 2]) 207 215 assert allclose(L[1], [1, 1]) … … 224 232 from mesh_factory import rectangular 225 233 from shallow_water import Domain 226 from Numeric import zeros, Float 234 ## from numpy.oldnumeric import zeros, Float 235 from numpy import zeros, float 227 236 228 237 #Create basic mesh … … 239 248 from mesh_factory import rectangular 240 249 from shallow_water import Domain 241 from Numeric import zeros, Float 250 ## from numpy.oldnumeric import zeros, Float 251 from numpy import zeros, float 242 252 243 253 #Create basic mesh … … 283 293 geo_reference = geo) 284 294 node = domain.get_node(2) 285 self.assert Equal(c, node)295 self.assertTrue(alltrue(c == node)) 286 296 287 297 node = domain.get_node(2, absolute=True) 288 self.assert Equal(nodes_absolute[2], node)298 self.assertTrue(alltrue(nodes_absolute[2] == node)) 289 299 290 300 node = domain.get_node(2, absolute=True) 291 self.assert Equal(nodes_absolute[2], node)301 self.assertTrue(alltrue(nodes_absolute[2] == node)) 292 302 293 303 … … 333 343 if __name__ == "__main__": 334 344 suite = unittest.makeSuite(Test_General_Mesh,'test') 335 #suite = unittest.makeSuite(Test_General_Mesh,'test_get_node')336 345 runner = unittest.TextTestRunner() 337 346 runner.run(suite) -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py
r5666 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 … … 6 8 from generic_boundary_conditions import * 7 9 from anuga.config import epsilon 8 from Numeric import allclose, array 10 #from numpy.oldnumeric import allclose, array 11 from numpy import allclose, array 9 12 10 13 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_ghost.py
r3514 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 … … 6 8 from domain import * 7 9 from anuga.config import epsilon 8 from Numeric import allclose, array, ones, Float10 from numpy.oldnumeric import allclose, array, ones, Float, alltrue 9 11 10 12 … … 43 45 44 46 45 assert domain.get_conserved_quantities(0, edge=1) == 0.47 assert alltrue(domain.get_conserved_quantities(0, edge=1) == 0.) 46 48 47 49 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py
r5288 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 … … 13 15 from mesh_factory import rectangular 14 16 from anuga.config import epsilon 15 from Numeric import allclose, array, Int 17 #from numpy.oldnumeric import allclose, array, Int 18 from numpy import allclose, array, int 16 19 17 20 from anuga.coordinate_transforms.geo_reference import Geo_reference … … 698 701 """ 699 702 from mesh_factory import rectangular 700 from Numeric import zeros, Float 703 # from numpy.oldnumeric import zeros, Float 704 from numpy import zeros, float 701 705 702 706 #Create basic mesh … … 717 721 from mesh_factory import rectangular 718 722 #from mesh import Mesh 719 from Numeric import zeros, Float 723 # from numpy.oldnumeric import zeros, Float 724 from numpy import zeros, float 720 725 721 726 #Create basic mesh … … 736 741 737 742 def test_boundary_polygon_II(self): 738 from Numeric import zeros, Float 743 # from numpy.oldnumeric import zeros, Float 744 from numpy import zeros, float 739 745 740 746 … … 775 781 """ 776 782 777 from Numeric import zeros, Float 783 # from numpy.oldnumeric import zeros, Float 784 from numpy import zeros, float 778 785 779 786 … … 816 823 """ 817 824 818 from Numeric import zeros, Float 825 # from numpy.oldnumeric import zeros, Float 826 from numpy import zeros, float 819 827 820 828 … … 864 872 """ 865 873 866 from Numeric import zeros, Float 874 # from numpy.oldnumeric import zeros, Float 875 from numpy import zeros, float 867 876 from mesh_factory import rectangular 868 877 … … 907 916 908 917 """ 909 from Numeric import zeros, Float 918 # from numpy.oldnumeric import zeros, Float 919 from numpy import zeros, float 910 920 911 921 … … 998 1008 [ 52341.70703125, 38563.39453125]] 999 1009 1000 ##points = ensure_numeric(points, Int)/1000 # Simplify for ease of interpretation1010 ##points = ensure_numeric(points, int)/1000 # Simplify for ease of interpretation 1001 1011 1002 1012 triangles = [[19, 0,15], … … 1101 1111 [ 35406.3359375 , 79332.9140625 ]] 1102 1112 1103 scaled_points = ensure_numeric(points, Int)/1000 # Simplify for ease of interpretation1113 scaled_points = ensure_numeric(points, int)/1000 # Simplify for ease of interpretation 1104 1114 1105 1115 triangles = [[ 0, 1, 2], -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py
r3563 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 # … … 4 6 import unittest 5 7 6 from Numeric import allclose, array 8 #from numpy.oldnumeric import allclose, array 9 from numpy import allclose, array 7 10 8 11 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_quantity.py
r5776 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 … … 7 9 from quantity import * 8 10 from anuga.config import epsilon 9 from Numeric import allclose, array, ones, Float 11 #from numpy.oldnumeric import allclose, array, ones, Float 12 from numpy import allclose, array, ones, float 10 13 11 14 from anuga.fit_interpolate.fit import fit_to_mesh … … 60 63 61 64 62 def test_creation(self):63 64 quantity = Quantity(self.mesh1, [[1,2,3]])65 assert allclose(quantity.vertex_values, [[1.,2.,3.]])66 67 try:68 quantity = Quantity()69 except:70 pass71 else:72 raise 'Should have raised empty quantity exception'73 74 75 try:76 quantity = Quantity([1,2,3])77 except AssertionError:78 pass79 except:80 raise 'Should have raised "mising mesh object" error'81 82 83 def test_creation_zeros(self):84 85 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.]])92 93 94 def test_interpolation(self):95 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]])99 100 101 def test_interpolation2(self):102 quantity = Quantity(self.mesh4,103 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])104 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid105 106 107 quantity.extrapolate_second_order()108 109 #print quantity.vertex_values110 assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],111 [3.+2./3, 6.+2./3, 4.+2./3],112 [4.6, 3.4, 1.],113 [-5.0, 1.0, 4.0]])114 115 #print quantity.edge_values116 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]])120 121 122 #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],123 # [5., 5., 5.],124 # [4.5, 4.5, 0.],125 # [3.0, -1.5, -1.5]])126 65 ## def test_creation(self): 66 ## 67 ## quantity = Quantity(self.mesh1, [[1,2,3]]) 68 ## assert allclose(quantity.vertex_values, [[1.,2.,3.]]) 69 ## 70 ## try: 71 ## quantity = Quantity() 72 ## except: 73 ## pass 74 ## else: 75 ## raise 'Should have raised empty quantity exception' 76 ## 77 ## 78 ## try: 79 ## quantity = Quantity([1,2,3]) 80 ## except AssertionError: 81 ## pass 82 ## except: 83 ## raise 'Should have raised "mising mesh object" error' 84 ## 85 ## 86 ## def test_creation_zeros(self): 87 ## 88 ## quantity = Quantity(self.mesh1) 89 ## assert allclose(quantity.vertex_values, [[0.,0.,0.]]) 90 ## 91 ## 92 ## quantity = Quantity(self.mesh4) 93 ## assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.], 94 ## [0.,0.,0.], [0.,0.,0.]]) 95 ## 96 ## 97 ## def test_interpolation(self): 98 ## quantity = Quantity(self.mesh1, [[1,2,3]]) 99 ## assert allclose(quantity.centroid_values, [2.0]) #Centroid 100 ## 101 ## assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 102 ## 103 ## 104 ## def test_interpolation2(self): 105 ## quantity = Quantity(self.mesh4, 106 ## [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 107 ## assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 108 ## 109 ## 110 ## quantity.extrapolate_second_order() 111 ## 112 ## #print quantity.vertex_values 113 ## assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5], 114 ## [3.+2./3, 6.+2./3, 4.+2./3], 115 ## [4.6, 3.4, 1.], 116 ## [-5.0, 1.0, 4.0]]) 117 ## 118 ## #print quantity.edge_values 119 ## assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25], 120 ## [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6], 121 ## [2.2, 2.8, 4.0], 122 ## [2.5, -0.5, -2.0]]) 123 ## 124 ## 125 ## #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 126 ## # [5., 5., 5.], 127 ## # [4.5, 4.5, 0.], 128 ## # [3.0, -1.5, -1.5]]) 129 ## 127 130 def test_get_extrema_1(self): 128 131 quantity = Quantity(self.mesh4, … … 156 159 157 160 158 def test_get_maximum_2(self): 159 160 a = [0.0, 0.0] 161 b = [0.0, 2.0] 162 c = [2.0,0.0] 163 d = [0.0, 4.0] 164 e = [2.0, 2.0] 165 f = [4.0,0.0] 166 167 points = [a, b, c, d, e, f] 168 #bac, bce, ecf, dbe 169 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 170 171 domain = Domain(points, vertices) 172 173 quantity = Quantity(domain) 174 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 175 176 v = quantity.get_maximum_value() 177 assert v == 6 178 179 v = quantity.get_minimum_value() 180 assert v == 2 181 182 i = quantity.get_maximum_index() 183 assert i == 3 184 185 i = quantity.get_minimum_index() 186 assert i == 0 187 188 x,y = quantity.get_maximum_location() 189 xref, yref = 2.0/3, 8.0/3 190 assert x == xref 191 assert y == yref 192 193 v = quantity.get_values(interpolation_points = [[x,y]]) 194 assert allclose(v, 6) 195 196 x,y = quantity.get_minimum_location() 197 v = quantity.get_values(interpolation_points = [[x,y]]) 198 assert allclose(v, 2) 199 200 #Multiple locations for maximum - 201 #Test that the algorithm picks the first occurrence 202 v = quantity.get_maximum_value(indices=[0,1,2]) 203 assert allclose(v, 4) 204 205 i = quantity.get_maximum_index(indices=[0,1,2]) 206 assert i == 1 207 208 x,y = quantity.get_maximum_location(indices=[0,1,2]) 209 xref, yref = 4.0/3, 4.0/3 210 assert x == xref 211 assert y == yref 212 213 v = quantity.get_values(interpolation_points = [[x,y]]) 214 assert allclose(v, 4) 215 216 # More test of indices...... 217 v = quantity.get_maximum_value(indices=[2,3]) 218 assert allclose(v, 6) 219 220 i = quantity.get_maximum_index(indices=[2,3]) 221 assert i == 3 222 223 x,y = quantity.get_maximum_location(indices=[2,3]) 224 xref, yref = 2.0/3, 8.0/3 225 assert x == xref 226 assert y == yref 227 228 v = quantity.get_values(interpolation_points = [[x,y]]) 229 assert allclose(v, 6) 230 231 232 233 def test_boundary_allocation(self): 234 quantity = Quantity(self.mesh4, 235 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 236 237 assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary) 238 239 240 def test_set_values(self): 241 quantity = Quantity(self.mesh4) 242 243 244 quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], 245 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.]) #Centroid 249 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]]) 253 254 255 # Test default 256 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.]) #Centroid 260 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]]) 264 265 # Test centroids 266 quantity.set_values([1,2,3,4], location = 'centroids') 267 assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid 268 269 # Test exceptions 270 try: 271 quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], 272 location = 'bas kamel tuba') 273 except: 274 pass 275 276 277 try: 278 quantity.set_values([[1,2,3], [0,0,9]]) 279 except AssertionError: 280 pass 281 except: 282 raise 'should have raised Assertionerror' 283 284 285 286 def test_set_values_const(self): 287 quantity = Quantity(self.mesh4) 288 289 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]) #Centroid 294 assert allclose(quantity.edge_values, [[1, 1, 1], 295 [1, 1, 1], 296 [1, 1, 1], 297 [1, 1, 1]]) 298 299 300 quantity.set_values(2.0, location = 'centroids') 301 assert allclose(quantity.centroid_values, [2, 2, 2, 2]) 302 303 304 def test_set_values_func(self): 305 quantity = Quantity(self.mesh4) 306 307 def f(x, y): 308 return x+y 309 310 quantity.set_values(f, location = 'vertices') 311 #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]]) 318 319 320 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 324 325 def test_integral(self): 326 quantity = Quantity(self.mesh4) 327 328 #Try constants first 329 const = 5 330 quantity.set_values(const, location = 'vertices') 331 #print 'Q', quantity.get_integral() 332 333 assert allclose(quantity.get_integral(), self.mesh4.get_area() * const) 334 335 #Try with a linear function 336 def f(x, y): 337 return x+y 338 339 quantity.set_values(f, location = 'vertices') 340 341 342 ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 343 344 assert allclose (quantity.get_integral(), ref_integral) 345 346 347 348 def test_set_vertex_values(self): 349 quantity = Quantity(self.mesh4) 350 quantity.set_vertex_values([0,1,2,3,4,5]) 351 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]) #Centroid 356 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]]) 360 361 362 def test_set_vertex_values_subset(self): 363 quantity = Quantity(self.mesh4) 364 quantity.set_vertex_values([0,1,2,3,4,5]) 365 quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5]) 366 367 assert allclose(quantity.vertex_values, 368 [[1,0,20], [1,20,4], [4,20,50], [30,1,4]]) 369 370 371 def test_set_vertex_values_using_general_interface(self): 372 quantity = Quantity(self.mesh4) 373 374 375 quantity.set_values([0,1,2,3,4,5]) 376 377 378 assert allclose(quantity.vertex_values, 379 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 380 381 #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]]) 388 389 390 391 def test_set_vertex_values_using_general_interface_with_subset(self): 392 """test_set_vertex_values_using_general_interface_with_subset(self): 393 394 Test that indices and polygon works (for constants values) 395 """ 396 397 quantity = Quantity(self.mesh4) 398 399 400 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 404 405 # Constant 406 quantity.set_values(0.0) 407 quantity.set_values(3.14, indices=[0,2], location='vertices') 408 409 # 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]]) 413 414 415 416 # Now try with polygon (pick points where y>2) 417 polygon = [[0,2.1], [4,2.1], [4,7], [0,7]] 418 quantity.set_values(0.0) 419 quantity.set_values(3.14, polygon=polygon) 420 421 assert allclose(quantity.vertex_values, 422 [[0,0,0], [0,0,0], [0,0,0], 423 [3.14,3.14,3.14]]) 424 425 426 # Another polygon (pick triangle 1 and 2 (rightmost triangles) 427 # using centroids 428 polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]] 429 quantity.set_values(0.0) 430 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]]) 436 437 438 # Same polygon now use vertices (default) 439 polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]] 440 quantity.set_values(0.0) 441 #print 'Here 2' 442 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]]) 448 449 450 # Test input checking 451 try: 452 quantity.set_values(3.14, polygon=polygon, indices = [0,2]) 453 except: 454 pass 455 else: 456 msg = 'Should have caught this' 457 raise msg 458 459 460 461 462 463 def test_set_vertex_values_using_general_interface_subset_and_geo(self): 464 """test_set_vertex_values_using_general_interface_with_subset(self): 465 Test that indices and polygon works using georeferencing 466 """ 467 468 quantity = Quantity(self.mesh4) 469 G = Geo_reference(56, 10, 100) 470 quantity.domain.geo_reference = G 471 472 #print quantity.domain.get_nodes(absolute=True) 473 474 475 # Constant 476 quantity.set_values(0.0) 477 quantity.set_values(3.14, indices=[0,2], location='vertices') 478 479 # 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]]) 483 484 485 486 # Now try with polygon (pick points where y>2) 487 polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]]) 488 polygon += [G.xllcorner, G.yllcorner] 489 490 quantity.set_values(0.0) 491 quantity.set_values(3.14, polygon=polygon, location='centroids') 492 493 assert allclose(quantity.vertex_values, 494 [[0,0,0], [0,0,0], [0,0,0], 495 [3.14,3.14,3.14]]) 496 497 498 # 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]]) 500 polygon += [G.xllcorner, G.yllcorner] 501 502 quantity.set_values(0.0) 503 quantity.set_values(3.14, polygon=polygon) 504 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]]) 510 511 512 513 def test_set_values_using_fit(self): 514 515 516 quantity = Quantity(self.mesh4) 517 518 #Get (enough) datapoints 519 data_points = [[ 0.66666667, 0.66666667], 520 [ 1.33333333, 1.33333333], 521 [ 2.66666667, 0.66666667], 522 [ 0.66666667, 2.66666667], 523 [ 0.0, 1.0], 524 [ 0.0, 3.0], 525 [ 1.0, 0.0], 526 [ 1.0, 1.0], 527 [ 1.0, 2.0], 528 [ 1.0, 3.0], 529 [ 2.0, 1.0], 530 [ 3.0, 0.0], 531 [ 3.0, 1.0]] 532 533 z = linear_function(data_points) 534 535 #Use built-in fit_interpolate.fit 536 quantity.set_values( Geospatial_data(data_points, z), alpha = 0 ) 537 #quantity.set_values(points = data_points, values = z, alpha = 0) 538 539 540 answer = linear_function(quantity.domain.get_vertex_coordinates()) 541 #print quantity.vertex_values, answer 542 assert allclose(quantity.vertex_values.flat, answer) 543 544 545 #Now try by setting the same values directly 546 vertex_attributes = fit_to_mesh(data_points, 547 quantity.domain.get_nodes(), 548 quantity.domain.triangles, #FIXME 549 point_attributes=z, 550 alpha = 0, 551 verbose=False) 552 553 #print vertex_attributes 554 quantity.set_values(vertex_attributes) 555 assert allclose(quantity.vertex_values.flat, answer) 556 557 558 559 560 561 def test_test_set_values_using_fit_w_geo(self): 562 563 564 #Mesh 565 vertex_coordinates = [[0.76, 0.76], 566 [0.76, 5.76], 567 [5.76, 0.76]] 568 triangles = [[0,2,1]] 569 570 mesh_georef = Geo_reference(56,-0.76,-0.76) 571 mesh1 = Domain(vertex_coordinates, triangles, 572 geo_reference = mesh_georef) 573 mesh1.check_integrity() 574 575 #Quantity 576 quantity = Quantity(mesh1) 577 578 #Data 579 data_points = [[ 201.0, 401.0], 580 [ 201.0, 403.0], 581 [ 203.0, 401.0]] 582 583 z = [2, 4, 4] 584 585 data_georef = Geo_reference(56,-200,-400) 586 587 588 #Reference 589 ref = fit_to_mesh(data_points, vertex_coordinates, triangles, 590 point_attributes=z, 591 data_origin = data_georef.get_origin(), 592 mesh_origin = mesh_georef.get_origin(), 593 alpha = 0) 594 595 assert allclose( ref, [0,5,5] ) 596 597 598 #Test set_values 599 600 quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 ) 601 602 #quantity.set_values(points = data_points, 603 # values = z, 604 # data_georef = data_georef, 605 # alpha = 0) 606 607 608 #quantity.set_values(points = data_points, 609 # values = z, 610 # data_georef = data_georef, 611 # alpha = 0) 612 assert allclose(quantity.vertex_values.flat, ref) 613 614 615 616 #Test set_values using geospatial data object 617 quantity.vertex_values[:] = 0.0 618 619 geo = Geospatial_data(data_points, z, data_georef) 620 621 622 quantity.set_values(geospatial_data = geo, alpha = 0) 623 assert allclose(quantity.vertex_values.flat, ref) 624 625 626 627 def test_set_values_from_file1(self): 628 quantity = Quantity(self.mesh4) 629 630 #Get (enough) datapoints 631 data_points = [[ 0.66666667, 0.66666667], 632 [ 1.33333333, 1.33333333], 633 [ 2.66666667, 0.66666667], 634 [ 0.66666667, 2.66666667], 635 [ 0.0, 1.0], 636 [ 0.0, 3.0], 637 [ 1.0, 0.0], 638 [ 1.0, 1.0], 639 [ 1.0, 2.0], 640 [ 1.0, 3.0], 641 [ 2.0, 1.0], 642 [ 3.0, 0.0], 643 [ 3.0, 1.0]] 644 645 data_geo_spatial = Geospatial_data(data_points, 646 geo_reference = Geo_reference(56, 0, 0)) 647 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 648 attributes = linear_function(data_points_absolute) 649 att = 'spam_and_eggs' 650 651 #Create .txt file 652 ptsfile = tempfile.mktemp(".txt") 653 file = open(ptsfile,"w") 654 file.write(" x,y," + att + " \n") 655 for data_point, attribute in map(None, data_points_absolute 656 ,attributes): 657 row = str(data_point[0]) + ',' + str(data_point[1]) \ 658 + ',' + str(attribute) 659 file.write(row + "\n") 660 file.close() 661 662 663 #Check that values can be set from file 664 quantity.set_values(filename = ptsfile, 665 attribute_name = att, alpha = 0) 666 answer = linear_function(quantity.domain.get_vertex_coordinates()) 667 668 #print quantity.vertex_values.flat 669 #print answer 670 671 672 assert allclose(quantity.vertex_values.flat, answer) 673 674 675 #Check that values can be set from file using default attribute 676 quantity.set_values(filename = ptsfile, alpha = 0) 677 assert allclose(quantity.vertex_values.flat, answer) 678 679 #Cleanup 680 import os 681 os.remove(ptsfile) 682 683 684 685 def Xtest_set_values_from_file_using_polygon(self): 686 """test_set_values_from_file_using_polygon(self): 687 688 Test that polygon restriction works for general points data 689 """ 690 691 quantity = Quantity(self.mesh4) 692 693 #Get (enough) datapoints 694 data_points = [[ 0.66666667, 0.66666667], 695 [ 1.33333333, 1.33333333], 696 [ 2.66666667, 0.66666667], 697 [ 0.66666667, 2.66666667], 698 [ 0.0, 1.0], 699 [ 0.0, 3.0], 700 [ 1.0, 0.0], 701 [ 1.0, 1.0], 702 [ 1.0, 2.0], 703 [ 1.0, 3.0], 704 [ 2.0, 1.0], 705 [ 3.0, 0.0], 706 [ 3.0, 1.0]] 707 708 data_geo_spatial = Geospatial_data(data_points, 709 geo_reference = Geo_reference(56, 0, 0)) 710 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 711 attributes = linear_function(data_points_absolute) 712 att = 'spam_and_eggs' 713 714 #Create .txt file 715 ptsfile = tempfile.mktemp(".txt") 716 file = open(ptsfile,"w") 717 file.write(" x,y," + att + " \n") 718 for data_point, attribute in map(None, data_points_absolute 719 ,attributes): 720 row = str(data_point[0]) + ',' + str(data_point[1]) \ 721 + ',' + str(attribute) 722 file.write(row + "\n") 723 file.close() 724 725 # Create restricting polygon (containing node #4 (2,2) and 726 # centroid of triangle #1 (bce) 727 polygon = [[1.0, 1.0], [4.0, 1.0], 728 [4.0, 4.0], [1.0, 4.0]] 729 730 #print self.mesh4.nodes 731 #print inside_polygon(self.mesh4.nodes, polygon) 732 assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 733 734 #print quantity.domain.get_vertex_coordinates() 735 #print quantity.domain.get_nodes() 736 737 # Check that values can be set from file 738 quantity.set_values(filename=ptsfile, 739 polygon=polygon, 740 location='unique vertices', 741 alpha=0) 742 743 # Get indices for vertex coordinates in polygon 744 indices = inside_polygon(quantity.domain.get_vertex_coordinates(), 745 polygon) 746 points = take(quantity.domain.get_vertex_coordinates(), indices) 747 748 answer = linear_function(points) 749 750 #print quantity.vertex_values.flat 751 #print answer 752 753 # Check vertices in polygon have been set 754 assert allclose(take(quantity.vertex_values.flat, indices), 755 answer) 756 757 # Check vertices outside polygon are zero 758 indices = outside_polygon(quantity.domain.get_vertex_coordinates(), 759 polygon) 760 assert allclose(take(quantity.vertex_values.flat, indices), 761 0.0) 762 763 #Cleanup 764 import os 765 os.remove(ptsfile) 766 767 768 769 770 def test_cache_test_set_values_from_file(self): 771 # FIXME (Ole): What is this about? 772 # I don't think it checks anything new 773 quantity = Quantity(self.mesh4) 774 775 #Get (enough) datapoints 776 data_points = [[ 0.66666667, 0.66666667], 777 [ 1.33333333, 1.33333333], 778 [ 2.66666667, 0.66666667], 779 [ 0.66666667, 2.66666667], 780 [ 0.0, 1.0], 781 [ 0.0, 3.0], 782 [ 1.0, 0.0], 783 [ 1.0, 1.0], 784 [ 1.0, 2.0], 785 [ 1.0, 3.0], 786 [ 2.0, 1.0], 787 [ 3.0, 0.0], 788 [ 3.0, 1.0]] 789 790 georef = Geo_reference(56, 0, 0) 791 data_geo_spatial = Geospatial_data(data_points, 792 geo_reference=georef) 793 794 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 795 attributes = linear_function(data_points_absolute) 796 att = 'spam_and_eggs' 797 798 # Create .txt file 799 ptsfile = tempfile.mktemp(".txt") 800 file = open(ptsfile,"w") 801 file.write(" x,y," + att + " \n") 802 for data_point, attribute in map(None, data_points_absolute 803 ,attributes): 804 row = str(data_point[0]) + ',' + str(data_point[1]) \ 805 + ',' + str(attribute) 806 file.write(row + "\n") 807 file.close() 808 809 810 # Check that values can be set from file 811 quantity.set_values(filename=ptsfile, 812 attribute_name=att, 813 alpha=0, 814 use_cache=True, 815 verbose=False) 816 answer = linear_function(quantity.domain.get_vertex_coordinates()) 817 assert allclose(quantity.vertex_values.flat, answer) 818 819 820 # Check that values can be set from file using default attribute 821 quantity.set_values(filename=ptsfile, 822 alpha=0) 823 assert allclose(quantity.vertex_values.flat, answer) 824 825 # Check cache 826 quantity.set_values(filename=ptsfile, 827 attribute_name=att, 828 alpha=0, 829 use_cache=True, 830 verbose=False) 831 832 833 #Cleanup 834 import os 835 os.remove(ptsfile) 836 837 def test_set_values_from_lat_long(self): 838 quantity = Quantity(self.mesh_onslow) 839 840 #Get (enough) datapoints 841 data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 842 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]] 843 844 data_geo_spatial = Geospatial_data(data_points, 845 points_are_lats_longs=True) 846 points_UTM = data_geo_spatial.get_data_points(absolute=True) 847 attributes = linear_function(points_UTM) 848 att = 'elevation' 849 850 #Create .txt file 851 txt_file = tempfile.mktemp(".txt") 852 file = open(txt_file,"w") 853 file.write(" lat,long," + att + " \n") 854 for data_point, attribute in map(None, data_points, attributes): 855 row = str(data_point[0]) + ',' + str(data_point[1]) \ 856 + ',' + str(attribute) 857 #print "row", row 858 file.write(row + "\n") 859 file.close() 860 861 862 #Check that values can be set from file 863 quantity.set_values(filename=txt_file, 864 attribute_name=att, 865 alpha=0) 866 answer = linear_function(quantity.domain.get_vertex_coordinates()) 867 868 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 869 #print "answer",answer 870 871 assert allclose(quantity.vertex_values.flat, answer) 872 873 874 #Check that values can be set from file using default attribute 875 quantity.set_values(filename=txt_file, alpha=0) 876 assert allclose(quantity.vertex_values.flat, answer) 877 878 #Cleanup 879 import os 880 os.remove(txt_file) 881 882 def test_set_values_from_lat_long(self): 883 quantity = Quantity(self.mesh_onslow) 884 885 #Get (enough) datapoints 886 data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 887 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]] 888 889 data_geo_spatial = Geospatial_data(data_points, 890 points_are_lats_longs=True) 891 points_UTM = data_geo_spatial.get_data_points(absolute=True) 892 attributes = linear_function(points_UTM) 893 att = 'elevation' 894 895 #Create .txt file 896 txt_file = tempfile.mktemp(".txt") 897 file = open(txt_file,"w") 898 file.write(" lat,long," + att + " \n") 899 for data_point, attribute in map(None, data_points, attributes): 900 row = str(data_point[0]) + ',' + str(data_point[1]) \ 901 + ',' + str(attribute) 902 #print "row", row 903 file.write(row + "\n") 904 file.close() 905 906 907 #Check that values can be set from file 908 quantity.set_values(filename=txt_file, 909 attribute_name=att, alpha=0) 910 answer = linear_function(quantity.domain.get_vertex_coordinates()) 911 912 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 913 #print "answer",answer 914 915 assert allclose(quantity.vertex_values.flat, answer) 916 917 918 #Check that values can be set from file using default attribute 919 quantity.set_values(filename=txt_file, alpha=0) 920 assert allclose(quantity.vertex_values.flat, answer) 921 922 #Cleanup 923 import os 924 os.remove(txt_file) 925 926 def test_set_values_from_UTM_pts(self): 927 quantity = Quantity(self.mesh_onslow) 928 929 #Get (enough) datapoints 930 data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 931 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]] 932 933 data_geo_spatial = Geospatial_data(data_points, 934 points_are_lats_longs=True) 935 points_UTM = data_geo_spatial.get_data_points(absolute=True) 936 attributes = linear_function(points_UTM) 937 att = 'elevation' 938 939 #Create .txt file 940 txt_file = tempfile.mktemp(".txt") 941 file = open(txt_file,"w") 942 file.write(" x,y," + att + " \n") 943 for data_point, attribute in map(None, points_UTM, attributes): 944 row = str(data_point[0]) + ',' + str(data_point[1]) \ 945 + ',' + str(attribute) 946 #print "row", row 947 file.write(row + "\n") 948 file.close() 949 950 951 pts_file = tempfile.mktemp(".pts") 952 convert = Geospatial_data(txt_file) 953 convert.export_points_file(pts_file) 954 955 #Check that values can be set from file 956 quantity.set_values_from_file(pts_file, att, 0, 957 'vertices', None) 958 answer = linear_function(quantity.domain.get_vertex_coordinates()) 959 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 960 #print "answer",answer 961 assert allclose(quantity.vertex_values.flat, answer) 962 963 #Check that values can be set from file 964 quantity.set_values(filename=pts_file, 965 attribute_name=att, alpha=0) 966 answer = linear_function(quantity.domain.get_vertex_coordinates()) 967 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 968 #print "answer",answer 969 assert allclose(quantity.vertex_values.flat, answer) 970 971 972 #Check that values can be set from file using default attribute 973 quantity.set_values(filename=txt_file, alpha=0) 974 assert allclose(quantity.vertex_values.flat, answer) 975 976 #Cleanup 977 import os 978 os.remove(txt_file) 979 os.remove(pts_file) 980 981 def verbose_test_set_values_from_UTM_pts(self): 982 quantity = Quantity(self.mesh_onslow) 983 984 #Get (enough) datapoints 985 data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 986 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 987 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 988 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 989 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 990 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 991 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 992 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 993 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 994 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 995 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 996 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 997 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 998 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 999 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1000 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1001 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1002 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1003 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1004 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1005 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1006 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1007 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1008 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1009 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1010 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1011 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1012 ] 1013 1014 data_geo_spatial = Geospatial_data(data_points, 1015 points_are_lats_longs=True) 1016 points_UTM = data_geo_spatial.get_data_points(absolute=True) 1017 attributes = linear_function(points_UTM) 1018 att = 'elevation' 1019 1020 #Create .txt file 1021 txt_file = tempfile.mktemp(".txt") 1022 file = open(txt_file,"w") 1023 file.write(" x,y," + att + " \n") 1024 for data_point, attribute in map(None, points_UTM, attributes): 1025 row = str(data_point[0]) + ',' + str(data_point[1]) \ 1026 + ',' + str(attribute) 1027 #print "row", row 1028 file.write(row + "\n") 1029 file.close() 1030 1031 1032 pts_file = tempfile.mktemp(".pts") 1033 convert = Geospatial_data(txt_file) 1034 convert.export_points_file(pts_file) 1035 1036 #Check that values can be set from file 1037 quantity.set_values_from_file(pts_file, att, 0, 1038 'vertices', None, verbose = True, 1039 max_read_lines=2) 1040 answer = linear_function(quantity.domain.get_vertex_coordinates()) 1041 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 1042 #print "answer",answer 1043 assert allclose(quantity.vertex_values.flat, answer) 1044 1045 #Check that values can be set from file 1046 quantity.set_values(filename=pts_file, 1047 attribute_name=att, alpha=0) 1048 answer = linear_function(quantity.domain.get_vertex_coordinates()) 1049 #print "quantity.vertex_values.flat", quantity.vertex_values.flat 1050 #print "answer",answer 1051 assert allclose(quantity.vertex_values.flat, answer) 1052 1053 1054 #Check that values can be set from file using default attribute 1055 quantity.set_values(filename=txt_file, alpha=0) 1056 assert allclose(quantity.vertex_values.flat, answer) 1057 1058 #Cleanup 1059 import os 1060 os.remove(txt_file) 1061 os.remove(pts_file) 1062 1063 def test_set_values_from_file_with_georef1(self): 1064 1065 #Mesh in zone 56 (absolute coords) 1066 1067 x0 = 314036.58727982 1068 y0 = 6224951.2960092 1069 1070 a = [x0+0.0, y0+0.0] 1071 b = [x0+0.0, y0+2.0] 1072 c = [x0+2.0, y0+0.0] 1073 d = [x0+0.0, y0+4.0] 1074 e = [x0+2.0, y0+2.0] 1075 f = [x0+4.0, y0+0.0] 1076 1077 points = [a, b, c, d, e, f] 1078 1079 #bac, bce, ecf, dbe 1080 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 1081 1082 #absolute going in .. 1083 mesh4 = Domain(points, elements, 1084 geo_reference = Geo_reference(56, 0, 0)) 1085 mesh4.check_integrity() 1086 quantity = Quantity(mesh4) 1087 1088 #Get (enough) datapoints (relative to georef) 1089 data_points_rel = [[ 0.66666667, 0.66666667], 1090 [ 1.33333333, 1.33333333], 1091 [ 2.66666667, 0.66666667], 1092 [ 0.66666667, 2.66666667], 1093 [ 0.0, 1.0], 1094 [ 0.0, 3.0], 1095 [ 1.0, 0.0], 1096 [ 1.0, 1.0], 1097 [ 1.0, 2.0], 1098 [ 1.0, 3.0], 1099 [ 2.0, 1.0], 1100 [ 3.0, 0.0], 1101 [ 3.0, 1.0]] 1102 1103 data_geo_spatial = Geospatial_data(data_points_rel, 1104 geo_reference = Geo_reference(56, x0, y0)) 1105 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 1106 attributes = linear_function(data_points_absolute) 1107 att = 'spam_and_eggs' 1108 1109 #Create .txt file 1110 ptsfile = tempfile.mktemp(".txt") 1111 file = open(ptsfile,"w") 1112 file.write(" x,y," + att + " \n") 1113 for data_point, attribute in map(None, data_points_absolute 1114 ,attributes): 1115 row = str(data_point[0]) + ',' + str(data_point[1]) \ 1116 + ',' + str(attribute) 1117 file.write(row + "\n") 1118 file.close() 1119 1120 #file = open(ptsfile, 'r') 1121 #lines = file.readlines() 1122 #file.close() 1123 1124 1125 #Check that values can be set from file 1126 quantity.set_values(filename=ptsfile, 1127 attribute_name=att, alpha=0) 1128 answer = linear_function(quantity.domain.get_vertex_coordinates()) 1129 1130 assert allclose(quantity.vertex_values.flat, answer) 1131 1132 1133 #Check that values can be set from file using default attribute 1134 quantity.set_values(filename=ptsfile, alpha=0) 1135 assert allclose(quantity.vertex_values.flat, answer) 1136 1137 #Cleanup 1138 import os 1139 os.remove(ptsfile) 1140 1141 1142 def test_set_values_from_file_with_georef2(self): 1143 1144 #Mesh in zone 56 (relative coords) 1145 1146 x0 = 314036.58727982 1147 y0 = 6224951.2960092 1148 #x0 = 0.0 1149 #y0 = 0.0 1150 1151 a = [0.0, 0.0] 1152 b = [0.0, 2.0] 1153 c = [2.0, 0.0] 1154 d = [0.0, 4.0] 1155 e = [2.0, 2.0] 1156 f = [4.0, 0.0] 1157 1158 points = [a, b, c, d, e, f] 1159 1160 #bac, bce, ecf, dbe 1161 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 1162 1163 mesh4 = Domain(points, elements, 1164 geo_reference = Geo_reference(56, x0, y0)) 1165 mesh4.check_integrity() 1166 quantity = Quantity(mesh4) 1167 1168 #Get (enough) datapoints 1169 data_points = [[ x0+0.66666667, y0+0.66666667], 1170 [ x0+1.33333333, y0+1.33333333], 1171 [ x0+2.66666667, y0+0.66666667], 1172 [ x0+0.66666667, y0+2.66666667], 1173 [ x0+0.0, y0+1.0], 1174 [ x0+0.0, y0+3.0], 1175 [ x0+1.0, y0+0.0], 1176 [ x0+1.0, y0+1.0], 1177 [ x0+1.0, y0+2.0], 1178 [ x0+1.0, y0+3.0], 1179 [ x0+2.0, y0+1.0], 1180 [ x0+3.0, y0+0.0], 1181 [ x0+3.0, y0+1.0]] 1182 1183 1184 data_geo_spatial = Geospatial_data(data_points, 1185 geo_reference = Geo_reference(56, 0, 0)) 1186 data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 1187 attributes = linear_function(data_points_absolute) 1188 att = 'spam_and_eggs' 1189 1190 #Create .txt file 1191 ptsfile = tempfile.mktemp(".txt") 1192 file = open(ptsfile,"w") 1193 file.write(" x,y," + att + " \n") 1194 for data_point, attribute in map(None, data_points_absolute 1195 ,attributes): 1196 row = str(data_point[0]) + ',' + str(data_point[1]) \ 1197 + ',' + str(attribute) 1198 file.write(row + "\n") 1199 file.close() 1200 1201 1202 #Check that values can be set from file 1203 quantity.set_values(filename=ptsfile, 1204 attribute_name=att, alpha=0) 1205 answer = linear_function(quantity.domain. \ 1206 get_vertex_coordinates(absolute=True)) 1207 1208 1209 assert allclose(quantity.vertex_values.flat, answer) 1210 1211 1212 #Check that values can be set from file using default attribute 1213 quantity.set_values(filename=ptsfile, alpha=0) 1214 assert allclose(quantity.vertex_values.flat, answer) 1215 1216 #Cleanup 1217 import os 1218 os.remove(ptsfile) 1219 1220 1221 1222 1223 def test_set_values_from_quantity(self): 1224 1225 quantity1 = Quantity(self.mesh4) 1226 quantity1.set_vertex_values([0,1,2,3,4,5]) 1227 1228 assert allclose(quantity1.vertex_values, 1229 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1230 1231 1232 quantity2 = Quantity(self.mesh4) 1233 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 1237 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 1241 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 1245 1246 #Check detection of quantity as first orgument 1247 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 1251 1252 1253 def Xtest_set_values_from_quantity_using_polygon(self): 1254 """test_set_values_from_quantity_using_polygon(self): 1255 1256 Check that polygon can be used to restrict set_values when 1257 using another quantity as argument. 1258 """ 1259 1260 # Create restricting polygon (containing node #4 (2,2) and 1261 # centroid of triangle #1 (bce) 1262 polygon = [[1.0, 1.0], [4.0, 1.0], 1263 [4.0, 4.0], [1.0, 4.0]] 1264 assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 1265 1266 quantity1 = Quantity(self.mesh4) 1267 quantity1.set_vertex_values([0,1,2,3,4,5]) 1268 1269 assert allclose(quantity1.vertex_values, 1270 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1271 1272 1273 quantity2 = Quantity(self.mesh4) 1274 quantity2.set_values(quantity=quantity1, 1275 polygon=polygon) 1276 1277 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]]), msg 1280 #bac, bce, ecf, dbe 1281 1282 1283 1284 def test_overloading(self): 1285 1286 quantity1 = Quantity(self.mesh4) 1287 quantity1.set_vertex_values([0,1,2,3,4,5]) 1288 1289 assert allclose(quantity1.vertex_values, 1290 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1291 1292 1293 quantity2 = Quantity(self.mesh4) 1294 quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], 1295 location = 'vertices') 1296 1297 1298 1299 quantity3 = Quantity(self.mesh4) 1300 quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]], 1301 location = 'vertices') 1302 1303 1304 # Negation 1305 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) 1309 1310 # Addition 1311 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) 1315 1316 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) 1320 1321 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) 1328 1329 1330 Q = quantity1 + quantity2 - 3 1331 assert allclose(Q.vertex_values, 1332 quantity1.vertex_values + quantity2.vertex_values - 3) 1333 1334 Q = quantity1 - quantity2 1335 assert allclose(Q.vertex_values, 1336 quantity1.vertex_values - quantity2.vertex_values) 1337 1338 #Scaling 1339 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) 1343 Q = 3*quantity1 1344 assert allclose(Q.vertex_values, quantity1.vertex_values*3) 1345 1346 #Multiplication 1347 Q = quantity1 * quantity2 1348 #print Q.vertex_values 1349 #print Q.centroid_values 1350 #print quantity1.centroid_values 1351 #print quantity2.centroid_values 1352 1353 assert allclose(Q.vertex_values, 1354 quantity1.vertex_values * quantity2.vertex_values) 1355 1356 #Linear combinations 1357 Q = 4*quantity1 + 2 1358 assert allclose(Q.vertex_values, 1359 4*quantity1.vertex_values + 2) 1360 1361 Q = quantity1*quantity2 + 2 1362 assert allclose(Q.vertex_values, 1363 quantity1.vertex_values * quantity2.vertex_values + 2) 1364 1365 Q = quantity1*quantity2 + quantity3 1366 assert allclose(Q.vertex_values, 1367 quantity1.vertex_values * quantity2.vertex_values + 1368 quantity3.vertex_values) 1369 Q = quantity1*quantity2 + 3*quantity3 1370 assert allclose(Q.vertex_values, 1371 quantity1.vertex_values * quantity2.vertex_values + 1372 3*quantity3.vertex_values) 1373 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) 1377 1378 Q = quantity1*quantity2 - quantity3 1379 assert allclose(Q.vertex_values, 1380 quantity1.vertex_values * quantity2.vertex_values - 1381 quantity3.vertex_values) 1382 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) 1386 1387 #Try combining quantities and arrays and scalars 1388 Q = 1.5*quantity1*quantity2.vertex_values -\ 1389 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) 1393 1394 1395 #Powers 1396 Q = quantity1**2 1397 assert allclose(Q.vertex_values, quantity1.vertex_values**2) 1398 1399 Q = quantity1**2 +quantity2**2 1400 assert allclose(Q.vertex_values, 1401 quantity1.vertex_values**2 + \ 1402 quantity2.vertex_values**2) 1403 1404 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) 1408 1409 1410 1411 1412 1413 1414 1415 def test_compute_gradient(self): 1416 quantity = Quantity(self.mesh4) 1417 1418 #Set up for a gradient of (2,0) at mid triangle 1419 quantity.set_values([2.0, 4.0, 6.0, 2.0], 1420 location = 'centroids') 1421 1422 1423 #Gradients 1424 quantity.compute_gradients() 1425 1426 a = quantity.x_gradient 1427 b = quantity.y_gradient 1428 #print self.mesh4.centroid_coordinates 1429 #print a, b 1430 1431 #The central triangle (1) 1432 #(using standard gradient based on neigbours controid values) 1433 assert allclose(a[1], 2.0) 1434 assert allclose(b[1], 0.0) 1435 1436 1437 #Left triangle (0) using two point gradient 1438 #q0 = q1 + a*(x0-x1) + b*(y0-y1) <=> 1439 #2 = 4 + a*(-2/3) + b*(-2/3) 1440 assert allclose(a[0] + b[0], 3) 1441 #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 1442 assert allclose(a[0] - b[0], 0) 1443 1444 1445 #Right triangle (2) using two point gradient 1446 #q2 = q1 + a*(x2-x1) + b*(y2-y1) <=> 1447 #6 = 4 + a*(4/3) + b*(-2/3) 1448 assert allclose(2*a[2] - b[2], 3) 1449 #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0) 1450 assert allclose(a[2] + 2*b[2], 0) 1451 1452 1453 #Top triangle (3) using two point gradient 1454 #q3 = q1 + a*(x3-x1) + b*(y3-y1) <=> 1455 #2 = 4 + a*(-2/3) + b*(4/3) 1456 assert allclose(a[3] - 2*b[3], 3) 1457 #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 1458 assert allclose(2*a[3] + b[3], 0) 1459 1460 1461 1462 #print a, b 1463 quantity.extrapolate_second_order() 1464 1465 #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 1469 1470 #a = 1.2, b=-0.6 1471 #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3) 1472 assert allclose(quantity.vertex_values[2,2], 8) 1473 1474 def test_get_gradients(self): 1475 quantity = Quantity(self.mesh4) 1476 1477 #Set up for a gradient of (2,0) at mid triangle 1478 quantity.set_values([2.0, 4.0, 6.0, 2.0], 1479 location = 'centroids') 1480 1481 1482 #Gradients 1483 quantity.compute_gradients() 1484 1485 a, b = quantity.get_gradients() 1486 #print self.mesh4.centroid_coordinates 1487 #print a, b 1488 1489 #The central triangle (1) 1490 #(using standard gradient based on neigbours controid values) 1491 assert allclose(a[1], 2.0) 1492 assert allclose(b[1], 0.0) 1493 1494 1495 #Left triangle (0) using two point gradient 1496 #q0 = q1 + a*(x0-x1) + b*(y0-y1) <=> 1497 #2 = 4 + a*(-2/3) + b*(-2/3) 1498 assert allclose(a[0] + b[0], 3) 1499 #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 1500 assert allclose(a[0] - b[0], 0) 1501 1502 1503 #Right triangle (2) using two point gradient 1504 #q2 = q1 + a*(x2-x1) + b*(y2-y1) <=> 1505 #6 = 4 + a*(4/3) + b*(-2/3) 1506 assert allclose(2*a[2] - b[2], 3) 1507 #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0) 1508 assert allclose(a[2] + 2*b[2], 0) 1509 1510 1511 #Top triangle (3) using two point gradient 1512 #q3 = q1 + a*(x3-x1) + b*(y3-y1) <=> 1513 #2 = 4 + a*(-2/3) + b*(4/3) 1514 assert allclose(a[3] - 2*b[3], 3) 1515 #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 1516 assert allclose(2*a[3] + b[3], 0) 1517 1518 1519 def test_second_order_extrapolation2(self): 1520 quantity = Quantity(self.mesh4) 1521 1522 #Set up for a gradient of (3,1), f(x) = 3x+y 1523 quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3], 1524 location = 'centroids') 1525 1526 #Gradients 1527 quantity.compute_gradients() 1528 1529 a = quantity.x_gradient 1530 b = quantity.y_gradient 1531 1532 #print a, b 1533 1534 assert allclose(a[1], 3.0) 1535 assert allclose(b[1], 1.0) 1536 1537 #Work out the others 1538 1539 quantity.extrapolate_second_order() 1540 1541 #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) 1545 1546 1547 1548 def test_backup_saxpy_centroid_values(self): 1549 quantity = Quantity(self.mesh4) 1550 1551 #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 quantity.set_values(c_values, location = 'centroids') 1555 1556 #Backup 1557 quantity.backup_centroid_values() 1558 1559 #print quantity.vertex_values 1560 assert allclose(quantity.centroid_values, quantity.centroid_backup_values) 1561 1562 1563 quantity.set_values(d_values, location = 'centroids') 1564 1565 quantity.saxpy_centroid_values(2.0, 3.0) 1566 1567 assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values) 1568 1569 1570 1571 def test_first_order_extrapolator(self): 1572 quantity = Quantity(self.mesh4) 1573 1574 #Test centroids 1575 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1576 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1577 1578 #Extrapolate 1579 quantity.extrapolate_first_order() 1580 1581 #Check that gradient is zero 1582 a,b = quantity.get_gradients() 1583 assert allclose(a, [0,0,0,0]) 1584 assert allclose(b, [0,0,0,0]) 1585 1586 #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 1590 1591 def test_second_order_extrapolator(self): 1592 quantity = Quantity(self.mesh4) 1593 1594 #Set up for a gradient of (3,0) at mid triangle 1595 quantity.set_values([2.0, 4.0, 8.0, 2.0], 1596 location = 'centroids') 1597 1598 1599 1600 quantity.extrapolate_second_order() 1601 quantity.limit() 1602 1603 1604 #Assert that central triangle is limited by neighbours 1605 assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 1606 assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1] 1607 1608 assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 1609 assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 1610 1611 assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 1612 assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1] 1613 1614 1615 #Assert that quantities are conserved 1616 from Numeric import sum 1617 for k in range(quantity.centroid_values.shape[0]): 1618 assert allclose (quantity.centroid_values[k], 1619 sum(quantity.vertex_values[k,:])/3) 1620 1621 1622 1623 1624 1625 def test_limit_vertices_by_all_neighbours(self): 1626 quantity = Quantity(self.mesh4) 1627 1628 #Create a deliberate overshoot (e.g. from gradient computation) 1629 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1630 1631 1632 #Limit 1633 quantity.limit_vertices_by_all_neighbours() 1634 1635 #Assert that central triangle is limited by neighbours 1636 assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 1637 assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1] 1638 1639 assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 1640 assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 1641 1642 assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 1643 assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1] 1644 1645 1646 1647 #Assert that quantities are conserved 1648 from Numeric import sum 1649 for k in range(quantity.centroid_values.shape[0]): 1650 assert allclose (quantity.centroid_values[k], 1651 sum(quantity.vertex_values[k,:])/3) 1652 1653 1654 1655 def test_limit_edges_by_all_neighbours(self): 1656 quantity = Quantity(self.mesh4) 1657 1658 #Create a deliberate overshoot (e.g. from gradient computation) 1659 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1660 1661 1662 #Limit 1663 quantity.limit_edges_by_all_neighbours() 1664 1665 #Assert that central triangle is limited by neighbours 1666 assert quantity.edge_values[1,0] <= quantity.centroid_values[2] 1667 assert quantity.edge_values[1,0] >= quantity.centroid_values[0] 1668 1669 assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1670 assert quantity.edge_values[1,1] >= quantity.centroid_values[0] 1671 1672 assert quantity.edge_values[1,2] <= quantity.centroid_values[2] 1673 assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1674 1675 1676 1677 #Assert that quantities are conserved 1678 from Numeric import sum 1679 for k in range(quantity.centroid_values.shape[0]): 1680 assert allclose (quantity.centroid_values[k], 1681 sum(quantity.vertex_values[k,:])/3) 1682 1683 1684 def test_limit_edges_by_neighbour(self): 1685 quantity = Quantity(self.mesh4) 1686 1687 #Create a deliberate overshoot (e.g. from gradient computation) 1688 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1689 1690 1691 #Limit 1692 quantity.limit_edges_by_neighbour() 1693 1694 #Assert that central triangle is limited by neighbours 1695 assert quantity.edge_values[1,0] <= quantity.centroid_values[3] 1696 assert quantity.edge_values[1,0] >= quantity.centroid_values[1] 1697 1698 assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1699 assert quantity.edge_values[1,1] >= quantity.centroid_values[1] 1700 1701 assert quantity.edge_values[1,2] <= quantity.centroid_values[1] 1702 assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1703 1704 1705 1706 #Assert that quantities are conserved 1707 from Numeric import sum 1708 for k in range(quantity.centroid_values.shape[0]): 1709 assert allclose (quantity.centroid_values[k], 1710 sum(quantity.vertex_values[k,:])/3) 1711 1712 def test_limiter2(self): 1713 """Taken from test_shallow_water 1714 """ 1715 quantity = Quantity(self.mesh4) 1716 quantity.domain.beta_w = 0.9 1717 1718 #Test centroids 1719 quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1720 assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1721 1722 1723 #Extrapolate 1724 quantity.extrapolate_second_order() 1725 1726 assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1727 1728 #Limit 1729 quantity.limit() 1730 1731 # limited value for beta_w = 0.9 1732 1733 assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9]) 1734 # limited values for beta_w = 0.5 1735 #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5]) 1736 1737 1738 #Assert that quantities are conserved 1739 from Numeric import sum 1740 for k in range(quantity.centroid_values.shape[0]): 1741 assert allclose (quantity.centroid_values[k], 1742 sum(quantity.vertex_values[k,:])/3) 1743 1744 1745 1746 1747 1748 def test_distribute_first_order(self): 1749 quantity = Quantity(self.mesh4) 1750 1751 #Test centroids 1752 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1753 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1754 1755 1756 #Extrapolate from centroid to vertices and edges 1757 quantity.extrapolate_first_order() 1758 1759 #Interpolate 1760 #quantity.interpolate_from_vertices_to_edges() 1761 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]]) 1766 1767 1768 def test_interpolate_from_vertices_to_edges(self): 1769 quantity = Quantity(self.mesh4) 1770 1771 quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float) 1772 1773 quantity.interpolate_from_vertices_to_edges() 1774 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]]) 1779 1780 1781 def test_interpolate_from_edges_to_vertices(self): 1782 quantity = Quantity(self.mesh4) 1783 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) 1788 1789 quantity.interpolate_from_edges_to_vertices() 1790 1791 assert allclose(quantity.vertex_values, 1792 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1793 1794 1795 1796 def test_distribute_second_order(self): 1797 quantity = Quantity(self.mesh4) 1798 1799 #Test centroids 1800 quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1801 assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1802 1803 1804 #Extrapolate 1805 quantity.extrapolate_second_order() 1806 1807 assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1808 1809 1810 def test_update_explicit(self): 1811 quantity = Quantity(self.mesh4) 1812 1813 #Test centroids 1814 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1815 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1816 1817 #Set explicit_update 1818 quantity.explicit_update = array( [1.,1.,1.,1.] ) 1819 1820 #Update with given timestep 1821 quantity.update(0.1) 1822 1823 x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] ) 1824 assert allclose( quantity.centroid_values, x) 1825 1826 def test_update_semi_implicit(self): 1827 quantity = Quantity(self.mesh4) 1828 1829 #Test centroids 1830 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1831 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1832 1833 #Set semi implicit update 1834 quantity.semi_implicit_update = array([1.,1.,1.,1.]) 1835 1836 #Update with given timestep 1837 timestep = 0.1 1838 quantity.update(timestep) 1839 1840 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 1841 denom = ones(4, Float)-timestep*sem 1842 1843 x = array([1, 2, 3, 4])/denom 1844 assert allclose( quantity.centroid_values, x) 1845 1846 1847 def test_both_updates(self): 1848 quantity = Quantity(self.mesh4) 1849 1850 #Test centroids 1851 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1852 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1853 1854 #Set explicit_update 1855 quantity.explicit_update = array( [4.,3.,2.,1.] ) 1856 1857 #Set semi implicit update 1858 quantity.semi_implicit_update = array( [1.,1.,1.,1.] ) 1859 1860 #Update with given timestep 1861 timestep = 0.1 1862 quantity.update(0.1) 1863 1864 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 1865 denom = ones(4, Float)-timestep*sem 1866 1867 x = array([1., 2., 3., 4.]) 1868 x /= denom 1869 x += timestep*array( [4.0, 3.0, 2.0, 1.0] ) 1870 1871 assert allclose( quantity.centroid_values, x) 1872 1873 1874 1875 1876 #Test smoothing 1877 def test_smoothing(self): 1878 1879 from mesh_factory import rectangular 1880 from shallow_water import Domain, Transmissive_boundary 1881 from Numeric import zeros, Float 1882 from anuga.utilities.numerical_tools import mean 1883 1884 #Create basic mesh 1885 points, vertices, boundary = rectangular(2, 2) 1886 1887 #Create shallow water domain 1888 domain = Domain(points, vertices, boundary) 1889 domain.default_order=2 1890 domain.reduction = mean 1891 1892 1893 #Set some field values 1894 domain.set_quantity('elevation', lambda x,y: x) 1895 domain.set_quantity('friction', 0.03) 1896 1897 1898 ###################### 1899 # Boundary conditions 1900 B = Transmissive_boundary(domain) 1901 domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B}) 1902 1903 1904 ###################### 1905 #Initial condition - with jumps 1906 1907 bed = domain.quantities['elevation'].vertex_values 1908 stage = zeros(bed.shape, Float) 1909 1910 h = 0.03 1911 for i in range(stage.shape[0]): 1912 if i % 2 == 0: 1913 stage[i,:] = bed[i,:] + h 1914 else: 1915 stage[i,:] = bed[i,:] 1916 1917 domain.set_quantity('stage', stage) 1918 1919 stage = domain.quantities['stage'] 1920 1921 #Get smoothed stage 1922 A, V = stage.get_vertex_values(xy=False, smooth=True) 1923 Q = stage.vertex_values 1924 1925 1926 assert A.shape[0] == 9 1927 assert V.shape[0] == 8 1928 assert V.shape[1] == 3 1929 1930 #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) 1935 1936 #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) 1939 1940 1941 #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]) 1950 1951 #Get smoothed stage with XY 1952 X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True) 1953 1954 assert allclose(A, A1) 1955 assert allclose(V, V1) 1956 1957 #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) 1963 1964 1965 1966 1967 def test_vertex_values_no_smoothing(self): 1968 1969 from mesh_factory import rectangular 1970 from shallow_water import Domain, Transmissive_boundary 1971 from Numeric import zeros, Float 1972 from anuga.utilities.numerical_tools import mean 1973 1974 1975 #Create basic mesh 1976 points, vertices, boundary = rectangular(2, 2) 1977 1978 #Create shallow water domain 1979 domain = Domain(points, vertices, boundary) 1980 domain.default_order=2 1981 domain.reduction = mean 1982 1983 1984 #Set some field values 1985 domain.set_quantity('elevation', lambda x,y: x) 1986 domain.set_quantity('friction', 0.03) 1987 1988 1989 ###################### 1990 #Initial condition - with jumps 1991 1992 bed = domain.quantities['elevation'].vertex_values 1993 stage = zeros(bed.shape, Float) 1994 1995 h = 0.03 1996 for i in range(stage.shape[0]): 1997 if i % 2 == 0: 1998 stage[i,:] = bed[i,:] + h 1999 else: 2000 stage[i,:] = bed[i,:] 2001 2002 domain.set_quantity('stage', stage) 2003 2004 #Get stage 2005 stage = domain.quantities['stage'] 2006 A, V = stage.get_vertex_values(xy=False, smooth=False) 2007 Q = stage.vertex_values.flat 2008 2009 for k in range(8): 2010 assert allclose(A[k], Q[k]) 2011 2012 2013 for k in range(8): 2014 assert V[k, 0] == 3*k 2015 assert V[k, 1] == 3*k+1 2016 assert V[k, 2] == 3*k+2 2017 2018 2019 2020 X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False) 2021 2022 2023 assert allclose(A, A1) 2024 assert allclose(V, V1) 2025 2026 #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) 2033 2034 2035 2036 def set_array_values_by_index(self): 2037 2038 from mesh_factory import rectangular 2039 from shallow_water import Domain 2040 from Numeric import zeros, Float 2041 2042 #Create basic mesh 2043 points, vertices, boundary = rectangular(1, 1) 2044 2045 #Create shallow water domain 2046 domain = Domain(points, vertices, boundary) 2047 #print "domain.number_of_elements ",domain.number_of_elements 2048 quantity = Quantity(domain,[[1,1,1],[2,2,2]]) 2049 value = [7] 2050 indices = [1] 2051 quantity.set_array_values_by_index(value, 2052 location = 'centroids', 2053 indices = indices) 2054 #print "quantity.centroid_values",quantity.centroid_values 2055 2056 assert allclose(quantity.centroid_values, [1,7]) 2057 2058 quantity.set_array_values([15,20,25], indices = indices) 2059 assert allclose(quantity.centroid_values, [1,20]) 2060 2061 quantity.set_array_values([15,20,25], indices = indices) 2062 assert allclose(quantity.centroid_values, [1,20]) 2063 2064 def test_setting_some_vertex_values(self): 2065 """ 2066 set values based on triangle lists. 2067 """ 2068 from mesh_factory import rectangular 2069 from shallow_water import Domain 2070 from Numeric import zeros, Float 2071 2072 #Create basic mesh 2073 points, vertices, boundary = rectangular(1, 3) 2074 #print "vertices",vertices 2075 #Create shallow water domain 2076 domain = Domain(points, vertices, boundary) 2077 #print "domain.number_of_elements ",domain.number_of_elements 2078 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2079 [4,4,4],[5,5,5],[6,6,6]]) 2080 2081 2082 # Check that constants work 2083 value = 7 2084 indices = [1] 2085 quantity.set_values(value, 2086 location = 'centroids', 2087 indices = indices) 2088 #print "quantity.centroid_values",quantity.centroid_values 2089 assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2090 2091 value = [7] 2092 indices = [1] 2093 quantity.set_values(value, 2094 location = 'centroids', 2095 indices = indices) 2096 #print "quantity.centroid_values",quantity.centroid_values 2097 assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2098 2099 value = [[15,20,25]] 2100 quantity.set_values(value, indices = indices) 2101 #print "1 quantity.vertex_values",quantity.vertex_values 2102 assert allclose(quantity.vertex_values[1], value[0]) 2103 2104 2105 #print "quantity",quantity.vertex_values 2106 values = [10,100,50] 2107 quantity.set_values(values, indices = [0,1,5], location = 'centroids') 2108 #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]) 2111 #quantity.interpolate() 2112 #print "quantity.centroid_values",quantity.centroid_values 2113 assert allclose(quantity.centroid_values, [10,100,3,4,5,50]) 2114 2115 2116 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2117 [4,4,4],[5,5,5],[6,6,6]]) 2118 values = [10,100,50] 2119 #this will be per unique vertex, indexing the vertices 2120 #print "quantity.vertex_values",quantity.vertex_values 2121 quantity.set_values(values, indices = [0,1,5]) 2122 #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]) 2126 2127 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2128 [4,4,4],[5,5,5],[6,6,6]]) 2129 values = [[31,30,29],[400,400,400],[1000,999,998]] 2130 quantity.set_values(values, indices = [3,3,5]) 2131 quantity.interpolate() 2132 assert allclose(quantity.centroid_values, [1,2,3,400,5,999]) 2133 2134 values = [[1,1,1],[2,2,2],[3,3,3], 2135 [4,4,4],[5,5,5],[6,6,6]] 2136 quantity.set_values(values) 2137 2138 # testing the standard set values by vertex 2139 # indexed by vertex_id in general_mesh.coordinates 2140 values = [0,1,2,3,4,5,6,7] 2141 2142 quantity.set_values(values) 2143 #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.]]) 2150 2151 def test_setting_unique_vertex_values(self): 2152 """ 2153 set values based on unique_vertex lists. 2154 """ 2155 from mesh_factory import rectangular 2156 from shallow_water import Domain 2157 from Numeric import zeros, Float 2158 2159 #Create basic mesh 2160 points, vertices, boundary = rectangular(1, 3) 2161 #print "vertices",vertices 2162 #Create shallow water domain 2163 domain = Domain(points, vertices, boundary) 2164 #print "domain.number_of_elements ",domain.number_of_elements 2165 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2166 [4,4,4],[5,5,5]]) 2167 value = 7 2168 indices = [1,5] 2169 quantity.set_values(value, 2170 location = 'unique vertices', 2171 indices = indices) 2172 #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]) 2176 2177 2178 def test_get_values(self): 2179 """ 2180 get values based on triangle lists. 2181 """ 2182 from mesh_factory import rectangular 2183 from shallow_water import Domain 2184 from Numeric import zeros, Float 2185 2186 #Create basic mesh 2187 points, vertices, boundary = rectangular(1, 3) 2188 2189 #print "points",points 2190 #print "vertices",vertices 2191 #print "boundary",boundary 2192 2193 #Create shallow water domain 2194 domain = Domain(points, vertices, boundary) 2195 #print "domain.number_of_elements ",domain.number_of_elements 2196 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2197 [4,4,4],[5,5,5]]) 2198 2199 #print "quantity.get_values(location = 'unique vertices')", \ 2200 # quantity.get_values(location = 'unique vertices') 2201 2202 #print "quantity.get_values(location = 'unique vertices')", \ 2203 # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2204 # location = 'unique vertices') 2205 2206 answer = [0.5,2,4,5,0,1,3,4.5] 2207 assert allclose(answer, 2208 quantity.get_values(location = 'unique vertices')) 2209 2210 indices = [0,5,3] 2211 answer = [0.5,1,5] 2212 assert allclose(answer, 2213 quantity.get_values(indices=indices, \ 2214 location = 'unique vertices')) 2215 #print "quantity.centroid_values",quantity.centroid_values 2216 #print "quantity.get_values(location = 'centroids') ",\ 2217 # quantity.get_values(location = 'centroids') 2218 2219 2220 2221 2222 def test_get_values_2(self): 2223 """Different mesh (working with domain object) - also check centroids. 2224 """ 2225 2226 2227 a = [0.0, 0.0] 2228 b = [0.0, 2.0] 2229 c = [2.0,0.0] 2230 d = [0.0, 4.0] 2231 e = [2.0, 2.0] 2232 f = [4.0,0.0] 2233 2234 points = [a, b, c, d, e, f] 2235 #bac, bce, ecf, dbe 2236 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2237 2238 domain = Domain(points, vertices) 2239 2240 quantity = Quantity(domain) 2241 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2242 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]]) 2263 2264 # Check averaging over vertices 2265 #a: 0 2266 #b: (4+4+4)/3 2267 #c: (2+2+2)/3 2268 #d: 8 2269 #e: (6+6+6)/3 2270 #f: 4 2271 assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4]) 2272 2273 2274 2275 2276 2277 2278 def test_get_interpolated_values(self): 2279 2280 from mesh_factory import rectangular 2281 from shallow_water import Domain 2282 from Numeric import zeros, Float 2283 2284 #Create basic mesh 2285 points, vertices, boundary = rectangular(1, 3) 2286 domain = Domain(points, vertices, boundary) 2287 2288 #Constant values 2289 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2290 [4,4,4],[5,5,5]]) 2291 2292 2293 2294 # Get interpolated values at centroids 2295 interpolation_points = domain.get_centroid_coordinates() 2296 answer = quantity.get_values(location='centroids') 2297 2298 2299 #print quantity.get_values(points=interpolation_points) 2300 assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points)) 2301 2302 2303 #Arbitrary values 2304 quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7], 2305 [1,4,-9],[2,5,0]]) 2306 2307 2308 # Get interpolated values at centroids 2309 interpolation_points = domain.get_centroid_coordinates() 2310 answer = quantity.get_values(location='centroids') 2311 #print answer 2312 #print quantity.get_values(interpolation_points=interpolation_points) 2313 assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points, 2314 verbose=False)) 2315 2316 2317 #FIXME TODO 2318 #indices = [0,5,3] 2319 #answer = [0.5,1,5] 2320 #assert allclose(answer, 2321 # quantity.get_values(indices=indices, \ 2322 # location = 'unique vertices')) 2323 2324 2325 2326 2327 def test_get_interpolated_values_2(self): 2328 a = [0.0, 0.0] 2329 b = [0.0, 2.0] 2330 c = [2.0,0.0] 2331 d = [0.0, 4.0] 2332 e = [2.0, 2.0] 2333 f = [4.0,0.0] 2334 2335 points = [a, b, c, d, e, f] 2336 #bac, bce, ecf, dbe 2337 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2338 2339 domain = Domain(points, vertices) 2340 2341 quantity = Quantity(domain) 2342 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2343 2344 #First pick one point 2345 x, y = 2.0/3, 8.0/3 2346 v = quantity.get_values(interpolation_points = [[x,y]]) 2347 assert allclose(v, 6) 2348 2349 # Then another to test that algorithm won't blindly 2350 # reuse interpolation matrix 2351 x, y = 4.0/3, 4.0/3 2352 v = quantity.get_values(interpolation_points = [[x,y]]) 2353 assert allclose(v, 4) 2354 2355 2356 2357 def test_get_interpolated_values_with_georef(self): 2358 2359 zone = 56 2360 xllcorner = 308500 2361 yllcorner = 6189000 2362 a = [0.0, 0.0] 2363 b = [0.0, 2.0] 2364 c = [2.0,0.0] 2365 d = [0.0, 4.0] 2366 e = [2.0, 2.0] 2367 f = [4.0,0.0] 2368 2369 points = [a, b, c, d, e, f] 2370 #bac, bce, ecf, dbe 2371 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2372 2373 domain = Domain(points, vertices, 2374 geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2375 2376 quantity = Quantity(domain) 2377 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2378 2379 #First pick one point (and turn it into absolute coordinates) 2380 x, y = 2.0/3, 8.0/3 2381 v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2382 assert allclose(v, 6) 2383 2384 2385 # Then another to test that algorithm won't blindly 2386 # reuse interpolation matrix 2387 x, y = 4.0/3, 4.0/3 2388 v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2389 assert allclose(v, 4) 2390 2391 # Try two points 2392 pts = [[2.0/3 + xllcorner, 8.0/3 + yllcorner], 2393 [4.0/3 + xllcorner, 4.0/3 + yllcorner]] 2394 v = quantity.get_values(interpolation_points=pts) 2395 assert allclose(v, [6, 4]) 2396 2397 # Test it using the geospatial data format with absolute input points and default georef 2398 pts = Geospatial_data(data_points=pts) 2399 v = quantity.get_values(interpolation_points=pts) 2400 assert allclose(v, [6, 4]) 2401 2402 2403 # Test it using the geospatial data format with relative input points 2404 pts = Geospatial_data(data_points=[[2.0/3, 8.0/3], [4.0/3, 4.0/3]], 2405 geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2406 v = quantity.get_values(interpolation_points=pts) 2407 assert allclose(v, [6, 4]) 2408 2409 2410 2411 2412 def test_getting_some_vertex_values(self): 2413 """ 2414 get values based on triangle lists. 2415 """ 2416 from mesh_factory import rectangular 2417 from shallow_water import Domain 2418 from Numeric import zeros, Float 2419 2420 #Create basic mesh 2421 points, vertices, boundary = rectangular(1, 3) 2422 2423 #print "points",points 2424 #print "vertices",vertices 2425 #print "boundary",boundary 2426 2427 #Create shallow water domain 2428 domain = Domain(points, vertices, boundary) 2429 #print "domain.number_of_elements ",domain.number_of_elements 2430 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2431 [4,4,4],[5,5,5],[6,6,6]]) 2432 value = [7] 2433 indices = [1] 2434 quantity.set_values(value, 2435 location = 'centroids', 2436 indices = indices) 2437 #print "quantity.centroid_values",quantity.centroid_values 2438 #print "quantity.get_values(location = 'centroids') ",\ 2439 # quantity.get_values(location = 'centroids') 2440 assert allclose(quantity.centroid_values, 2441 quantity.get_values(location = 'centroids')) 2442 2443 2444 value = [[15,20,25]] 2445 quantity.set_values(value, indices = indices) 2446 #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')) 2451 2452 # get a subset of elements 2453 subset = quantity.get_values(location='centroids', indices=[0,5]) 2454 answer = [quantity.centroid_values[0],quantity.centroid_values[5]] 2455 assert allclose(subset, answer) 2456 2457 2458 subset = quantity.get_values(location='edges', indices=[0,5]) 2459 answer = [quantity.edge_values[0],quantity.edge_values[5]] 2460 #print "subset",subset 2461 #print "answer",answer 2462 assert allclose(subset, answer) 2463 2464 subset = quantity.get_values( indices=[1,5]) 2465 answer = [quantity.vertex_values[1],quantity.vertex_values[5]] 2466 #print "subset",subset 2467 #print "answer",answer 2468 assert allclose(subset, answer) 2469 2470 def test_smooth_vertex_values(self): 2471 """ 2472 get values based on triangle lists. 2473 """ 2474 from mesh_factory import rectangular 2475 from shallow_water import Domain 2476 from Numeric import zeros, Float 2477 2478 #Create basic mesh 2479 points, vertices, boundary = rectangular(2, 2) 2480 2481 #print "points",points 2482 #print "vertices",vertices 2483 #print "boundary",boundary 2484 2485 #Create shallow water domain 2486 domain = Domain(points, vertices, boundary) 2487 #print "domain.number_of_elements ",domain.number_of_elements 2488 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2489 [4,4,4],[5,5,5],[6,6,6],[7,7,7]]) 2490 2491 #print "quantity.get_values(location = 'unique vertices')", \ 2492 # quantity.get_values(location = 'unique vertices') 2493 2494 #print "quantity.get_values(location = 'unique vertices')", \ 2495 # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2496 # location = 'unique vertices') 2497 2498 #print quantity.get_values(location = 'unique vertices') 2499 #print quantity.domain.number_of_triangles_per_node 2500 #print quantity.vertex_values 2501 2502 #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5] 2503 #assert allclose(answer, 2504 # quantity.get_values(location = 'unique vertices')) 2505 2506 quantity.smooth_vertex_values() 2507 2508 #print quantity.vertex_values 2509 2510 2511 answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4], 2512 [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]] 2513 2514 assert allclose(answer_vertex_values, 2515 quantity.vertex_values) 2516 #print "quantity.centroid_values",quantity.centroid_values 2517 #print "quantity.get_values(location = 'centroids') ",\ 2518 # quantity.get_values(location = 'centroids') 161 ## def test_get_maximum_2(self): 162 ## 163 ## a = [0.0, 0.0] 164 ## b = [0.0, 2.0] 165 ## c = [2.0,0.0] 166 ## d = [0.0, 4.0] 167 ## e = [2.0, 2.0] 168 ## f = [4.0,0.0] 169 ## 170 ## points = [a, b, c, d, e, f] 171 ## #bac, bce, ecf, dbe 172 ## vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 173 ## 174 ## domain = Domain(points, vertices) 175 ## 176 ## quantity = Quantity(domain) 177 ## quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 178 ## 179 ## v = quantity.get_maximum_value() 180 ## assert v == 6 181 ## 182 ## v = quantity.get_minimum_value() 183 ## assert v == 2 184 ## 185 ## i = quantity.get_maximum_index() 186 ## assert i == 3 187 ## 188 ## i = quantity.get_minimum_index() 189 ## assert i == 0 190 ## 191 ## x,y = quantity.get_maximum_location() 192 ## xref, yref = 2.0/3, 8.0/3 193 ## assert x == xref 194 ## assert y == yref 195 ## 196 ## v = quantity.get_values(interpolation_points = [[x,y]]) 197 ## assert allclose(v, 6) 198 ## 199 ## x,y = quantity.get_minimum_location() 200 ## v = quantity.get_values(interpolation_points = [[x,y]]) 201 ## assert allclose(v, 2) 202 ## 203 ## #Multiple locations for maximum - 204 ## #Test that the algorithm picks the first occurrence 205 ## v = quantity.get_maximum_value(indices=[0,1,2]) 206 ## assert allclose(v, 4) 207 ## 208 ## i = quantity.get_maximum_index(indices=[0,1,2]) 209 ## assert i == 1 210 ## 211 ## x,y = quantity.get_maximum_location(indices=[0,1,2]) 212 ## xref, yref = 4.0/3, 4.0/3 213 ## assert x == xref 214 ## assert y == yref 215 ## 216 ## v = quantity.get_values(interpolation_points = [[x,y]]) 217 ## assert allclose(v, 4) 218 ## 219 ## # More test of indices...... 220 ## v = quantity.get_maximum_value(indices=[2,3]) 221 ## assert allclose(v, 6) 222 ## 223 ## i = quantity.get_maximum_index(indices=[2,3]) 224 ## assert i == 3 225 ## 226 ## x,y = quantity.get_maximum_location(indices=[2,3]) 227 ## xref, yref = 2.0/3, 8.0/3 228 ## assert x == xref 229 ## assert y == yref 230 ## 231 ## v = quantity.get_values(interpolation_points = [[x,y]]) 232 ## assert allclose(v, 6) 233 ## 234 ## 235 ## 236 ## def test_boundary_allocation(self): 237 ## quantity = Quantity(self.mesh4, 238 ## [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 239 ## 240 ## assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary) 241 ## 242 ## 243 ## def test_set_values(self): 244 ## quantity = Quantity(self.mesh4) 245 ## 246 ## 247 ## quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], 248 ## location = 'vertices') 249 ## assert allclose(quantity.vertex_values, 250 ## [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 251 ## assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 252 ## assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 253 ## [5., 5., 5.], 254 ## [4.5, 4.5, 0.], 255 ## [3.0, -1.5, -1.5]]) 256 ## 257 ## 258 ## # Test default 259 ## quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 260 ## assert allclose(quantity.vertex_values, 261 ## [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 262 ## assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 263 ## assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 264 ## [5., 5., 5.], 265 ## [4.5, 4.5, 0.], 266 ## [3.0, -1.5, -1.5]]) 267 ## 268 ## # Test centroids 269 ## quantity.set_values([1,2,3,4], location = 'centroids') 270 ## assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid 271 ## 272 ## # Test exceptions 273 ## try: 274 ## quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], 275 ## location = 'bas kamel tuba') 276 ## except: 277 ## pass 278 ## 279 ## 280 ## try: 281 ## quantity.set_values([[1,2,3], [0,0,9]]) 282 ## except AssertionError: 283 ## pass 284 ## except: 285 ## raise 'should have raised Assertionerror' 286 ## 287 ## 288 ## 289 ## def test_set_values_const(self): 290 ## quantity = Quantity(self.mesh4) 291 ## 292 ## quantity.set_values(1.0, location = 'vertices') 293 ## assert allclose(quantity.vertex_values, 294 ## [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]]) 295 ## 296 ## assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid 297 ## assert allclose(quantity.edge_values, [[1, 1, 1], 298 ## [1, 1, 1], 299 ## [1, 1, 1], 300 ## [1, 1, 1]]) 301 ## 302 ## 303 ## quantity.set_values(2.0, location = 'centroids') 304 ## assert allclose(quantity.centroid_values, [2, 2, 2, 2]) 305 ## 306 ## 307 ## def test_set_values_func(self): 308 ## quantity = Quantity(self.mesh4) 309 ## 310 ## def f(x, y): 311 ## return x+y 312 ## 313 ## quantity.set_values(f, location = 'vertices') 314 ## #print "quantity.vertex_values",quantity.vertex_values 315 ## assert allclose(quantity.vertex_values, 316 ## [[2,0,2], [2,2,4], [4,2,4], [4,2,4]]) 317 ## assert allclose(quantity.centroid_values, 318 ## [4.0/3, 8.0/3, 10.0/3, 10.0/3]) 319 ## assert allclose(quantity.edge_values, 320 ## [[1,2,1], [3,3,2], [3,4,3], [3,4,3]]) 321 ## 322 ## 323 ## quantity.set_values(f, location = 'centroids') 324 ## assert allclose(quantity.centroid_values, 325 ## [4.0/3, 8.0/3, 10.0/3, 10.0/3]) 326 ## 327 ## 328 ## def test_integral(self): 329 ## quantity = Quantity(self.mesh4) 330 ## 331 ## #Try constants first 332 ## const = 5 333 ## quantity.set_values(const, location = 'vertices') 334 ## #print 'Q', quantity.get_integral() 335 ## 336 ## assert allclose(quantity.get_integral(), self.mesh4.get_area() * const) 337 ## 338 ## #Try with a linear function 339 ## def f(x, y): 340 ## return x+y 341 ## 342 ## quantity.set_values(f, location = 'vertices') 343 ## 344 ## 345 ## ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 346 ## 347 ## assert allclose (quantity.get_integral(), ref_integral) 348 ## 349 ## 350 ## 351 ## def test_set_vertex_values(self): 352 ## quantity = Quantity(self.mesh4) 353 ## quantity.set_vertex_values([0,1,2,3,4,5]) 354 ## 355 ## assert allclose(quantity.vertex_values, 356 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 357 ## assert allclose(quantity.centroid_values, 358 ## [1., 7./3, 11./3, 8./3]) #Centroid 359 ## assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 360 ## [3., 2.5, 1.5], 361 ## [3.5, 4.5, 3.], 362 ## [2.5, 3.5, 2]]) 363 ## 364 ## 365 ## def test_set_vertex_values_subset(self): 366 ## quantity = Quantity(self.mesh4) 367 ## quantity.set_vertex_values([0,1,2,3,4,5]) 368 ## quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5]) 369 ## 370 ## assert allclose(quantity.vertex_values, 371 ## [[1,0,20], [1,20,4], [4,20,50], [30,1,4]]) 372 ## 373 ## 374 ## def test_set_vertex_values_using_general_interface(self): 375 ## quantity = Quantity(self.mesh4) 376 ## 377 ## 378 ## quantity.set_values([0,1,2,3,4,5]) 379 ## 380 ## 381 ## assert allclose(quantity.vertex_values, 382 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 383 ## 384 ## #Centroid 385 ## assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) 386 ## 387 ## assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 388 ## [3., 2.5, 1.5], 389 ## [3.5, 4.5, 3.], 390 ## [2.5, 3.5, 2]]) 391 ## 392 ## 393 ## 394 ## def test_set_vertex_values_using_general_interface_with_subset(self): 395 ## """test_set_vertex_values_using_general_interface_with_subset(self): 396 ## 397 ## Test that indices and polygon works (for constants values) 398 ## """ 399 ## 400 ## quantity = Quantity(self.mesh4) 401 ## 402 ## 403 ## quantity.set_values([0,2,3,5], indices=[0,2,3,5]) 404 ## assert allclose(quantity.vertex_values, 405 ## [[0,0,2], [0,2,0], [0,2,5], [3,0,0]]) 406 ## 407 ## 408 ## # Constant 409 ## quantity.set_values(0.0) 410 ## quantity.set_values(3.14, indices=[0,2], location='vertices') 411 ## 412 ## # Indices refer to triangle numbers 413 ## assert allclose(quantity.vertex_values, 414 ## [[3.14,3.14,3.14], [0,0,0], 415 ## [3.14,3.14,3.14], [0,0,0]]) 416 ## 417 ## 418 ## 419 ## # Now try with polygon (pick points where y>2) 420 ## polygon = [[0,2.1], [4,2.1], [4,7], [0,7]] 421 ## quantity.set_values(0.0) 422 ## quantity.set_values(3.14, polygon=polygon) 423 ## 424 ## assert allclose(quantity.vertex_values, 425 ## [[0,0,0], [0,0,0], [0,0,0], 426 ## [3.14,3.14,3.14]]) 427 ## 428 ## 429 ## # Another polygon (pick triangle 1 and 2 (rightmost triangles) 430 ## # using centroids 431 ## polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]] 432 ## quantity.set_values(0.0) 433 ## quantity.set_values(3.14, location='centroids', polygon=polygon) 434 ## assert allclose(quantity.vertex_values, 435 ## [[0,0,0], 436 ## [3.14,3.14,3.14], 437 ## [3.14,3.14,3.14], 438 ## [0,0,0]]) 439 ## 440 ## 441 ## # Same polygon now use vertices (default) 442 ## polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]] 443 ## quantity.set_values(0.0) 444 ## #print 'Here 2' 445 ## quantity.set_values(3.14, polygon=polygon) 446 ## assert allclose(quantity.vertex_values, 447 ## [[0,0,0], 448 ## [3.14,3.14,3.14], 449 ## [3.14,3.14,3.14], 450 ## [0,0,0]]) 451 ## 452 ## 453 ## # Test input checking 454 ## try: 455 ## quantity.set_values(3.14, polygon=polygon, indices = [0,2]) 456 ## except: 457 ## pass 458 ## else: 459 ## msg = 'Should have caught this' 460 ## raise msg 461 ## 462 ## 463 ## 464 ## 465 ## 466 ## def test_set_vertex_values_using_general_interface_subset_and_geo(self): 467 ## """test_set_vertex_values_using_general_interface_with_subset(self): 468 ## Test that indices and polygon works using georeferencing 469 ## """ 470 ## 471 ## quantity = Quantity(self.mesh4) 472 ## G = Geo_reference(56, 10, 100) 473 ## quantity.domain.geo_reference = G 474 ## 475 ## #print quantity.domain.get_nodes(absolute=True) 476 ## 477 ## 478 ## # Constant 479 ## quantity.set_values(0.0) 480 ## quantity.set_values(3.14, indices=[0,2], location='vertices') 481 ## 482 ## # Indices refer to triangle numbers here - not vertices (why?) 483 ## assert allclose(quantity.vertex_values, 484 ## [[3.14,3.14,3.14], [0,0,0], 485 ## [3.14,3.14,3.14], [0,0,0]]) 486 ## 487 ## 488 ## 489 ## # Now try with polygon (pick points where y>2) 490 ## polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]]) 491 ## polygon += [G.xllcorner, G.yllcorner] 492 ## 493 ## quantity.set_values(0.0) 494 ## quantity.set_values(3.14, polygon=polygon, location='centroids') 495 ## 496 ## assert allclose(quantity.vertex_values, 497 ## [[0,0,0], [0,0,0], [0,0,0], 498 ## [3.14,3.14,3.14]]) 499 ## 500 ## 501 ## # Another polygon (pick triangle 1 and 2 (rightmost triangles) 502 ## polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]) 503 ## polygon += [G.xllcorner, G.yllcorner] 504 ## 505 ## quantity.set_values(0.0) 506 ## quantity.set_values(3.14, polygon=polygon) 507 ## 508 ## assert allclose(quantity.vertex_values, 509 ## [[0,0,0], 510 ## [3.14,3.14,3.14], 511 ## [3.14,3.14,3.14], 512 ## [0,0,0]]) 513 ## 514 ## 515 ## 516 ## def test_set_values_using_fit(self): 517 ## 518 ## 519 ## quantity = Quantity(self.mesh4) 520 ## 521 ## #Get (enough) datapoints 522 ## data_points = [[ 0.66666667, 0.66666667], 523 ## [ 1.33333333, 1.33333333], 524 ## [ 2.66666667, 0.66666667], 525 ## [ 0.66666667, 2.66666667], 526 ## [ 0.0, 1.0], 527 ## [ 0.0, 3.0], 528 ## [ 1.0, 0.0], 529 ## [ 1.0, 1.0], 530 ## [ 1.0, 2.0], 531 ## [ 1.0, 3.0], 532 ## [ 2.0, 1.0], 533 ## [ 3.0, 0.0], 534 ## [ 3.0, 1.0]] 535 ## 536 ## z = linear_function(data_points) 537 ## 538 ## #Use built-in fit_interpolate.fit 539 ## quantity.set_values( Geospatial_data(data_points, z), alpha = 0 ) 540 ## #quantity.set_values(points = data_points, values = z, alpha = 0) 541 ## 542 ## 543 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 544 ## #print quantity.vertex_values, answer 545 ## assert allclose(quantity.vertex_values.ravel(), answer) 546 ## 547 ## 548 ## #Now try by setting the same values directly 549 ## vertex_attributes = fit_to_mesh(data_points, 550 ## quantity.domain.get_nodes(), 551 ## quantity.domain.triangles, #FIXME 552 ## point_attributes=z, 553 ## alpha = 0, 554 ## verbose=False) 555 ## 556 ## #print vertex_attributes 557 ## quantity.set_values(vertex_attributes) 558 ## assert allclose(quantity.vertex_values.ravel(), answer) 559 ## 560 ## 561 ## 562 ## 563 ## 564 ## def test_test_set_values_using_fit_w_geo(self): 565 ## 566 ## 567 ## #Mesh 568 ## vertex_coordinates = [[0.76, 0.76], 569 ## [0.76, 5.76], 570 ## [5.76, 0.76]] 571 ## triangles = [[0,2,1]] 572 ## 573 ## mesh_georef = Geo_reference(56,-0.76,-0.76) 574 ## mesh1 = Domain(vertex_coordinates, triangles, 575 ## geo_reference = mesh_georef) 576 ## mesh1.check_integrity() 577 ## 578 ## #Quantity 579 ## quantity = Quantity(mesh1) 580 ## 581 ## #Data 582 ## data_points = [[ 201.0, 401.0], 583 ## [ 201.0, 403.0], 584 ## [ 203.0, 401.0]] 585 ## 586 ## z = [2, 4, 4] 587 ## 588 ## data_georef = Geo_reference(56,-200,-400) 589 ## 590 ## 591 ## #Reference 592 ## ref = fit_to_mesh(data_points, vertex_coordinates, triangles, 593 ## point_attributes=z, 594 ## data_origin = data_georef.get_origin(), 595 ## mesh_origin = mesh_georef.get_origin(), 596 ## alpha = 0) 597 ## 598 ## assert allclose( ref, [0,5,5] ) 599 ## 600 ## 601 ## #Test set_values 602 ## 603 ## quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 ) 604 ## 605 ## #quantity.set_values(points = data_points, 606 ## # values = z, 607 ## # data_georef = data_georef, 608 ## # alpha = 0) 609 ## 610 ## 611 ## #quantity.set_values(points = data_points, 612 ## # values = z, 613 ## # data_georef = data_georef, 614 ## # alpha = 0) 615 ## assert allclose(quantity.vertex_values.ravel(), ref) 616 ## 617 ## 618 ## 619 ## #Test set_values using geospatial data object 620 ## quantity.vertex_values[:] = 0.0 621 ## 622 ## geo = Geospatial_data(data_points, z, data_georef) 623 ## 624 ## 625 ## quantity.set_values(geospatial_data = geo, alpha = 0) 626 ## assert allclose(quantity.vertex_values.ravel(), ref) 627 ## 628 ## 629 ## 630 ## def test_set_values_from_file1(self): 631 ## quantity = Quantity(self.mesh4) 632 ## 633 ## #Get (enough) datapoints 634 ## data_points = [[ 0.66666667, 0.66666667], 635 ## [ 1.33333333, 1.33333333], 636 ## [ 2.66666667, 0.66666667], 637 ## [ 0.66666667, 2.66666667], 638 ## [ 0.0, 1.0], 639 ## [ 0.0, 3.0], 640 ## [ 1.0, 0.0], 641 ## [ 1.0, 1.0], 642 ## [ 1.0, 2.0], 643 ## [ 1.0, 3.0], 644 ## [ 2.0, 1.0], 645 ## [ 3.0, 0.0], 646 ## [ 3.0, 1.0]] 647 ## 648 ## data_geo_spatial = Geospatial_data(data_points, 649 ## geo_reference = Geo_reference(56, 0, 0)) 650 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 651 ## attributes = linear_function(data_points_absolute) 652 ## att = 'spam_and_eggs' 653 ## 654 ## #Create .txt file 655 ## ptsfile = tempfile.mktemp(".txt") 656 ## file = open(ptsfile,"w") 657 ## file.write(" x,y," + att + " \n") 658 ## for data_point, attribute in map(None, data_points_absolute 659 ## ,attributes): 660 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 661 ## + ',' + str(attribute) 662 ## file.write(row + "\n") 663 ## file.close() 664 ## 665 ## 666 ## #Check that values can be set from file 667 ## quantity.set_values(filename = ptsfile, 668 ## attribute_name = att, alpha = 0) 669 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 670 ## 671 ## #print quantity.vertex_values.ravel() 672 ## #print answer 673 ## 674 ## 675 ## assert allclose(quantity.vertex_values.ravel(), answer) 676 ## 677 ## 678 ## #Check that values can be set from file using default attribute 679 ## quantity.set_values(filename = ptsfile, alpha = 0) 680 ## assert allclose(quantity.vertex_values.ravel(), answer) 681 ## 682 ## #Cleanup 683 ## import os 684 ## os.remove(ptsfile) 685 ## 686 ## 687 ## 688 ## def Xtest_set_values_from_file_using_polygon(self): 689 ## """test_set_values_from_file_using_polygon(self): 690 ## 691 ## Test that polygon restriction works for general points data 692 ## """ 693 ## 694 ## quantity = Quantity(self.mesh4) 695 ## 696 ## #Get (enough) datapoints 697 ## data_points = [[ 0.66666667, 0.66666667], 698 ## [ 1.33333333, 1.33333333], 699 ## [ 2.66666667, 0.66666667], 700 ## [ 0.66666667, 2.66666667], 701 ## [ 0.0, 1.0], 702 ## [ 0.0, 3.0], 703 ## [ 1.0, 0.0], 704 ## [ 1.0, 1.0], 705 ## [ 1.0, 2.0], 706 ## [ 1.0, 3.0], 707 ## [ 2.0, 1.0], 708 ## [ 3.0, 0.0], 709 ## [ 3.0, 1.0]] 710 ## 711 ## data_geo_spatial = Geospatial_data(data_points, 712 ## geo_reference = Geo_reference(56, 0, 0)) 713 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 714 ## attributes = linear_function(data_points_absolute) 715 ## att = 'spam_and_eggs' 716 ## 717 ## #Create .txt file 718 ## ptsfile = tempfile.mktemp(".txt") 719 ## file = open(ptsfile,"w") 720 ## file.write(" x,y," + att + " \n") 721 ## for data_point, attribute in map(None, data_points_absolute 722 ## ,attributes): 723 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 724 ## + ',' + str(attribute) 725 ## file.write(row + "\n") 726 ## file.close() 727 ## 728 ## # Create restricting polygon (containing node #4 (2,2) and 729 ## # centroid of triangle #1 (bce) 730 ## polygon = [[1.0, 1.0], [4.0, 1.0], 731 ## [4.0, 4.0], [1.0, 4.0]] 732 ## 733 ## #print self.mesh4.nodes 734 ## #print inside_polygon(self.mesh4.nodes, polygon) 735 ## assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 736 ## 737 ## #print quantity.domain.get_vertex_coordinates() 738 ## #print quantity.domain.get_nodes() 739 ## 740 ## # Check that values can be set from file 741 ## quantity.set_values(filename=ptsfile, 742 ## polygon=polygon, 743 ## location='unique vertices', 744 ## alpha=0) 745 ## 746 ## # Get indices for vertex coordinates in polygon 747 ## indices = inside_polygon(quantity.domain.get_vertex_coordinates(), 748 ## polygon) 749 ## points = take(quantity.domain.get_vertex_coordinates(), indices) 750 ## 751 ## answer = linear_function(points) 752 ## 753 ## #print quantity.vertex_values.ravel() 754 ## #print answer 755 ## 756 ## # Check vertices in polygon have been set 757 ## assert allclose(take(quantity.vertex_values.ravel(), indices), 758 ## answer) 759 ## 760 ## # Check vertices outside polygon are zero 761 ## indices = outside_polygon(quantity.domain.get_vertex_coordinates(), 762 ## polygon) 763 ## assert allclose(take(quantity.vertex_values.ravel(), indices), 764 ## 0.0) 765 ## 766 ## #Cleanup 767 ## import os 768 ## os.remove(ptsfile) 769 ## 770 ## 771 ## 772 ## 773 ## def test_cache_test_set_values_from_file(self): 774 ## # FIXME (Ole): What is this about? 775 ## # I don't think it checks anything new 776 ## quantity = Quantity(self.mesh4) 777 ## 778 ## #Get (enough) datapoints 779 ## data_points = [[ 0.66666667, 0.66666667], 780 ## [ 1.33333333, 1.33333333], 781 ## [ 2.66666667, 0.66666667], 782 ## [ 0.66666667, 2.66666667], 783 ## [ 0.0, 1.0], 784 ## [ 0.0, 3.0], 785 ## [ 1.0, 0.0], 786 ## [ 1.0, 1.0], 787 ## [ 1.0, 2.0], 788 ## [ 1.0, 3.0], 789 ## [ 2.0, 1.0], 790 ## [ 3.0, 0.0], 791 ## [ 3.0, 1.0]] 792 ## 793 ## georef = Geo_reference(56, 0, 0) 794 ## data_geo_spatial = Geospatial_data(data_points, 795 ## geo_reference=georef) 796 ## 797 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 798 ## attributes = linear_function(data_points_absolute) 799 ## att = 'spam_and_eggs' 800 ## 801 ## # Create .txt file 802 ## ptsfile = tempfile.mktemp(".txt") 803 ## file = open(ptsfile,"w") 804 ## file.write(" x,y," + att + " \n") 805 ## for data_point, attribute in map(None, data_points_absolute 806 ## ,attributes): 807 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 808 ## + ',' + str(attribute) 809 ## file.write(row + "\n") 810 ## file.close() 811 ## 812 ## 813 ## # Check that values can be set from file 814 ## quantity.set_values(filename=ptsfile, 815 ## attribute_name=att, 816 ## alpha=0, 817 ## use_cache=True, 818 ## verbose=False) 819 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 820 ## assert allclose(quantity.vertex_values.ravel(), answer) 821 ## 822 ## 823 ## # Check that values can be set from file using default attribute 824 ## quantity.set_values(filename=ptsfile, 825 ## alpha=0) 826 ## assert allclose(quantity.vertex_values.ravel(), answer) 827 ## 828 ## # Check cache 829 ## quantity.set_values(filename=ptsfile, 830 ## attribute_name=att, 831 ## alpha=0, 832 ## use_cache=True, 833 ## verbose=False) 834 ## 835 ## 836 ## #Cleanup 837 ## import os 838 ## os.remove(ptsfile) 839 ## 840 ## def test_set_values_from_lat_long(self): 841 ## quantity = Quantity(self.mesh_onslow) 842 ## 843 ## #Get (enough) datapoints 844 ## data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 845 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]] 846 ## 847 ## data_geo_spatial = Geospatial_data(data_points, 848 ## points_are_lats_longs=True) 849 ## points_UTM = data_geo_spatial.get_data_points(absolute=True) 850 ## attributes = linear_function(points_UTM) 851 ## att = 'elevation' 852 ## 853 ## #Create .txt file 854 ## txt_file = tempfile.mktemp(".txt") 855 ## file = open(txt_file,"w") 856 ## file.write(" lat,long," + att + " \n") 857 ## for data_point, attribute in map(None, data_points, attributes): 858 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 859 ## + ',' + str(attribute) 860 ## #print "row", row 861 ## file.write(row + "\n") 862 ## file.close() 863 ## 864 ## 865 ## #Check that values can be set from file 866 ## quantity.set_values(filename=txt_file, 867 ## attribute_name=att, 868 ## alpha=0) 869 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 870 ## 871 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 872 ## #print "answer",answer 873 ## 874 ## assert allclose(quantity.vertex_values.ravel(), answer) 875 ## 876 ## 877 ## #Check that values can be set from file using default attribute 878 ## quantity.set_values(filename=txt_file, alpha=0) 879 ## assert allclose(quantity.vertex_values.ravel(), answer) 880 ## 881 ## #Cleanup 882 ## import os 883 ## os.remove(txt_file) 884 ## 885 ## def test_set_values_from_lat_long(self): 886 ## quantity = Quantity(self.mesh_onslow) 887 ## 888 ## #Get (enough) datapoints 889 ## data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 890 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]] 891 ## 892 ## data_geo_spatial = Geospatial_data(data_points, 893 ## points_are_lats_longs=True) 894 ## points_UTM = data_geo_spatial.get_data_points(absolute=True) 895 ## attributes = linear_function(points_UTM) 896 ## att = 'elevation' 897 ## 898 ## #Create .txt file 899 ## txt_file = tempfile.mktemp(".txt") 900 ## file = open(txt_file,"w") 901 ## file.write(" lat,long," + att + " \n") 902 ## for data_point, attribute in map(None, data_points, attributes): 903 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 904 ## + ',' + str(attribute) 905 ## #print "row", row 906 ## file.write(row + "\n") 907 ## file.close() 908 ## 909 ## 910 ## #Check that values can be set from file 911 ## quantity.set_values(filename=txt_file, 912 ## attribute_name=att, alpha=0) 913 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 914 ## 915 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 916 ## #print "answer",answer 917 ## 918 ## assert allclose(quantity.vertex_values.ravel(), answer) 919 ## 920 ## 921 ## #Check that values can be set from file using default attribute 922 ## quantity.set_values(filename=txt_file, alpha=0) 923 ## assert allclose(quantity.vertex_values.ravel(), answer) 924 ## 925 ## #Cleanup 926 ## import os 927 ## os.remove(txt_file) 928 ## 929 ## def test_set_values_from_UTM_pts(self): 930 ## quantity = Quantity(self.mesh_onslow) 931 ## 932 ## #Get (enough) datapoints 933 ## data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 934 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]] 935 ## 936 ## data_geo_spatial = Geospatial_data(data_points, 937 ## points_are_lats_longs=True) 938 ## points_UTM = data_geo_spatial.get_data_points(absolute=True) 939 ## attributes = linear_function(points_UTM) 940 ## att = 'elevation' 941 ## 942 ## #Create .txt file 943 ## txt_file = tempfile.mktemp(".txt") 944 ## file = open(txt_file,"w") 945 ## file.write(" x,y," + att + " \n") 946 ## for data_point, attribute in map(None, points_UTM, attributes): 947 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 948 ## + ',' + str(attribute) 949 ## #print "row", row 950 ## file.write(row + "\n") 951 ## file.close() 952 ## 953 ## 954 ## pts_file = tempfile.mktemp(".pts") 955 ## convert = Geospatial_data(txt_file) 956 ## convert.export_points_file(pts_file) 957 ## 958 ## #Check that values can be set from file 959 ## quantity.set_values_from_file(pts_file, att, 0, 960 ## 'vertices', None) 961 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 962 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 963 ## #print "answer",answer 964 ## assert allclose(quantity.vertex_values.ravel(), answer) 965 ## 966 ## #Check that values can be set from file 967 ## quantity.set_values(filename=pts_file, 968 ## attribute_name=att, alpha=0) 969 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 970 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 971 ## #print "answer",answer 972 ## assert allclose(quantity.vertex_values.ravel(), answer) 973 ## 974 ## 975 ## #Check that values can be set from file using default attribute 976 ## quantity.set_values(filename=txt_file, alpha=0) 977 ## assert allclose(quantity.vertex_values.ravel(), answer) 978 ## 979 ## #Cleanup 980 ## import os 981 ## os.remove(txt_file) 982 ## os.remove(pts_file) 983 ## 984 ## def verbose_test_set_values_from_UTM_pts(self): 985 ## quantity = Quantity(self.mesh_onslow) 986 ## 987 ## #Get (enough) datapoints 988 ## data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 989 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 990 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 991 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 992 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 993 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 994 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 995 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 996 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 997 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 998 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 999 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1000 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1001 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1002 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1003 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1004 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1005 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1006 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1007 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1008 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1009 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1010 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1011 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1012 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1013 ## [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65], 1014 ## [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6], 1015 ## ] 1016 ## 1017 ## data_geo_spatial = Geospatial_data(data_points, 1018 ## points_are_lats_longs=True) 1019 ## points_UTM = data_geo_spatial.get_data_points(absolute=True) 1020 ## attributes = linear_function(points_UTM) 1021 ## att = 'elevation' 1022 ## 1023 ## #Create .txt file 1024 ## txt_file = tempfile.mktemp(".txt") 1025 ## file = open(txt_file,"w") 1026 ## file.write(" x,y," + att + " \n") 1027 ## for data_point, attribute in map(None, points_UTM, attributes): 1028 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 1029 ## + ',' + str(attribute) 1030 ## #print "row", row 1031 ## file.write(row + "\n") 1032 ## file.close() 1033 ## 1034 ## 1035 ## pts_file = tempfile.mktemp(".pts") 1036 ## convert = Geospatial_data(txt_file) 1037 ## convert.export_points_file(pts_file) 1038 ## 1039 ## #Check that values can be set from file 1040 ## quantity.set_values_from_file(pts_file, att, 0, 1041 ## 'vertices', None, verbose = True, 1042 ## max_read_lines=2) 1043 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 1044 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 1045 ## #print "answer",answer 1046 ## assert allclose(quantity.vertex_values.ravel(), answer) 1047 ## 1048 ## #Check that values can be set from file 1049 ## quantity.set_values(filename=pts_file, 1050 ## attribute_name=att, alpha=0) 1051 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 1052 ## #print "quantity.vertex_values.ravel()", quantity.vertex_values.ravel() 1053 ## #print "answer",answer 1054 ## assert allclose(quantity.vertex_values.ravel(), answer) 1055 ## 1056 ## 1057 ## #Check that values can be set from file using default attribute 1058 ## quantity.set_values(filename=txt_file, alpha=0) 1059 ## assert allclose(quantity.vertex_values.ravel(), answer) 1060 ## 1061 ## #Cleanup 1062 ## import os 1063 ## os.remove(txt_file) 1064 ## os.remove(pts_file) 1065 ## 1066 ## def test_set_values_from_file_with_georef1(self): 1067 ## 1068 ## #Mesh in zone 56 (absolute coords) 1069 ## 1070 ## x0 = 314036.58727982 1071 ## y0 = 6224951.2960092 1072 ## 1073 ## a = [x0+0.0, y0+0.0] 1074 ## b = [x0+0.0, y0+2.0] 1075 ## c = [x0+2.0, y0+0.0] 1076 ## d = [x0+0.0, y0+4.0] 1077 ## e = [x0+2.0, y0+2.0] 1078 ## f = [x0+4.0, y0+0.0] 1079 ## 1080 ## points = [a, b, c, d, e, f] 1081 ## 1082 ## #bac, bce, ecf, dbe 1083 ## elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 1084 ## 1085 ## #absolute going in .. 1086 ## mesh4 = Domain(points, elements, 1087 ## geo_reference = Geo_reference(56, 0, 0)) 1088 ## mesh4.check_integrity() 1089 ## quantity = Quantity(mesh4) 1090 ## 1091 ## #Get (enough) datapoints (relative to georef) 1092 ## data_points_rel = [[ 0.66666667, 0.66666667], 1093 ## [ 1.33333333, 1.33333333], 1094 ## [ 2.66666667, 0.66666667], 1095 ## [ 0.66666667, 2.66666667], 1096 ## [ 0.0, 1.0], 1097 ## [ 0.0, 3.0], 1098 ## [ 1.0, 0.0], 1099 ## [ 1.0, 1.0], 1100 ## [ 1.0, 2.0], 1101 ## [ 1.0, 3.0], 1102 ## [ 2.0, 1.0], 1103 ## [ 3.0, 0.0], 1104 ## [ 3.0, 1.0]] 1105 ## 1106 ## data_geo_spatial = Geospatial_data(data_points_rel, 1107 ## geo_reference = Geo_reference(56, x0, y0)) 1108 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 1109 ## attributes = linear_function(data_points_absolute) 1110 ## att = 'spam_and_eggs' 1111 ## 1112 ## #Create .txt file 1113 ## ptsfile = tempfile.mktemp(".txt") 1114 ## file = open(ptsfile,"w") 1115 ## file.write(" x,y," + att + " \n") 1116 ## for data_point, attribute in map(None, data_points_absolute 1117 ## ,attributes): 1118 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 1119 ## + ',' + str(attribute) 1120 ## file.write(row + "\n") 1121 ## file.close() 1122 ## 1123 ## #file = open(ptsfile, 'r') 1124 ## #lines = file.readlines() 1125 ## #file.close() 1126 ## 1127 ## 1128 ## #Check that values can be set from file 1129 ## quantity.set_values(filename=ptsfile, 1130 ## attribute_name=att, alpha=0) 1131 ## answer = linear_function(quantity.domain.get_vertex_coordinates()) 1132 ## 1133 ## assert allclose(quantity.vertex_values.ravel(), answer) 1134 ## 1135 ## 1136 ## #Check that values can be set from file using default attribute 1137 ## quantity.set_values(filename=ptsfile, alpha=0) 1138 ## assert allclose(quantity.vertex_values.ravel(), answer) 1139 ## 1140 ## #Cleanup 1141 ## import os 1142 ## os.remove(ptsfile) 1143 ## 1144 ## 1145 ## def test_set_values_from_file_with_georef2(self): 1146 ## 1147 ## #Mesh in zone 56 (relative coords) 1148 ## 1149 ## x0 = 314036.58727982 1150 ## y0 = 6224951.2960092 1151 ## #x0 = 0.0 1152 ## #y0 = 0.0 1153 ## 1154 ## a = [0.0, 0.0] 1155 ## b = [0.0, 2.0] 1156 ## c = [2.0, 0.0] 1157 ## d = [0.0, 4.0] 1158 ## e = [2.0, 2.0] 1159 ## f = [4.0, 0.0] 1160 ## 1161 ## points = [a, b, c, d, e, f] 1162 ## 1163 ## #bac, bce, ecf, dbe 1164 ## elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 1165 ## 1166 ## mesh4 = Domain(points, elements, 1167 ## geo_reference = Geo_reference(56, x0, y0)) 1168 ## mesh4.check_integrity() 1169 ## quantity = Quantity(mesh4) 1170 ## 1171 ## #Get (enough) datapoints 1172 ## data_points = [[ x0+0.66666667, y0+0.66666667], 1173 ## [ x0+1.33333333, y0+1.33333333], 1174 ## [ x0+2.66666667, y0+0.66666667], 1175 ## [ x0+0.66666667, y0+2.66666667], 1176 ## [ x0+0.0, y0+1.0], 1177 ## [ x0+0.0, y0+3.0], 1178 ## [ x0+1.0, y0+0.0], 1179 ## [ x0+1.0, y0+1.0], 1180 ## [ x0+1.0, y0+2.0], 1181 ## [ x0+1.0, y0+3.0], 1182 ## [ x0+2.0, y0+1.0], 1183 ## [ x0+3.0, y0+0.0], 1184 ## [ x0+3.0, y0+1.0]] 1185 ## 1186 ## 1187 ## data_geo_spatial = Geospatial_data(data_points, 1188 ## geo_reference = Geo_reference(56, 0, 0)) 1189 ## data_points_absolute = data_geo_spatial.get_data_points(absolute=True) 1190 ## attributes = linear_function(data_points_absolute) 1191 ## att = 'spam_and_eggs' 1192 ## 1193 ## #Create .txt file 1194 ## ptsfile = tempfile.mktemp(".txt") 1195 ## file = open(ptsfile,"w") 1196 ## file.write(" x,y," + att + " \n") 1197 ## for data_point, attribute in map(None, data_points_absolute 1198 ## ,attributes): 1199 ## row = str(data_point[0]) + ',' + str(data_point[1]) \ 1200 ## + ',' + str(attribute) 1201 ## file.write(row + "\n") 1202 ## file.close() 1203 ## 1204 ## 1205 ## #Check that values can be set from file 1206 ## quantity.set_values(filename=ptsfile, 1207 ## attribute_name=att, alpha=0) 1208 ## answer = linear_function(quantity.domain. \ 1209 ## get_vertex_coordinates(absolute=True)) 1210 ## 1211 ## 1212 ## assert allclose(quantity.vertex_values.ravel(), answer) 1213 ## 1214 ## 1215 ## #Check that values can be set from file using default attribute 1216 ## quantity.set_values(filename=ptsfile, alpha=0) 1217 ## assert allclose(quantity.vertex_values.ravel(), answer) 1218 ## 1219 ## #Cleanup 1220 ## import os 1221 ## os.remove(ptsfile) 1222 ## 1223 ## 1224 ## 1225 ## 1226 ## def test_set_values_from_quantity(self): 1227 ## 1228 ## quantity1 = Quantity(self.mesh4) 1229 ## quantity1.set_vertex_values([0,1,2,3,4,5]) 1230 ## 1231 ## assert allclose(quantity1.vertex_values, 1232 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1233 ## 1234 ## 1235 ## quantity2 = Quantity(self.mesh4) 1236 ## quantity2.set_values(quantity=quantity1) 1237 ## assert allclose(quantity2.vertex_values, 1238 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1239 ## 1240 ## quantity2.set_values(quantity = 2*quantity1) 1241 ## assert allclose(quantity2.vertex_values, 1242 ## [[2,0,4], [2,4,8], [8,4,10], [6,2,8]]) 1243 ## 1244 ## quantity2.set_values(quantity = 2*quantity1 + 3) 1245 ## assert allclose(quantity2.vertex_values, 1246 ## [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 1247 ## 1248 ## 1249 ## #Check detection of quantity as first orgument 1250 ## quantity2.set_values(2*quantity1 + 3) 1251 ## assert allclose(quantity2.vertex_values, 1252 ## [[5,3,7], [5,7,11], [11,7,13], [9,5,11]]) 1253 ## 1254 ## 1255 ## 1256 ## def Xtest_set_values_from_quantity_using_polygon(self): 1257 ## """test_set_values_from_quantity_using_polygon(self): 1258 ## 1259 ## Check that polygon can be used to restrict set_values when 1260 ## using another quantity as argument. 1261 ## """ 1262 ## 1263 ## # Create restricting polygon (containing node #4 (2,2) and 1264 ## # centroid of triangle #1 (bce) 1265 ## polygon = [[1.0, 1.0], [4.0, 1.0], 1266 ## [4.0, 4.0], [1.0, 4.0]] 1267 ## assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 1268 ## 1269 ## quantity1 = Quantity(self.mesh4) 1270 ## quantity1.set_vertex_values([0,1,2,3,4,5]) 1271 ## 1272 ## assert allclose(quantity1.vertex_values, 1273 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1274 ## 1275 ## 1276 ## quantity2 = Quantity(self.mesh4) 1277 ## quantity2.set_values(quantity=quantity1, 1278 ## polygon=polygon) 1279 ## 1280 ## msg = 'Only node #4(e) at (2,2) should have values applied ' 1281 ## assert allclose(quantity2.vertex_values, 1282 ## [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg 1283 ## #bac, bce, ecf, dbe 1284 ## 1285 ## 1286 ## 1287 ## def test_overloading(self): 1288 ## 1289 ## quantity1 = Quantity(self.mesh4) 1290 ## quantity1.set_vertex_values([0,1,2,3,4,5]) 1291 ## 1292 ## assert allclose(quantity1.vertex_values, 1293 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1294 ## 1295 ## 1296 ## quantity2 = Quantity(self.mesh4) 1297 ## quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], 1298 ## location = 'vertices') 1299 ## 1300 ## 1301 ## 1302 ## quantity3 = Quantity(self.mesh4) 1303 ## quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]], 1304 ## location = 'vertices') 1305 ## 1306 ## 1307 ## # Negation 1308 ## Q = -quantity1 1309 ## assert allclose(Q.vertex_values, -quantity1.vertex_values) 1310 ## assert allclose(Q.centroid_values, -quantity1.centroid_values) 1311 ## assert allclose(Q.edge_values, -quantity1.edge_values) 1312 ## 1313 ## # Addition 1314 ## Q = quantity1 + 7 1315 ## assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 1316 ## assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 1317 ## assert allclose(Q.edge_values, quantity1.edge_values + 7) 1318 ## 1319 ## Q = 7 + quantity1 1320 ## assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 1321 ## assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 1322 ## assert allclose(Q.edge_values, quantity1.edge_values + 7) 1323 ## 1324 ## Q = quantity1 + quantity2 1325 ## assert allclose(Q.vertex_values, 1326 ## quantity1.vertex_values + quantity2.vertex_values) 1327 ## assert allclose(Q.centroid_values, 1328 ## quantity1.centroid_values + quantity2.centroid_values) 1329 ## assert allclose(Q.edge_values, 1330 ## quantity1.edge_values + quantity2.edge_values) 1331 ## 1332 ## 1333 ## Q = quantity1 + quantity2 - 3 1334 ## assert allclose(Q.vertex_values, 1335 ## quantity1.vertex_values + quantity2.vertex_values - 3) 1336 ## 1337 ## Q = quantity1 - quantity2 1338 ## assert allclose(Q.vertex_values, 1339 ## quantity1.vertex_values - quantity2.vertex_values) 1340 ## 1341 ## #Scaling 1342 ## Q = quantity1*3 1343 ## assert allclose(Q.vertex_values, quantity1.vertex_values*3) 1344 ## assert allclose(Q.centroid_values, quantity1.centroid_values*3) 1345 ## assert allclose(Q.edge_values, quantity1.edge_values*3) 1346 ## Q = 3*quantity1 1347 ## assert allclose(Q.vertex_values, quantity1.vertex_values*3) 1348 ## 1349 ## #Multiplication 1350 ## Q = quantity1 * quantity2 1351 ## #print Q.vertex_values 1352 ## #print Q.centroid_values 1353 ## #print quantity1.centroid_values 1354 ## #print quantity2.centroid_values 1355 ## 1356 ## assert allclose(Q.vertex_values, 1357 ## quantity1.vertex_values * quantity2.vertex_values) 1358 ## 1359 ## #Linear combinations 1360 ## Q = 4*quantity1 + 2 1361 ## assert allclose(Q.vertex_values, 1362 ## 4*quantity1.vertex_values + 2) 1363 ## 1364 ## Q = quantity1*quantity2 + 2 1365 ## assert allclose(Q.vertex_values, 1366 ## quantity1.vertex_values * quantity2.vertex_values + 2) 1367 ## 1368 ## Q = quantity1*quantity2 + quantity3 1369 ## assert allclose(Q.vertex_values, 1370 ## quantity1.vertex_values * quantity2.vertex_values + 1371 ## quantity3.vertex_values) 1372 ## Q = quantity1*quantity2 + 3*quantity3 1373 ## assert allclose(Q.vertex_values, 1374 ## quantity1.vertex_values * quantity2.vertex_values + 1375 ## 3*quantity3.vertex_values) 1376 ## Q = quantity1*quantity2 + 3*quantity3 + 5.0 1377 ## assert allclose(Q.vertex_values, 1378 ## quantity1.vertex_values * quantity2.vertex_values + 1379 ## 3*quantity3.vertex_values + 5) 1380 ## 1381 ## Q = quantity1*quantity2 - quantity3 1382 ## assert allclose(Q.vertex_values, 1383 ## quantity1.vertex_values * quantity2.vertex_values - 1384 ## quantity3.vertex_values) 1385 ## Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0 1386 ## assert allclose(Q.vertex_values, 1387 ## 1.5*quantity1.vertex_values * quantity2.vertex_values - 1388 ## 3*quantity3.vertex_values + 5) 1389 ## 1390 ## #Try combining quantities and arrays and scalars 1391 ## Q = 1.5*quantity1*quantity2.vertex_values -\ 1392 ## 3*quantity3.vertex_values + 5.0 1393 ## assert allclose(Q.vertex_values, 1394 ## 1.5*quantity1.vertex_values * quantity2.vertex_values - 1395 ## 3*quantity3.vertex_values + 5) 1396 ## 1397 ## 1398 ## #Powers 1399 ## Q = quantity1**2 1400 ## assert allclose(Q.vertex_values, quantity1.vertex_values**2) 1401 ## 1402 ## Q = quantity1**2 +quantity2**2 1403 ## assert allclose(Q.vertex_values, 1404 ## quantity1.vertex_values**2 + \ 1405 ## quantity2.vertex_values**2) 1406 ## 1407 ## Q = (quantity1**2 +quantity2**2)**0.5 1408 ## assert allclose(Q.vertex_values, 1409 ## (quantity1.vertex_values**2 + \ 1410 ## quantity2.vertex_values**2)**0.5) 1411 ## 1412 ## 1413 ## 1414 ## 1415 ## 1416 ## 1417 ## 1418 ## def test_compute_gradient(self): 1419 ## quantity = Quantity(self.mesh4) 1420 ## 1421 ## #Set up for a gradient of (2,0) at mid triangle 1422 ## quantity.set_values([2.0, 4.0, 6.0, 2.0], 1423 ## location = 'centroids') 1424 ## 1425 ## 1426 ## #Gradients 1427 ## quantity.compute_gradients() 1428 ## 1429 ## a = quantity.x_gradient 1430 ## b = quantity.y_gradient 1431 ## #print self.mesh4.centroid_coordinates 1432 ## #print a, b 1433 ## 1434 ## #The central triangle (1) 1435 ## #(using standard gradient based on neigbours controid values) 1436 ## assert allclose(a[1], 2.0) 1437 ## assert allclose(b[1], 0.0) 1438 ## 1439 ## 1440 ## #Left triangle (0) using two point gradient 1441 ## #q0 = q1 + a*(x0-x1) + b*(y0-y1) <=> 1442 ## #2 = 4 + a*(-2/3) + b*(-2/3) 1443 ## assert allclose(a[0] + b[0], 3) 1444 ## #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 1445 ## assert allclose(a[0] - b[0], 0) 1446 ## 1447 ## 1448 ## #Right triangle (2) using two point gradient 1449 ## #q2 = q1 + a*(x2-x1) + b*(y2-y1) <=> 1450 ## #6 = 4 + a*(4/3) + b*(-2/3) 1451 ## assert allclose(2*a[2] - b[2], 3) 1452 ## #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0) 1453 ## assert allclose(a[2] + 2*b[2], 0) 1454 ## 1455 ## 1456 ## #Top triangle (3) using two point gradient 1457 ## #q3 = q1 + a*(x3-x1) + b*(y3-y1) <=> 1458 ## #2 = 4 + a*(-2/3) + b*(4/3) 1459 ## assert allclose(a[3] - 2*b[3], 3) 1460 ## #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 1461 ## assert allclose(2*a[3] + b[3], 0) 1462 ## 1463 ## 1464 ## 1465 ## #print a, b 1466 ## quantity.extrapolate_second_order() 1467 ## 1468 ## #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc) 1469 ## assert allclose(quantity.vertex_values[0,:], [3., 0., 3.]) 1470 ## assert allclose(quantity.vertex_values[1,:], [4./3, 16./3, 16./3]) 1471 ## 1472 ## 1473 ## #a = 1.2, b=-0.6 1474 ## #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3) 1475 ## assert allclose(quantity.vertex_values[2,2], 8) 1476 ## 1477 ## def test_get_gradients(self): 1478 ## quantity = Quantity(self.mesh4) 1479 ## 1480 ## #Set up for a gradient of (2,0) at mid triangle 1481 ## quantity.set_values([2.0, 4.0, 6.0, 2.0], 1482 ## location = 'centroids') 1483 ## 1484 ## 1485 ## #Gradients 1486 ## quantity.compute_gradients() 1487 ## 1488 ## a, b = quantity.get_gradients() 1489 ## #print self.mesh4.centroid_coordinates 1490 ## #print a, b 1491 ## 1492 ## #The central triangle (1) 1493 ## #(using standard gradient based on neigbours controid values) 1494 ## assert allclose(a[1], 2.0) 1495 ## assert allclose(b[1], 0.0) 1496 ## 1497 ## 1498 ## #Left triangle (0) using two point gradient 1499 ## #q0 = q1 + a*(x0-x1) + b*(y0-y1) <=> 1500 ## #2 = 4 + a*(-2/3) + b*(-2/3) 1501 ## assert allclose(a[0] + b[0], 3) 1502 ## #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 1503 ## assert allclose(a[0] - b[0], 0) 1504 ## 1505 ## 1506 ## #Right triangle (2) using two point gradient 1507 ## #q2 = q1 + a*(x2-x1) + b*(y2-y1) <=> 1508 ## #6 = 4 + a*(4/3) + b*(-2/3) 1509 ## assert allclose(2*a[2] - b[2], 3) 1510 ## #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0) 1511 ## assert allclose(a[2] + 2*b[2], 0) 1512 ## 1513 ## 1514 ## #Top triangle (3) using two point gradient 1515 ## #q3 = q1 + a*(x3-x1) + b*(y3-y1) <=> 1516 ## #2 = 4 + a*(-2/3) + b*(4/3) 1517 ## assert allclose(a[3] - 2*b[3], 3) 1518 ## #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 1519 ## assert allclose(2*a[3] + b[3], 0) 1520 ## 1521 ## 1522 ## def test_second_order_extrapolation2(self): 1523 ## quantity = Quantity(self.mesh4) 1524 ## 1525 ## #Set up for a gradient of (3,1), f(x) = 3x+y 1526 ## quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3], 1527 ## location = 'centroids') 1528 ## 1529 ## #Gradients 1530 ## quantity.compute_gradients() 1531 ## 1532 ## a = quantity.x_gradient 1533 ## b = quantity.y_gradient 1534 ## 1535 ## #print a, b 1536 ## 1537 ## assert allclose(a[1], 3.0) 1538 ## assert allclose(b[1], 1.0) 1539 ## 1540 ## #Work out the others 1541 ## 1542 ## quantity.extrapolate_second_order() 1543 ## 1544 ## #print quantity.vertex_values 1545 ## assert allclose(quantity.vertex_values[1,0], 2.0) 1546 ## assert allclose(quantity.vertex_values[1,1], 6.0) 1547 ## assert allclose(quantity.vertex_values[1,2], 8.0) 1548 ## 1549 ## 1550 ## 1551 ## def test_backup_saxpy_centroid_values(self): 1552 ## quantity = Quantity(self.mesh4) 1553 ## 1554 ## #Set up for a gradient of (3,1), f(x) = 3x+y 1555 ## c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3]) 1556 ## d_values = array([1.0, 2.0, 3.0, 4.0]) 1557 ## quantity.set_values(c_values, location = 'centroids') 1558 ## 1559 ## #Backup 1560 ## quantity.backup_centroid_values() 1561 ## 1562 ## #print quantity.vertex_values 1563 ## assert allclose(quantity.centroid_values, quantity.centroid_backup_values) 1564 ## 1565 ## 1566 ## quantity.set_values(d_values, location = 'centroids') 1567 ## 1568 ## quantity.saxpy_centroid_values(2.0, 3.0) 1569 ## 1570 ## assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values) 1571 ## 1572 ## 1573 ## 1574 ## def test_first_order_extrapolator(self): 1575 ## quantity = Quantity(self.mesh4) 1576 ## 1577 ## #Test centroids 1578 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1579 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1580 ## 1581 ## #Extrapolate 1582 ## quantity.extrapolate_first_order() 1583 ## 1584 ## #Check that gradient is zero 1585 ## a,b = quantity.get_gradients() 1586 ## assert allclose(a, [0,0,0,0]) 1587 ## assert allclose(b, [0,0,0,0]) 1588 ## 1589 ## #Check vertices but not edge values 1590 ## assert allclose(quantity.vertex_values, 1591 ## [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 1592 ## 1593 ## 1594 ## def test_second_order_extrapolator(self): 1595 ## quantity = Quantity(self.mesh4) 1596 ## 1597 ## #Set up for a gradient of (3,0) at mid triangle 1598 ## quantity.set_values([2.0, 4.0, 8.0, 2.0], 1599 ## location = 'centroids') 1600 ## 1601 ## 1602 ## 1603 ## quantity.extrapolate_second_order() 1604 ## quantity.limit() 1605 ## 1606 ## 1607 ## #Assert that central triangle is limited by neighbours 1608 ## assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 1609 ## assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1] 1610 ## 1611 ## assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 1612 ## assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 1613 ## 1614 ## assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 1615 ## assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1] 1616 ## 1617 ## 1618 ## #Assert that quantities are conserved 1619 ### from numpy.oldnumeric import sum 1620 ## from numpy import sum 1621 ## for k in range(quantity.centroid_values.shape[0]): 1622 ## assert allclose (quantity.centroid_values[k], 1623 ## sum(quantity.vertex_values[k,:])/3) 1624 ## 1625 ## 1626 ## 1627 ## 1628 ## 1629 ## def test_limit_vertices_by_all_neighbours(self): 1630 ## quantity = Quantity(self.mesh4) 1631 ## 1632 ## #Create a deliberate overshoot (e.g. from gradient computation) 1633 ## quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1634 ## 1635 ## 1636 ## #Limit 1637 ## quantity.limit_vertices_by_all_neighbours() 1638 ## 1639 ## #Assert that central triangle is limited by neighbours 1640 ## assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 1641 ## assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1] 1642 ## 1643 ## assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 1644 ## assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 1645 ## 1646 ## assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 1647 ## assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1] 1648 ## 1649 ## 1650 ## 1651 ## #Assert that quantities are conserved 1652 ### from numpy.oldnumeric import sum 1653 ## from numpy import sum 1654 ## for k in range(quantity.centroid_values.shape[0]): 1655 ## assert allclose (quantity.centroid_values[k], 1656 ## sum(quantity.vertex_values[k,:])/3) 1657 ## 1658 ## 1659 ## 1660 ## def test_limit_edges_by_all_neighbours(self): 1661 ## quantity = Quantity(self.mesh4) 1662 ## 1663 ## #Create a deliberate overshoot (e.g. from gradient computation) 1664 ## quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1665 ## 1666 ## 1667 ## #Limit 1668 ## quantity.limit_edges_by_all_neighbours() 1669 ## 1670 ## #Assert that central triangle is limited by neighbours 1671 ## assert quantity.edge_values[1,0] <= quantity.centroid_values[2] 1672 ## assert quantity.edge_values[1,0] >= quantity.centroid_values[0] 1673 ## 1674 ## assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1675 ## assert quantity.edge_values[1,1] >= quantity.centroid_values[0] 1676 ## 1677 ## assert quantity.edge_values[1,2] <= quantity.centroid_values[2] 1678 ## assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1679 ## 1680 ## 1681 ## 1682 ## #Assert that quantities are conserved 1683 ### from numpy.oldnumeric import sum 1684 ## from numpy import sum 1685 ## for k in range(quantity.centroid_values.shape[0]): 1686 ## assert allclose (quantity.centroid_values[k], 1687 ## sum(quantity.vertex_values[k,:])/3) 1688 ## 1689 ## 1690 ## def test_limit_edges_by_neighbour(self): 1691 ## quantity = Quantity(self.mesh4) 1692 ## 1693 ## #Create a deliberate overshoot (e.g. from gradient computation) 1694 ## quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1695 ## 1696 ## 1697 ## #Limit 1698 ## quantity.limit_edges_by_neighbour() 1699 ## 1700 ## #Assert that central triangle is limited by neighbours 1701 ## assert quantity.edge_values[1,0] <= quantity.centroid_values[3] 1702 ## assert quantity.edge_values[1,0] >= quantity.centroid_values[1] 1703 ## 1704 ## assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1705 ## assert quantity.edge_values[1,1] >= quantity.centroid_values[1] 1706 ## 1707 ## assert quantity.edge_values[1,2] <= quantity.centroid_values[1] 1708 ## assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1709 ## 1710 ## 1711 ## 1712 ## #Assert that quantities are conserved 1713 ### from numpy.oldnumeric import sum 1714 ## from numpy import sum 1715 ## for k in range(quantity.centroid_values.shape[0]): 1716 ## assert allclose (quantity.centroid_values[k], 1717 ## sum(quantity.vertex_values[k,:])/3) 1718 ## 1719 ## def test_limiter2(self): 1720 ## """Taken from test_shallow_water 1721 ## """ 1722 ## quantity = Quantity(self.mesh4) 1723 ## quantity.domain.beta_w = 0.9 1724 ## 1725 ## #Test centroids 1726 ## quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1727 ## assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1728 ## 1729 ## 1730 ## #Extrapolate 1731 ## quantity.extrapolate_second_order() 1732 ## 1733 ## assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1734 ## 1735 ## #Limit 1736 ## quantity.limit() 1737 ## 1738 ## # limited value for beta_w = 0.9 1739 ## 1740 ## assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9]) 1741 ## # limited values for beta_w = 0.5 1742 ## #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5]) 1743 ## 1744 ## 1745 ## #Assert that quantities are conserved 1746 ### from numpy.oldnumeric import sum 1747 ## from numpy import sum 1748 ## for k in range(quantity.centroid_values.shape[0]): 1749 ## assert allclose (quantity.centroid_values[k], 1750 ## sum(quantity.vertex_values[k,:])/3) 1751 ## 1752 ## 1753 ## 1754 ## 1755 ## 1756 ## def test_distribute_first_order(self): 1757 ## quantity = Quantity(self.mesh4) 1758 ## 1759 ## #Test centroids 1760 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1761 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1762 ## 1763 ## 1764 ## #Extrapolate from centroid to vertices and edges 1765 ## quantity.extrapolate_first_order() 1766 ## 1767 ## #Interpolate 1768 ## #quantity.interpolate_from_vertices_to_edges() 1769 ## 1770 ## assert allclose(quantity.vertex_values, 1771 ## [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 1772 ## assert allclose(quantity.edge_values, [[1,1,1], [2,2,2], 1773 ## [3,3,3], [4, 4, 4]]) 1774 ## 1775 ## 1776 ## def test_interpolate_from_vertices_to_edges(self): 1777 ## quantity = Quantity(self.mesh4) 1778 ## 1779 ## quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],float) 1780 ## 1781 ## quantity.interpolate_from_vertices_to_edges() 1782 ## 1783 ## assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 1784 ## [3., 2.5, 1.5], 1785 ## [3.5, 4.5, 3.], 1786 ## [2.5, 3.5, 2]]) 1787 ## 1788 ## 1789 ## def test_interpolate_from_edges_to_vertices(self): 1790 ## quantity = Quantity(self.mesh4) 1791 ## 1792 ## quantity.edge_values = array([[1., 1.5, 0.5], 1793 ## [3., 2.5, 1.5], 1794 ## [3.5, 4.5, 3.], 1795 ## [2.5, 3.5, 2]],float) 1796 ## 1797 ## quantity.interpolate_from_edges_to_vertices() 1798 ## 1799 ## assert allclose(quantity.vertex_values, 1800 ## [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1801 ## 1802 ## 1803 ## 1804 ## def test_distribute_second_order(self): 1805 ## quantity = Quantity(self.mesh4) 1806 ## 1807 ## #Test centroids 1808 ## quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1809 ## assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1810 ## 1811 ## 1812 ## #Extrapolate 1813 ## quantity.extrapolate_second_order() 1814 ## 1815 ## assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1816 ## 1817 ## 1818 ## def test_update_explicit(self): 1819 ## quantity = Quantity(self.mesh4) 1820 ## 1821 ## #Test centroids 1822 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1823 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1824 ## 1825 ## #Set explicit_update 1826 ## quantity.explicit_update = array( [1.,1.,1.,1.] ) 1827 ## 1828 ## #Update with given timestep 1829 ## quantity.update(0.1) 1830 ## 1831 ## x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] ) 1832 ## assert allclose( quantity.centroid_values, x) 1833 ## 1834 ## def test_update_semi_implicit(self): 1835 ## quantity = Quantity(self.mesh4) 1836 ## 1837 ## #Test centroids 1838 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1839 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1840 ## 1841 ## #Set semi implicit update 1842 ## quantity.semi_implicit_update = array([1.,1.,1.,1.]) 1843 ## 1844 ## #Update with given timestep 1845 ## timestep = 0.1 1846 ## quantity.update(timestep) 1847 ## 1848 ## sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 1849 ## denom = ones(4, float)-timestep*sem 1850 ## 1851 ## x = array([1, 2, 3, 4])/denom 1852 ## assert allclose( quantity.centroid_values, x) 1853 ## 1854 ## 1855 ## def test_both_updates(self): 1856 ## quantity = Quantity(self.mesh4) 1857 ## 1858 ## #Test centroids 1859 ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1860 ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1861 ## 1862 ## #Set explicit_update 1863 ## quantity.explicit_update = array( [4.,3.,2.,1.] ) 1864 ## 1865 ## #Set semi implicit update 1866 ## quantity.semi_implicit_update = array( [1.,1.,1.,1.] ) 1867 ## 1868 ## #Update with given timestep 1869 ## timestep = 0.1 1870 ## quantity.update(0.1) 1871 ## 1872 ## sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 1873 ## denom = ones(4, float)-timestep*sem 1874 ## 1875 ## x = array([1., 2., 3., 4.]) 1876 ## x /= denom 1877 ## x += timestep*array( [4.0, 3.0, 2.0, 1.0] ) 1878 ## 1879 ## assert allclose( quantity.centroid_values, x) 1880 ## 1881 ## 1882 ## 1883 ## 1884 ## #Test smoothing 1885 ## def test_smoothing(self): 1886 ## 1887 ## from mesh_factory import rectangular 1888 ## from shallow_water import Domain, Transmissive_boundary 1889 ### from numpy.oldnumeric import zeros, Float 1890 ## from numpy import zeros, float 1891 ## from anuga.utilities.numerical_tools import mean 1892 ## 1893 ## #Create basic mesh 1894 ## points, vertices, boundary = rectangular(2, 2) 1895 ## 1896 ## #Create shallow water domain 1897 ## domain = Domain(points, vertices, boundary) 1898 ## domain.default_order=2 1899 ## domain.reduction = mean 1900 ## 1901 ## 1902 ## #Set some field values 1903 ## domain.set_quantity('elevation', lambda x,y: x) 1904 ## domain.set_quantity('friction', 0.03) 1905 ## 1906 ## 1907 ## ###################### 1908 ## # Boundary conditions 1909 ## B = Transmissive_boundary(domain) 1910 ## domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B}) 1911 ## 1912 ## 1913 ## ###################### 1914 ## #Initial condition - with jumps 1915 ## 1916 ## bed = domain.quantities['elevation'].vertex_values 1917 ## stage = zeros(bed.shape, float) 1918 ## 1919 ## h = 0.03 1920 ## for i in range(stage.shape[0]): 1921 ## if i % 2 == 0: 1922 ## stage[i,:] = bed[i,:] + h 1923 ## else: 1924 ## stage[i,:] = bed[i,:] 1925 ## 1926 ## domain.set_quantity('stage', stage) 1927 ## 1928 ## stage = domain.quantities['stage'] 1929 ## 1930 ## #Get smoothed stage 1931 ## A, V = stage.get_vertex_values(xy=False, smooth=True) 1932 ## Q = stage.vertex_values 1933 ## 1934 ## 1935 ## assert A.shape[0] == 9 1936 ## assert V.shape[0] == 8 1937 ## assert V.shape[1] == 3 1938 ## 1939 ## #First four points 1940 ## assert allclose(A[0], (Q[0,2] + Q[1,1])/2) 1941 ## assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3) 1942 ## assert allclose(A[2], Q[3,0]) 1943 ## assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3) 1944 ## 1945 ## #Center point 1946 ## assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\ 1947 ## Q[5,0] + Q[6,2] + Q[7,1])/6) 1948 ## 1949 ## 1950 ## #Check V 1951 ## assert allclose(V[0,:], [3,4,0]) 1952 ## assert allclose(V[1,:], [1,0,4]) 1953 ## assert allclose(V[2,:], [4,5,1]) 1954 ## assert allclose(V[3,:], [2,1,5]) 1955 ## assert allclose(V[4,:], [6,7,3]) 1956 ## assert allclose(V[5,:], [4,3,7]) 1957 ## assert allclose(V[6,:], [7,8,4]) 1958 ## assert allclose(V[7,:], [5,4,8]) 1959 ## 1960 ## #Get smoothed stage with XY 1961 ## X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True) 1962 ## 1963 ## assert allclose(A, A1) 1964 ## assert allclose(V, V1) 1965 ## 1966 ## #Check XY 1967 ## assert allclose(X[4], 0.5) 1968 ## assert allclose(Y[4], 0.5) 1969 ## 1970 ## assert allclose(X[7], 1.0) 1971 ## assert allclose(Y[7], 0.5) 1972 ## 1973 ## 1974 ## 1975 ## 1976 ## def test_vertex_values_no_smoothing(self): 1977 ## 1978 ## from mesh_factory import rectangular 1979 ## from shallow_water import Domain, Transmissive_boundary 1980 ### from numpy.oldnumeric import zeros, Float 1981 ## from numpy import zeros, float 1982 ## from anuga.utilities.numerical_tools import mean 1983 ## 1984 ## 1985 ## #Create basic mesh 1986 ## points, vertices, boundary = rectangular(2, 2) 1987 ## 1988 ## #Create shallow water domain 1989 ## domain = Domain(points, vertices, boundary) 1990 ## domain.default_order=2 1991 ## domain.reduction = mean 1992 ## 1993 ## 1994 ## #Set some field values 1995 ## domain.set_quantity('elevation', lambda x,y: x) 1996 ## domain.set_quantity('friction', 0.03) 1997 ## 1998 ## 1999 ## ###################### 2000 ## #Initial condition - with jumps 2001 ## 2002 ## bed = domain.quantities['elevation'].vertex_values 2003 ## stage = zeros(bed.shape, float) 2004 ## 2005 ## h = 0.03 2006 ## for i in range(stage.shape[0]): 2007 ## if i % 2 == 0: 2008 ## stage[i,:] = bed[i,:] + h 2009 ## else: 2010 ## stage[i,:] = bed[i,:] 2011 ## 2012 ## domain.set_quantity('stage', stage) 2013 ## 2014 ## #Get stage 2015 ## stage = domain.quantities['stage'] 2016 ## A, V = stage.get_vertex_values(xy=False, smooth=False) 2017 ## Q = stage.vertex_values.ravel() 2018 ## 2019 ## for k in range(8): 2020 ## assert allclose(A[k], Q[k]) 2021 ## 2022 ## 2023 ## for k in range(8): 2024 ## assert V[k, 0] == 3*k 2025 ## assert V[k, 1] == 3*k+1 2026 ## assert V[k, 2] == 3*k+2 2027 ## 2028 ## 2029 ## 2030 ## X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False) 2031 ## 2032 ## 2033 ## assert allclose(A, A1) 2034 ## assert allclose(V, V1) 2035 ## 2036 ## #Check XY 2037 ## assert allclose(X[1], 0.5) 2038 ## assert allclose(Y[1], 0.5) 2039 ## assert allclose(X[4], 0.0) 2040 ## assert allclose(Y[4], 0.0) 2041 ## assert allclose(X[12], 1.0) 2042 ## assert allclose(Y[12], 0.0) 2043 ## 2044 ## 2045 ## 2046 ## def set_array_values_by_index(self): 2047 ## 2048 ## from mesh_factory import rectangular 2049 ## from shallow_water import Domain 2050 ### from numpy.oldnumeric import zeros, Float 2051 ## from numpy import zeros, float 2052 ## 2053 ## #Create basic mesh 2054 ## points, vertices, boundary = rectangular(1, 1) 2055 ## 2056 ## #Create shallow water domain 2057 ## domain = Domain(points, vertices, boundary) 2058 ## #print "domain.number_of_elements ",domain.number_of_elements 2059 ## quantity = Quantity(domain,[[1,1,1],[2,2,2]]) 2060 ## value = [7] 2061 ## indices = [1] 2062 ## quantity.set_array_values_by_index(value, 2063 ## location = 'centroids', 2064 ## indices = indices) 2065 ## #print "quantity.centroid_values",quantity.centroid_values 2066 ## 2067 ## assert allclose(quantity.centroid_values, [1,7]) 2068 ## 2069 ## quantity.set_array_values([15,20,25], indices = indices) 2070 ## assert allclose(quantity.centroid_values, [1,20]) 2071 ## 2072 ## quantity.set_array_values([15,20,25], indices = indices) 2073 ## assert allclose(quantity.centroid_values, [1,20]) 2074 ## 2075 ## def test_setting_some_vertex_values(self): 2076 ## """ 2077 ## set values based on triangle lists. 2078 ## """ 2079 ## from mesh_factory import rectangular 2080 ## from shallow_water import Domain 2081 ### from numpy.oldnumeric import zeros, Float 2082 ## from numpy import zeros, float 2083 ## 2084 ## #Create basic mesh 2085 ## points, vertices, boundary = rectangular(1, 3) 2086 ## #print "vertices",vertices 2087 ## #Create shallow water domain 2088 ## domain = Domain(points, vertices, boundary) 2089 ## #print "domain.number_of_elements ",domain.number_of_elements 2090 ## quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2091 ## [4,4,4],[5,5,5],[6,6,6]]) 2092 ## 2093 ## 2094 ## # Check that constants work 2095 ## value = 7 2096 ## indices = [1] 2097 ## quantity.set_values(value, 2098 ## location = 'centroids', 2099 ## indices = indices) 2100 ## #print "quantity.centroid_values",quantity.centroid_values 2101 ## assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2102 ## 2103 ## value = [7] 2104 ## indices = [1] 2105 ## quantity.set_values(value, 2106 ## location = 'centroids', 2107 ## indices = indices) 2108 ## #print "quantity.centroid_values",quantity.centroid_values 2109 ## assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2110 ## 2111 ## value = [[15,20,25]] 2112 ## quantity.set_values(value, indices = indices) 2113 ## #print "1 quantity.vertex_values",quantity.vertex_values 2114 ## assert allclose(quantity.vertex_values[1], value[0]) 2115 ## 2116 ## 2117 ## #print "quantity",quantity.vertex_values 2118 ## values = [10,100,50] 2119 ## quantity.set_values(values, indices = [0,1,5], location = 'centroids') 2120 ## #print "2 quantity.vertex_values",quantity.vertex_values 2121 ## assert allclose(quantity.vertex_values[0], [10,10,10]) 2122 ## assert allclose(quantity.vertex_values[5], [50,50,50]) 2123 ## #quantity.interpolate() 2124 ## #print "quantity.centroid_values",quantity.centroid_values 2125 ## assert allclose(quantity.centroid_values, [10,100,3,4,5,50]) 2126 ## 2127 ## 2128 ## quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2129 ## [4,4,4],[5,5,5],[6,6,6]]) 2130 ## values = [10,100,50] 2131 ## #this will be per unique vertex, indexing the vertices 2132 ## #print "quantity.vertex_values",quantity.vertex_values 2133 ## quantity.set_values(values, indices = [0,1,5]) 2134 ## #print "quantity.vertex_values",quantity.vertex_values 2135 ## assert allclose(quantity.vertex_values[0], [1,50,10]) 2136 ## assert allclose(quantity.vertex_values[5], [6,6,6]) 2137 ## assert allclose(quantity.vertex_values[1], [100,10,50]) 2138 ## 2139 ## quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2140 ## [4,4,4],[5,5,5],[6,6,6]]) 2141 ## values = [[31,30,29],[400,400,400],[1000,999,998]] 2142 ## quantity.set_values(values, indices = [3,3,5]) 2143 ## quantity.interpolate() 2144 ## assert allclose(quantity.centroid_values, [1,2,3,400,5,999]) 2145 ## 2146 ## values = [[1,1,1],[2,2,2],[3,3,3], 2147 ## [4,4,4],[5,5,5],[6,6,6]] 2148 ## quantity.set_values(values) 2149 ## 2150 ## # testing the standard set values by vertex 2151 ## # indexed by vertex_id in general_mesh.coordinates 2152 ## values = [0,1,2,3,4,5,6,7] 2153 ## 2154 ## quantity.set_values(values) 2155 ## #print "1 quantity.vertex_values",quantity.vertex_values 2156 ## assert allclose(quantity.vertex_values,[[ 4., 5., 0.], 2157 ## [ 1., 0., 5.], 2158 ## [ 5., 6., 1.], 2159 ## [ 2., 1., 6.], 2160 ## [ 6., 7., 2.], 2161 ## [ 3., 2., 7.]]) 2162 ## 2163 ## def test_setting_unique_vertex_values(self): 2164 ## """ 2165 ## set values based on unique_vertex lists. 2166 ## """ 2167 ## from mesh_factory import rectangular 2168 ## from shallow_water import Domain 2169 ### from numpy.oldnumeric import zeros, Float 2170 ## from numpy import zeros, float 2171 ## 2172 ## #Create basic mesh 2173 ## points, vertices, boundary = rectangular(1, 3) 2174 ## #print "vertices",vertices 2175 ## #Create shallow water domain 2176 ## domain = Domain(points, vertices, boundary) 2177 ## #print "domain.number_of_elements ",domain.number_of_elements 2178 ## quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2179 ## [4,4,4],[5,5,5]]) 2180 ## value = 7 2181 ## indices = [1,5] 2182 ## quantity.set_values(value, 2183 ## location = 'unique vertices', 2184 ## indices = indices) 2185 ## #print "quantity.centroid_values",quantity.centroid_values 2186 ## assert allclose(quantity.vertex_values[0], [0,7,0]) 2187 ## assert allclose(quantity.vertex_values[1], [7,1,7]) 2188 ## assert allclose(quantity.vertex_values[2], [7,2,7]) 2189 ## 2190 ## 2191 ## def test_get_values(self): 2192 ## """ 2193 ## get values based on triangle lists. 2194 ## """ 2195 ## from mesh_factory import rectangular 2196 ## from shallow_water import Domain 2197 ### from numpy.oldnumeric import zeros, Float 2198 ## from numpy import zeros, float 2199 ## 2200 ## #Create basic mesh 2201 ## points, vertices, boundary = rectangular(1, 3) 2202 ## 2203 ## #print "points",points 2204 ## #print "vertices",vertices 2205 ## #print "boundary",boundary 2206 ## 2207 ## #Create shallow water domain 2208 ## domain = Domain(points, vertices, boundary) 2209 ## #print "domain.number_of_elements ",domain.number_of_elements 2210 ## quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2211 ## [4,4,4],[5,5,5]]) 2212 ## 2213 ## #print "quantity.get_values(location = 'unique vertices')", \ 2214 ## # quantity.get_values(location = 'unique vertices') 2215 ## 2216 ## #print "quantity.get_values(location = 'unique vertices')", \ 2217 ## # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2218 ## # location = 'unique vertices') 2219 ## 2220 ## answer = [0.5,2,4,5,0,1,3,4.5] 2221 ## assert allclose(answer, 2222 ## quantity.get_values(location = 'unique vertices')) 2223 ## 2224 ## indices = [0,5,3] 2225 ## answer = [0.5,1,5] 2226 ## assert allclose(answer, 2227 ## quantity.get_values(indices=indices, \ 2228 ## location = 'unique vertices')) 2229 ## #print "quantity.centroid_values",quantity.centroid_values 2230 ## #print "quantity.get_values(location = 'centroids') ",\ 2231 ## # quantity.get_values(location = 'centroids') 2232 ## 2233 ## 2234 ## 2235 ## 2236 ## def test_get_values_2(self): 2237 ## """Different mesh (working with domain object) - also check centroids. 2238 ## """ 2239 ## 2240 ## 2241 ## a = [0.0, 0.0] 2242 ## b = [0.0, 2.0] 2243 ## c = [2.0,0.0] 2244 ## d = [0.0, 4.0] 2245 ## e = [2.0, 2.0] 2246 ## f = [4.0,0.0] 2247 ## 2248 ## points = [a, b, c, d, e, f] 2249 ## #bac, bce, ecf, dbe 2250 ## vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2251 ## 2252 ## domain = Domain(points, vertices) 2253 ## 2254 ## quantity = Quantity(domain) 2255 ## quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2256 ## 2257 ## assert allclose(quantity.get_values(location='centroids'), [2,4,4,6]) 2258 ## assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6]) 2259 ## 2260 ## 2261 ## assert allclose(quantity.get_values(location='vertices'), [[4,0,2], 2262 ## [4,2,6], 2263 ## [6,2,4], 2264 ## [8,4,6]]) 2265 ## 2266 ## assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6], 2267 ## [8,4,6]]) 2268 ## 2269 ## 2270 ## assert allclose(quantity.get_values(location='edges'), [[1,3,2], 2271 ## [4,5,3], 2272 ## [3,5,4], 2273 ## [5,7,6]]) 2274 ## assert allclose(quantity.get_values(location='edges', indices=[1,3]), 2275 ## [[4,5,3], 2276 ## [5,7,6]]) 2277 ## 2278 ## # Check averaging over vertices 2279 ## #a: 0 2280 ## #b: (4+4+4)/3 2281 ## #c: (2+2+2)/3 2282 ## #d: 8 2283 ## #e: (6+6+6)/3 2284 ## #f: 4 2285 ## assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4]) 2286 ## 2287 ## 2288 ## 2289 ## 2290 ## 2291 ## 2292 ## def test_get_interpolated_values(self): 2293 ## 2294 ## from mesh_factory import rectangular 2295 ## from shallow_water import Domain 2296 ### from numpy.oldnumeric import zeros, Float 2297 ## from numpy import zeros, float 2298 ## 2299 ## #Create basic mesh 2300 ## points, vertices, boundary = rectangular(1, 3) 2301 ## domain = Domain(points, vertices, boundary) 2302 ## 2303 ## #Constant values 2304 ## quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2305 ## [4,4,4],[5,5,5]]) 2306 ## 2307 ## 2308 ## 2309 ## # Get interpolated values at centroids 2310 ## interpolation_points = domain.get_centroid_coordinates() 2311 ## answer = quantity.get_values(location='centroids') 2312 ## 2313 ## 2314 ## #print quantity.get_values(points=interpolation_points) 2315 ## assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points)) 2316 ## 2317 ## 2318 ## #Arbitrary values 2319 ## quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7], 2320 ## [1,4,-9],[2,5,0]]) 2321 ## 2322 ## 2323 ## # Get interpolated values at centroids 2324 ## interpolation_points = domain.get_centroid_coordinates() 2325 ## answer = quantity.get_values(location='centroids') 2326 ## #print answer 2327 ## #print quantity.get_values(interpolation_points=interpolation_points) 2328 ## assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points, 2329 ## verbose=False)) 2330 ## 2331 ## 2332 ## #FIXME TODO 2333 ## #indices = [0,5,3] 2334 ## #answer = [0.5,1,5] 2335 ## #assert allclose(answer, 2336 ## # quantity.get_values(indices=indices, \ 2337 ## # location = 'unique vertices')) 2338 ## 2339 ## 2340 ## 2341 ## 2342 ## def test_get_interpolated_values_2(self): 2343 ## a = [0.0, 0.0] 2344 ## b = [0.0, 2.0] 2345 ## c = [2.0,0.0] 2346 ## d = [0.0, 4.0] 2347 ## e = [2.0, 2.0] 2348 ## f = [4.0,0.0] 2349 ## 2350 ## points = [a, b, c, d, e, f] 2351 ## #bac, bce, ecf, dbe 2352 ## vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2353 ## 2354 ## domain = Domain(points, vertices) 2355 ## 2356 ## quantity = Quantity(domain) 2357 ## quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2358 ## 2359 ## #First pick one point 2360 ## x, y = 2.0/3, 8.0/3 2361 ## v = quantity.get_values(interpolation_points = [[x,y]]) 2362 ## assert allclose(v, 6) 2363 ## 2364 ## # Then another to test that algorithm won't blindly 2365 ## # reuse interpolation matrix 2366 ## x, y = 4.0/3, 4.0/3 2367 ## v = quantity.get_values(interpolation_points = [[x,y]]) 2368 ## assert allclose(v, 4) 2369 ## 2370 ## 2371 ## 2372 ## def test_get_interpolated_values_with_georef(self): 2373 ## 2374 ## zone = 56 2375 ## xllcorner = 308500 2376 ## yllcorner = 6189000 2377 ## a = [0.0, 0.0] 2378 ## b = [0.0, 2.0] 2379 ## c = [2.0,0.0] 2380 ## d = [0.0, 4.0] 2381 ## e = [2.0, 2.0] 2382 ## f = [4.0,0.0] 2383 ## 2384 ## points = [a, b, c, d, e, f] 2385 ## #bac, bce, ecf, dbe 2386 ## vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2387 ## 2388 ## domain = Domain(points, vertices, 2389 ## geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2390 ## 2391 ## quantity = Quantity(domain) 2392 ## quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2393 ## 2394 ## #First pick one point (and turn it into absolute coordinates) 2395 ## x, y = 2.0/3, 8.0/3 2396 ## v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2397 ## assert allclose(v, 6) 2398 ## 2399 ## 2400 ## # Then another to test that algorithm won't blindly 2401 ## # reuse interpolation matrix 2402 ## x, y = 4.0/3, 4.0/3 2403 ## v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2404 ## assert allclose(v, 4) 2405 ## 2406 ## # Try two points 2407 ## pts = [[2.0/3 + xllcorner, 8.0/3 + yllcorner], 2408 ## [4.0/3 + xllcorner, 4.0/3 + yllcorner]] 2409 ## v = quantity.get_values(interpolation_points=pts) 2410 ## assert allclose(v, [6, 4]) 2411 ## 2412 ## # Test it using the geospatial data format with absolute input points and default georef 2413 ## pts = Geospatial_data(data_points=pts) 2414 ## v = quantity.get_values(interpolation_points=pts) 2415 ## assert allclose(v, [6, 4]) 2416 ## 2417 ## 2418 ## # Test it using the geospatial data format with relative input points 2419 ## pts = Geospatial_data(data_points=[[2.0/3, 8.0/3], [4.0/3, 4.0/3]], 2420 ## geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2421 ## v = quantity.get_values(interpolation_points=pts) 2422 ## assert allclose(v, [6, 4]) 2423 ## 2424 ## 2425 ## 2426 ## 2427 ## def test_getting_some_vertex_values(self): 2428 ## """ 2429 ## get values based on triangle lists. 2430 ## """ 2431 ## from mesh_factory import rectangular 2432 ## from shallow_water import Domain 2433 ### from numpy.oldnumeric import zeros, Float 2434 ## from numpy import zeros, float 2435 ## 2436 ## #Create basic mesh 2437 ## points, vertices, boundary = rectangular(1, 3) 2438 ## 2439 ## #print "points",points 2440 ## #print "vertices",vertices 2441 ## #print "boundary",boundary 2442 ## 2443 ## #Create shallow water domain 2444 ## domain = Domain(points, vertices, boundary) 2445 ## #print "domain.number_of_elements ",domain.number_of_elements 2446 ## quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2447 ## [4,4,4],[5,5,5],[6,6,6]]) 2448 ## value = [7] 2449 ## indices = [1] 2450 ## quantity.set_values(value, 2451 ## location = 'centroids', 2452 ## indices = indices) 2453 ## #print "quantity.centroid_values",quantity.centroid_values 2454 ## #print "quantity.get_values(location = 'centroids') ",\ 2455 ## # quantity.get_values(location = 'centroids') 2456 ## assert allclose(quantity.centroid_values, 2457 ## quantity.get_values(location = 'centroids')) 2458 ## 2459 ## 2460 ## value = [[15,20,25]] 2461 ## quantity.set_values(value, indices = indices) 2462 ## #print "1 quantity.vertex_values",quantity.vertex_values 2463 ## assert allclose(quantity.vertex_values, quantity.get_values()) 2464 ## 2465 ## assert allclose(quantity.edge_values, 2466 ## quantity.get_values(location = 'edges')) 2467 ## 2468 ## # get a subset of elements 2469 ## subset = quantity.get_values(location='centroids', indices=[0,5]) 2470 ## answer = [quantity.centroid_values[0],quantity.centroid_values[5]] 2471 ## assert allclose(subset, answer) 2472 ## 2473 ## 2474 ## subset = quantity.get_values(location='edges', indices=[0,5]) 2475 ## answer = [quantity.edge_values[0],quantity.edge_values[5]] 2476 ## #print "subset",subset 2477 ## #print "answer",answer 2478 ## assert allclose(subset, answer) 2479 ## 2480 ## subset = quantity.get_values( indices=[1,5]) 2481 ## answer = [quantity.vertex_values[1],quantity.vertex_values[5]] 2482 ## #print "subset",subset 2483 ## #print "answer",answer 2484 ## assert allclose(subset, answer) 2485 ## 2486 ## def test_smooth_vertex_values(self): 2487 ## """ 2488 ## get values based on triangle lists. 2489 ## """ 2490 ## from mesh_factory import rectangular 2491 ## from shallow_water import Domain 2492 ### from numpy.oldnumeric import zeros, Float 2493 ## from numpy import zeros, float 2494 ## 2495 ## #Create basic mesh 2496 ## points, vertices, boundary = rectangular(2, 2) 2497 ## 2498 ## #print "points",points 2499 ## #print "vertices",vertices 2500 ## #print "boundary",boundary 2501 ## 2502 ## #Create shallow water domain 2503 ## domain = Domain(points, vertices, boundary) 2504 ## #print "domain.number_of_elements ",domain.number_of_elements 2505 ## quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2506 ## [4,4,4],[5,5,5],[6,6,6],[7,7,7]]) 2507 ## 2508 ## #print "quantity.get_values(location = 'unique vertices')", \ 2509 ## # quantity.get_values(location = 'unique vertices') 2510 ## 2511 ## #print "quantity.get_values(location = 'unique vertices')", \ 2512 ## # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2513 ## # location = 'unique vertices') 2514 ## 2515 ## #print quantity.get_values(location = 'unique vertices') 2516 ## #print quantity.domain.number_of_triangles_per_node 2517 ## #print quantity.vertex_values 2518 ## 2519 ## #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5] 2520 ## #assert allclose(answer, 2521 ## # quantity.get_values(location = 'unique vertices')) 2522 ## 2523 ## quantity.smooth_vertex_values() 2524 ## 2525 ## #print quantity.vertex_values 2526 ## 2527 ## 2528 ## answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4], 2529 ## [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]] 2530 ## 2531 ## assert allclose(answer_vertex_values, 2532 ## quantity.vertex_values) 2533 ## #print "quantity.centroid_values",quantity.centroid_values 2534 ## #print "quantity.get_values(location = 'centroids') ",\ 2535 ## # quantity.get_values(location = 'centroids') 2519 2536 2520 2537 -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_region.py
r5211 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 … … 7 9 from region import * 8 10 #from anuga.config import epsilon 9 from Numeric import allclose, average #, array, ones, Float11 from numpy.oldnumeric import allclose, average #, array, ones, Float 10 12 """ 11 13 This is what the mesh in these tests look like; … … 44 46 from mesh_factory import rectangular 45 47 from shallow_water import Domain 46 from Numeric import zeros, Float48 from numpy.oldnumeric import zeros, Float 47 49 48 50 #Create basic mesh … … 129 131 from mesh_factory import rectangular 130 132 from shallow_water import Domain 131 from Numeric import zeros, Float133 from numpy.oldnumeric import zeros, Float 132 134 133 135 #Create basic mesh … … 162 164 from mesh_factory import rectangular 163 165 from shallow_water import Domain 164 from Numeric import zeros, Float166 from numpy.oldnumeric import zeros, Float 165 167 166 168 #Create basic mesh … … 194 196 from mesh_factory import rectangular 195 197 from shallow_water import Domain 196 from Numeric import zeros, Float198 from numpy.oldnumeric import zeros, Float 197 199 198 200 #Create basic mesh … … 229 231 from mesh_factory import rectangular 230 232 from shallow_water import Domain 231 from Numeric import zeros, Float233 from numpy.oldnumeric import zeros, Float 232 234 233 235 #Create basic mesh … … 264 266 #------------------------------------------------------------- 265 267 if __name__ == "__main__": 266 suite = unittest.makeSuite(Test_Region,'test') 268 ## suite = unittest.makeSuite(Test_Region,'test') 269 suite = unittest.makeSuite(Test_Region,'test_unique_vertices_average_loc_unique_vert') 267 270 runner = unittest.TextTestRunner() 268 271 runner.run(suite) -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/test_util.py
r5657 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 #!/usr/bin/env python 2 4 3 5 4 6 import unittest 5 from Numeric import zeros, array, allclose, Float 7 ##from numpy.oldnumeric import zeros, array, allclose, Float, alltrue 8 from numpy import zeros, array, allclose, float, alltrue 6 9 from math import sqrt, pi 7 10 import tempfile, os … … 125 128 from shallow_water import Domain, Dirichlet_boundary 126 129 from mesh_factory import rectangular 127 from Numeric import take, concatenate, reshape130 from numpy.oldnumeric import take, concatenate, reshape 128 131 129 132 #Create basic mesh and shallow water domain … … 326 329 from shallow_water import Domain, Dirichlet_boundary 327 330 from mesh_factory import rectangular 328 from Numeric import take, concatenate, reshape331 from numpy.oldnumeric import take, concatenate, reshape 329 332 330 333 … … 542 545 import os, time 543 546 from anuga.config import time_format 544 from Numeric import sin, pi, exp547 from numpy.oldnumeric import sin, pi, exp 545 548 from mesh_factory import rectangular 546 549 from shallow_water import Domain … … 703 706 import os, time 704 707 from anuga.config import time_format 705 from Numeric import sin, pi, exp708 from numpy.oldnumeric import sin, pi, exp 706 709 from mesh_factory import rectangular 707 710 from shallow_water import Domain … … 1092 1095 #This must be caught. 1093 1096 foo = array([[1,2,3], 1094 [4,5,6]], Float)1097 [4,5,6]], float) 1095 1098 1096 1099 bar = array([[-1,0,5], 1097 [6,1,1]], Float)1100 [6,1,1]], float) 1098 1101 1099 1102 D = {'X': foo, 'Y': bar} … … 1268 1271 tris = [[0,1,2]] 1269 1272 new_verts, new_tris = remove_lone_verts(verts, tris) 1270 assert new_verts == verts1271 assert new_tris == tris1273 assert alltrue(new_verts == verts) 1274 assert alltrue(new_tris == tris) 1272 1275 1273 1276 … … 1284 1287 new_verts, new_tris = remove_lone_verts(verts, tris) 1285 1288 #print "new_verts", new_verts 1286 assert new_verts == [[0,0],[1,0],[0,1]]1287 assert new_tris == [[0,1,2]]1289 assert alltrue(new_verts == [[0,0],[1,0],[0,1]]) 1290 assert alltrue(new_tris == [[0,1,2]]) 1288 1291 1289 1292 def test_remove_lone_verts_c(self): … … 1291 1294 tris = [[0,1,3]] 1292 1295 new_verts, new_tris = remove_lone_verts(verts, tris) 1293 #print "new_verts", new_verts1294 assert new_verts == [[0,0],[1,0],[0,1]]1295 assert new_tris == [[0,1,2]]1296 print "new_verts", new_verts 1297 assert alltrue(new_verts == [[0,0],[1,0],[0,1]]) 1298 assert alltrue(new_tris == [[0,1,2]]) 1296 1299 1297 1300 def test_remove_lone_verts_b(self): … … 1299 1302 tris = [[0,1,2]] 1300 1303 new_verts, new_tris = remove_lone_verts(verts, tris) 1301 assert new_verts == verts[0:3]1302 assert new_tris == tris1304 assert alltrue(new_verts == verts[0:3]) 1305 assert alltrue(new_tris == tris) 1303 1306 1304 1307 … … 1307 1310 tris = [[0,1,2]] 1308 1311 new_verts, new_tris = remove_lone_verts(verts, tris) 1309 assert new_verts == verts[0:3]1310 assert new_tris == tris1312 assert alltrue(new_verts == verts[0:3]) 1313 assert alltrue(new_tris == tris) 1311 1314 1312 1315 def test_get_min_max_values(self): -
anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/util.py
r5732 r5899 1 ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py 2 1 3 """This module contains various auxiliary function used by pyvolution. 2 4 … … 15 17 16 18 from anuga.utilities.numerical_tools import ensure_numeric 17 from Numeric import arange, choose, zeros, Float, array, allclose, take, compress19 from numpy.oldnumeric import arange, choose, zeros, Float, array, allclose, take, compress 18 20 19 21 from anuga.geospatial_data.geospatial_data import ensure_absolute … … 237 239 from anuga.config import time_format 238 240 from Scientific.IO.NetCDF import NetCDFFile 239 from Numeric import array, zeros, Float, alltrue, concatenate, reshape241 from numpy.oldnumeric import array, zeros, Float, alltrue, concatenate, reshape 240 242 241 243 # Open NetCDF file … … 303 305 # Get variables 304 306 # if verbose: print 'Get variables' 305 time = fid.variables['time'][:]307 time = array(fid.variables['time'][:] ) 306 308 307 309 # Get time independent stuff … … 373 375 # move time back - relative to domain's time 374 376 if domain_starttime > starttime: 377 ## print 'type(time)=%s' % type(time) 378 ## print 'type(domain_starttime)=%s' % type(domain_starttime) 379 ## print 'type(starttime)=%s' % type(starttime) 380 ## a = time - domain_starttime 375 381 time = time - domain_starttime + starttime 376 382 … … 1036 1042 """ 1037 1043 # from math import sqrt, atan, degrees 1038 from Numeric import ones, allclose, zeros, Float, ravel1044 from numpy.oldnumeric import ones, allclose, zeros, Float, ravel 1039 1045 from os import sep, altsep, getcwd, mkdir, access, F_OK, environ 1040 1046 … … 2124 2130 from csv import reader,writer 2125 2131 from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN 2126 from Numeric import array, resize, shape, Float, zeros, take, argsort, argmin2132 from numpy.oldnumeric import array, resize, shape, Float, zeros, take, argsort, argmin 2127 2133 import string 2128 2134 from anuga.shallow_water.data_manager import get_all_swwfiles
Note: See TracChangeset
for help on using the changeset viewer.
