# Changeset 6428

Ignore:
Timestamp:
Feb 27, 2009, 11:54:09 AM (15 years ago)
Message:

numpy changes

Location:
branches/numpy/anuga
Files:
16 edited

Unmodified
Removed
• ## branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py

 r6410 triangles = [ [1,0,2], [1,2,3] ]   # bac, bce # Create mesh with two triangles: bac and bce # Create mesh with two triangles: bac and bce mesh = Mesh(nodes, triangles) """ #FIXME: It would be a good idea to use geospatial data as an alternative #input def __init__(self, nodes, triangles, geo_reference=None, # FIXME: It would be a good idea to use geospatial data as an alternative #        input def __init__(self, nodes, triangles, geo_reference=None, number_of_full_nodes=None, number_of_full_triangles=None, number_of_full_triangles=None, verbose=False): """Build triangular 2d mesh from nodes and triangle information Input: nodes: x,y coordinates represented as a sequence of 2-tuples or a Nx2 numeric array of floats. triangles: sequence of 3-tuples or Mx3 numeric array of non-negative integers representing indices into the nodes array. georeference (optional): If specified coordinates are assumed to be relative to this origin. In this case it is usefull to specify the number of real (called full) nodes and triangles. If omitted they will default to all. """ if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' """ if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' self.triangles = num.array(triangles, num.int) self.nodes = num.array(nodes, num.float) # Register number of elements and nodes # Register number of elements and nodes self.number_of_triangles = N = self.triangles.shape[0] self.number_of_nodes = self.nodes.shape[0] self.number_of_nodes = self.nodes.shape[0] if number_of_full_nodes is None: else: assert int(number_of_full_nodes) self.number_of_full_nodes = number_of_full_nodes self.number_of_full_nodes = number_of_full_nodes if number_of_full_triangles is None: self.number_of_full_triangles = self.number_of_triangles self.number_of_full_triangles = self.number_of_triangles else: assert int(number_of_full_triangles) assert int(number_of_full_triangles) self.number_of_full_triangles = number_of_full_triangles #print self.number_of_full_nodes, self.number_of_nodes #print self.number_of_full_triangles, self.number_of_triangles # FIXME: this stores a geo_reference, but when coords are returned # This geo_ref is not taken into account! if geo_reference is None: self.geo_reference = Geo_reference() #Use defaults self.geo_reference = Geo_reference()    # Use defaults else: self.geo_reference = geo_reference # Input checks msg = 'Triangles must an Mx3 numeric array or a sequence of 3-tuples. ' msg += 'The supplied array has the shape: %s'\ %str(self.triangles.shape) msg = ('Triangles must an Mx3 numeric array or a sequence of 3-tuples. ' 'The supplied array has the shape: %s' % str(self.triangles.shape)) assert len(self.triangles.shape) == 2, msg msg = 'Nodes must an Nx2 numeric array or a sequence of 2-tuples' msg += 'The supplied array has the shape: %s'\ %str(self.nodes.shape) msg = ('Nodes must an Nx2 numeric array or a sequence of 2-tuples' 'The supplied array has the shape: %s' % str(self.nodes.shape)) assert len(self.nodes.shape) == 2, msg assert max(self.triangles.flat) < self.nodes.shape[0], msg # FIXME: Maybe move to statistics? # Or use with get_extent xy_extent = [ min(self.nodes[:,0]), min(self.nodes[:,1]) , max(self.nodes[:,0]), max(self.nodes[:,1]) ] xy_extent = [min(self.nodes[:,0]), min(self.nodes[:,1]), max(self.nodes[:,0]), max(self.nodes[:,1])] self.xy_extent = num.array(xy_extent, num.float) # Allocate space for geometric quantities self.vertex_coordinates = V = self.compute_vertex_coordinates() # Initialise each triangle if verbose: print 'General_mesh: Computing areas, normals and edgelenghts' print 'General_mesh: Computing areas, normals and edgelengths' for i in range(N): if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N) if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' % (i, N) x0, y0 = V[3*i, :] x1, y1 = V[3*i+1, :] x2, y2 = V[3*i+2, :] x2, y2 = V[3*i+2, :] # Area self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2 msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2) msg += ' is degenerate:  area == %f' %self.areas[i] self.areas[i] = abs((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2 msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' % (x0,y0,x1,y1,x2,y2) msg += ' is degenerate:  area == %f' % self.areas[i] assert self.areas[i] > 0.0, msg # Normals #     the first vertex, etc) #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle n0 = num.array([x2 - x1, y2 - y1], num.float) n0 = num.array([x2-x1, y2-y1], num.float) l0 = num.sqrt(num.sum(n0**2)) n1 = num.array([x0 - x2, y0 - y2], num.float) n1 = num.array([x0-x2, y0-y2], num.float) l1 = num.sqrt(num.sum(n1**2)) n2 = num.array([x1 - x0, y1 - y0], num.float) n2 = num.array([x1-x0, y1-y0], num.float) l2 = num.sqrt(num.sum(n2**2)) self.edgelengths[i, :] = [l0, l1, l2] # Build structure listing which trianglse belong to which node. if verbose: print 'Building inverted triangle structure' if verbose: print 'Building inverted triangle structure' self.build_inverted_triangle_structure() def __len__(self): return self.number_of_triangles def __repr__(self): return 'Mesh: %d vertices, %d triangles'\ %(self.nodes.shape[0], len(self)) return ('Mesh: %d vertices, %d triangles' % (self.nodes.shape[0], len(self))) def get_normals(self): """Return all normal vectors. Return normal vectors for all triangles as an Nx6 array (ordered as x0, y0, x1, y1, x2, y2 for each triangle) """ return self.normals def get_normal(self, i, j): """Return normal vector j of the i'th triangle. Return value is the numeric array slice [x, y] """ return self.normals[i, 2*j:2*j+2] def get_number_of_nodes(self): return self.number_of_nodes def get_nodes(self, absolute=False): """Return all nodes in mesh. are to be made absolute by taking georeference into account Default is False as many parts of ANUGA expects relative coordinates. (To see which, switch to default absolute=True and run tests). (To see which, switch to default absolute=True and run tests). """ if not self.geo_reference.is_absolute(): V = self.geo_reference.get_absolute(V) return V def get_node(self, i, absolute=False): def get_node(self, i, absolute=False): """Return node coordinates for triangle i. are to be made absolute by taking georeference into account Default is False as many parts of ANUGA expects relative coordinates. (To see which, switch to default absolute=True and run tests). """ (To see which, switch to default absolute=True and run tests). """ V = self.nodes[i,:] if absolute is True: if not self.geo_reference.is_absolute(): return V + num.array([self.geo_reference.get_xllcorner(), self.geo_reference.get_yllcorner()], num.float) else: return V else: return V def get_vertex_coordinates(self, triangle_id=None, absolute=False): """Return vertex coordinates for all triangles. V += num.array([self.geo_reference.get_xllcorner(), self.geo_reference.get_yllcorner()], num.float) return V def get_vertex_coordinates(self, triangle_id=None, absolute=False): """Return vertex coordinates for all triangles. Return all vertex coordinates for all triangles as a 3*M x 2 array where the jth vertex of the ith triangle is located in row 3*i+j and if triangle_id is specified (an integer) the 3 vertex coordinates for triangle_id are returned. Boolean keyword argument absolute determines whether coordinates are to be made absolute by taking georeference into account Default is False as many parts of ANUGA expects relative coordinates. """ V = self.vertex_coordinates if triangle_id is None: if triangle_id is None: if absolute is True: if not self.geo_reference.is_absolute(): V = self.geo_reference.get_absolute(V) return V else: assert int(i) == i, msg assert 0 <= i < self.number_of_triangles i3 = 3*i i3 = 3*i if absolute is True and not self.geo_reference.is_absolute(): offset=num.array([self.geo_reference.get_xllcorner(), else: return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.float) def get_vertex_coordinate(self, i, j, absolute=False): msg = 'vertex id j must be an integer in [0,1,2]' assert j in [0,1,2], msg V = self.get_vertex_coordinates(triangle_id=i, absolute=absolute) V = self.get_vertex_coordinates(triangle_id=i, absolute=absolute) return V[j,:] def compute_vertex_coordinates(self): return vertex_coordinates def get_triangles(self, indices=None): """Get mesh triangles. if indices is None: return self.triangles #indices = range(M) return num.take(self.triangles, indices, axis=0) def get_disconnected_triangles(self): The triangles created will have the format [[0,1,2], [3,4,5], [6,7,8], ... [3*M-3 3*M-2 3*M-1]] [3*M-3 3*M-2 3*M-1]] """ M = len(self) # Number of triangles K = 3*M       # Total number of unique vertices T = num.reshape(num.arange(K, dtype=num.int), (M,3)) return T def get_unique_vertices(self,  indices=None): return num.reshape(num.arange(K, dtype=num.int), (M,3)) def get_unique_vertices(self, indices=None): """FIXME(Ole): This function needs a docstring""" unique_verts = {} for triangle in triangles: #print 'triangle(%s)=%s' % (type(triangle), str(triangle)) unique_verts[triangle[0]] = 0 unique_verts[triangle[1]] = 0 return unique_verts.keys() def get_triangles_and_vertices_per_node(self, node=None): """Get triangles associated with given node. # Get index for this node first = num.sum(self.number_of_triangles_per_node[:node]) # Get number of triangles for this node count = self.number_of_triangles_per_node[node] # working directly with the inverted triangle structure for i in range(self.number_of_full_nodes): L = self.get_triangles_and_vertices_per_node(node=i) triangle_list.append(L) return triangle_list def build_inverted_triangle_structure(self): listing for each node how many triangles use it. N is the number of nodes in mesh. vertex_value_indices: An array of length M listing indices into triangles ordered by node number. The (triangle_id, vertex_id) used to average vertex values efficiently. Example: a = [0.0, 0.0] # node 0 b = [0.0, 2.0] # node 1 e = [2.0, 2.0] # node 4 f = [4.0, 0.0] # node 5 nodes = array([a, b, c, d, e, f]) #bac, bce, ecf, dbe, daf, dae triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) For this structure #                    bac,     bce,     ecf,     dbe triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) For this structure: number_of_triangles_per_node = [1 3 3 1 3 1] which means that node a has 1 triangle associated with it, node b has 3, node has 3 and so on. vertex_value_indices = [ 1  0  3 10  2  4  7  9  5  6 11  8] which reflects the fact that and by triangle 2, vertex 0 (index = 6) and by triangle 3, vertex 2 (index = 11) node 5 is used by triangle 2, vertex 2 (index = 8) node 5 is used by triangle 2, vertex 2 (index = 8) Preconditions: Postcondition: self.number_of_triangles_per_node is built self.vertex_value_indices is built self.vertex_value_indices is built """ # Register (triangle, vertex) indices for each node vertexlist = [None]*self.number_of_full_nodes vertexlist = [None] * self.number_of_full_nodes for volume_id in range(self.number_of_full_triangles): a = self.triangles[volume_id, 0] b = self.triangles[volume_id, 1] c = self.triangles[volume_id, 2] for vertex_id, node_id in enumerate([a,b,c]): for vertex_id, node_id in enumerate([a, b, c]): if vertexlist[node_id] is None: vertexlist[node_id] = [] vertexlist[node_id].append( (volume_id, vertex_id) ) vertexlist[node_id].append((volume_id, vertex_id)) # Build inverted triangle index array for volume_id, vertex_id in vertices: vertex_value_indices[k] = 3*volume_id + vertex_id k += 1 self.vertex_value_indices = vertex_value_indices def get_extent(self, absolute=False): """Return min and max of all x and y coordinates are to be made absolute by taking georeference into account """ C = self.get_vertex_coordinates(absolute=absolute) ymin = min(Y.flat) ymax = max(Y.flat) #print "C",C return xmin, xmax, ymin, ymax def get_areas(self): """Get areas of all individual triangles. """ return self.areas """Get areas of all individual triangles.""" return self.areas def get_area(self): """Return total area of mesh """ """Return total area of mesh""" return num.sum(self.areas) def set_georeference(self, g): self.geo_reference = g def get_georeference(self): return self.geo_reference

• ## branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py

 r6410 e = [2.0, 2.0] f = [4.0, 0.0] nodes = num.array([a, b, c, d, e, f]) nodes_absolute = geo.get_absolute(nodes) #bac, bce, ecf, dbe, daf, dae #                        bac,     bce,     ecf,     dbe triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int) domain = General_mesh(nodes, triangles, geo_reference = geo) verts = domain.get_vertex_coordinates(triangle_id=0) msg = ("num.array([b,a,c])=\n%s\nshould be the same as 'verts'=\n%s" domain = General_mesh(nodes, triangles, geo_reference=geo) verts = domain.get_vertex_coordinates(triangle_id=0)    # bac msg = ("num.array([b,a,c])=\n%s\nshould be close to 'verts'=\n%s" % (str(num.array([b,a,c])), str(verts))) #self.assert_(num.allclose(num.array([b,a,c]), verts)) assert num.allclose(num.array([b,a,c]), verts), msg self.failUnless(num.allclose(num.array([b,a,c]), verts), msg) verts = domain.get_vertex_coordinates(triangle_id=0) self.assert_(num.allclose(num.array([b,a,c]), verts)) msg = ("num.array([b,a,c])=\n%s\nshould be close to 'verts'=\n%s" % (str(num.array([b,a,c])), str(verts))) self.assert_(num.allclose(num.array([b,a,c]), verts), msg) verts = domain.get_vertex_coordinates(triangle_id=0, absolute=True) msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s" % (str(num.array([nodes_absolute[1], nodes_absolute[0], nodes_absolute[2]])), str(verts))) self.assert_(num.allclose(num.array([nodes_absolute[1], nodes_absolute[0], nodes_absolute[2]]), verts), msg) verts = domain.get_vertex_coordinates(triangle_id=0, absolute=True) msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s" % (str(num.array([nodes_absolute[1], nodes_absolute[0], nodes_absolute[2]])), str(verts))) self.assert_(num.allclose(num.array([nodes_absolute[1], nodes_absolute[0], nodes_absolute[2]]), verts)) verts = domain.get_vertex_coordinates(triangle_id=0, absolute=True) self.assert_(num.allclose(num.array([nodes_absolute[1], nodes_absolute[0], nodes_absolute[2]]), verts)) nodes_absolute[0], nodes_absolute[2]]), verts), msg) def test_get_vertex_coordinates_triangle_id(self): nodes_absolute = geo.get_absolute(nodes) #bac, bce, ecf, dbe, daf, dae #                        bac,     bce,     ecf,     dbe triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) domain = General_mesh(nodes, triangles, geo_reference = geo) domain = General_mesh(nodes, triangles, geo_reference = geo) node = domain.get_node(2) #self.assertEqual(c, node) msg = ('\nc=%s\nmode=%s' % (str(c), str(node))) msg = ('\nc=%s\nnode=%s' % (str(c), str(node))) self.failUnless(num.alltrue(c == node), msg) # repeat get_node(), see if result same node = domain.get_node(2) msg = ('\nc=%s\nnode=%s' % (str(c), str(node))) self.failUnless(num.alltrue(c == node), msg) node = domain.get_node(2, absolute=True) #self.assertEqual(nodes_absolute[2], node) self.failUnless(num.alltrue(nodes_absolute[2] == node)) msg = ('\nnodes_absolute[2]=%s\nnode=%s' % (str(nodes_absolute[2]), str(node))) self.failUnless(num.alltrue(nodes_absolute[2] == node), msg) # repeat get_node(absolute=True), see if result same node = domain.get_node(2, absolute=True) #self.assertEqual(nodes_absolute[2], node) self.failUnless(num.alltrue(nodes_absolute[2] == node)) msg = ('\nnodes_absolute[2]=%s\nnode=%s' % (str(nodes_absolute[2]), str(node))) self.failUnless(num.alltrue(nodes_absolute[2] == node), msg) self.failUnlessRaises(AssertionError, General_mesh, nodes, triangles, geo_reference=geo) def test_raw_change_points_geo_ref(self): x0 = 1000.0 y0 = 2000.0 geo = Geo_reference(56, x0, y0) if __name__ == "__main__": suite = unittest.makeSuite(Test_General_Mesh,'test') suite = unittest.makeSuite(Test_General_Mesh, 'test') #suite = unittest.makeSuite(Test_General_Mesh, 'test_get_vertex_coordinates_with_geo_ref') runner = unittest.TextTestRunner() runner.run(suite)
• ## branches/numpy/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

 r6410 """ # test a = [0.0, 0.0] absolute_points = [a, b, c] vertices = [[0,1,2]] geo = Geo_reference(56,67,-56) vertices = [[0, 1, 2]] geo = Geo_reference(56, 67, -56) relative_points = geo.change_points_geo_ref(absolute_points) #print 'Relative', relative_points #print 'Absolute', absolute_points mesh = Mesh(relative_points, vertices, geo_reference=geo) if __name__ == "__main__": suite = unittest.makeSuite(Test_Mesh,'test') suite = unittest.makeSuite(Test_Mesh, 'test') runner = unittest.TextTestRunner() runner.run(suite)
• ## branches/numpy/anuga/coordinate_transforms/geo_reference.py

 r6360 self.read_ASCII(ASCIIFile, read_title=read_title) #    # Might be better to have this method instead of the following 3. #    def get_origin(self): #        return (self.zone, self.xllcorner, self.yllcorner) ## # @brief Get the X coordinate of the origin of this georef. ################################################################################ ## # @brief Change points to be absolute wrt new georef 'points_geo_ref'. # @param points The points to change. # @param points_geo_ref The new georef to make points absolute wrt. # @return The changed points. # @note If 'points' is a list then a changed list is returned. # @note The input points data is changed. def change_points_geo_ref(self, points, points_geo_ref=None): """Change the geo reference of a list or numeric array of points to be this reference.(The reference used for this object) If the points do not have a geo ref, assume 'absolute' input values """ Change the geo reference of a list or numeric array of points to be this reference.(The reference used for this object) If the points do not have a geo ref, assume 'absolute' values """ is_list = False if type(points) == types.ListType: is_list = True # remember if we got a list is_list = isinstance(points, list) points = ensure_numeric(points, num.float) # sanity checks if len(points.shape) == 1: # One point has been passed # @return True if ??? def is_absolute(self): """Return True if xllcorner==yllcorner==0 indicating that points in question are absolute. """ """Return True if xllcorner==yllcorner==0""" return num.allclose([self.xllcorner, self.yllcorner], 0) # @param points # @return # @note # @note This method changes the input points! def get_absolute(self, points): """Given a set of points geo referenced to this instance, """Given a set of points geo referenced to this instance return the points as absolute values. """ #if self.is_absolute: #    return points # remember if we got a list is_list = False if type(points) == types.ListType: is_list = True is_list = isinstance(points, list) # convert to numeric array points = ensure_numeric(points, num.float) # sanity checks if len(points.shape) == 1: #One point has been passed # Add geo ref to points #if not self.is_absolute: if not self.is_absolute(): points[:,0] += self.xllcorner
• ## branches/numpy/anuga/coordinate_transforms/test_geo_reference.py

 r6410 'bad shape did not raise error!') os.remove(point_file) def test_functionality_get_absolute(self): x0 = 1000.0 y0 = 2000.0 geo = Geo_reference(56, x0, y0) # iterable points (*not* num.array()) points = ((2,3), (3,1), (5,2)) abs_points = geo.get_absolute(points) # check we haven't changed 'points' itself self.failIf(num.alltrue(abs_points == points)) new_points = abs_points.copy() new_points[:,0] -= x0 new_points[:,1] -= y0 self.failUnless(num.alltrue(new_points == points)) # points in num.array() points = num.array(((2,3), (3,1), (5,2)), num.float) abs_points = geo.get_absolute(points) # check we haven't changed 'points' itself self.failIf(num.alltrue(abs_points == points)) new_points = abs_points.copy() new_points[:,0] -= x0 new_points[:,1] -= y0 self.failUnless(num.alltrue(new_points == points)) #------------------------------------------------------------- if __name__ == "__main__": suite = unittest.makeSuite(geo_referenceTestCase,'test') suite = unittest.makeSuite(geo_referenceTestCase, 'test') #suite = unittest.makeSuite(geo_referenceTestCase, 'test_functionality_get_absolute') runner = unittest.TextTestRunner() #verbosity=2) runner.run(suite)
• ## branches/numpy/anuga/fit_interpolate/interpolate.py

 r6410 t = self.interpolate_block(f, point_coordinates[start:end], verbose=verbose) z = num.concatenate((z, t)) z = num.concatenate((z, t), axis=0)    #??default# start = end t = self.interpolate_block(f, point_coordinates[start:end], verbose=verbose) z = num.concatenate((z, t)) z = num.concatenate((z, t), axis=0)    #??default# return z #Add the x and y together vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]),axis=1) vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]), axis=1) #Will return the quantity values at the specified times and locations

• ## branches/numpy/anuga/geospatial_data/test_geospatial_data.py

 r6410 P = G.get_data_points(absolute=True) assert num.allclose(P, num.concatenate( (P1,P2) )) assert num.allclose(P, num.concatenate((P1,P2), axis=0))    #??default# msg = ('P=\n%s\nshould be close to\n' '[[2.0, 4.1], [4.0, 7.3],\n' P = G.get_data_points(absolute=True) assert num.allclose(P, num.concatenate((points1, points2))) assert num.allclose(P, num.concatenate((points1, points2), axis=0)) #??default# def test_add_with_None(self):

 r6410 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a from anuga.config import netcdf_float, netcdf_char, netcdf_int from anuga.utilities.system_tools import * from Scientific.IO.NetCDF import NetCDFFile num.array(mesh['vertex_attributes'], num.float) mesh['vertex_attribute_titles'] = \ num.array(mesh['vertex_attribute_titles'], num.character) num.array(string_to_char(mesh['vertex_attribute_titles']), num.character) mesh['segments'] = num.array(mesh['segments'], IntType) mesh['segment_tags'] = num.array(mesh['segment_tags'], num.character) mesh['segment_tags'] = num.array(string_to_char(mesh['segment_tags']), num.character) mesh['triangles'] = num.array(mesh['triangles'], IntType) mesh['triangle_tags'] = num.array(mesh['triangle_tags']) mesh['triangle_tags'] = num.array(string_to_char(mesh['triangle_tags']), num.character) mesh['triangle_neighbors'] = \ num.array(mesh['triangle_neighbors'], IntType) mesh['outline_segments'] = num.array(mesh['outline_segments'], IntType) mesh['outline_segment_tags'] = \ num.array(mesh['outline_segment_tags'], num.character) num.array(string_to_char(mesh['outline_segment_tags']), num.character) mesh['holes'] = num.array(mesh['holes'], num.float) mesh['regions'] = num.array(mesh['regions'], num.float) mesh['region_tags'] = num.array(mesh['region_tags'], num.character) mesh['region_tags'] = num.array(string_to_char(mesh['region_tags']), num.character) mesh['region_max_areas'] = num.array(mesh['region_max_areas'], num.float) ('num_of_segments', 'num_of_segment_ends')) outfile.variables['segments'][:] = mesh['segments'] if mesh['segment_tags'].shape[1] > 0: if mesh['segment_tags'].shape[0] > 0: outfile.createDimension('num_of_segment_tag_chars', mesh['segment_tags'].shape[1]) 'num_of_segment_ends')) outfile.variables['outline_segments'][:] = mesh['outline_segments'] if (mesh['outline_segment_tags'].shape[1] > 0): if mesh['outline_segment_tags'].shape[1] > 0: outfile.createDimension('num_of_outline_segment_tag_chars', mesh['outline_segment_tags'].shape[1]) try: titles = fid.variables['vertex_attribute_titles'][:] for i, title in enumerate(titles): mesh['vertex_attribute_titles'].append(titles[i].tostring().strip()) except KeyError: mesh['vertex_attribute_titles'] = [x.tostring().strip() for x in titles] except: pass mesh['segments'] = num.array([], num.int)      #array default# mesh['segment_tags'] =[] mesh['segment_tags'] = [] try: tags = fid.variables['segment_tags'][:] for i, tag in enumerate(tags): mesh['segment_tags'].append(tags[i].tostring().strip()) except KeyError: for ob in mesh['segments']: mesh['segment_tags'].append('') mesh['segment_tags'] = [x.tostring().strip() for x in tags] except: pass try: try: tags = fid.variables['triangle_tags'][:] for i, tag in enumerate(tags): mesh['triangle_tags'].append(tags[i].tostring().strip()) mesh['triangle_tags'] = [x.tostring().strip() for x in tags] except KeyError: for ob in mesh['triangles']:

• ## branches/numpy/anuga/shallow_water/test_data_manager.py

 r6410 # Invoke interpolation for vertex points points = num.concatenate( (x[:,num.newaxis],y[:,num.newaxis]), axis=1 ) points = num.ascontiguousarray(points) sww2pts(self.domain.get_name(), quantity = 'elevation', domain.set_quantity('xmomentum', uh) domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br}) for t in domain.evolve(yieldstep=1, finaltime = t_end): pass if __name__ == "__main__": suite = unittest.makeSuite(Test_Data_Manager,'test') #suite = unittest.makeSuite(Test_Data_Manager,'test_get_maximum_inundation') # FIXME (Ole): This is the test that fails under Windows
• ## branches/numpy/anuga/shallow_water/test_smf.py

 r6410 slide = slide_tsunami(length=len, depth=dep, slope=th, x0=x0, \ width = wid, thickness=thk, kappa=kappa, kappad=kappad, \ verbose=False) verbose=False) assert num.allclose(slide.a3D, 0.07775819) slide = slide_tsunami(length, dep, th, x0, y0, \ wid, thk, kappa, kappad, \ domain=domain,verbose=False) domain=domain,verbose=False) domain.set_quantity('stage', slide) stage = domain.get_quantity('stage') w = stage.get_values() stage = domain.get_quantity('stage') w = stage.get_values() ##      check = [[-0.0 -0.0 -0.0], ##        check = [[-0.0 -0.0 -0.0], ##                 [-.189709745 -517.877716 -0.0], ##                 [-0.0 -0.0 -2.7695931e-08], ##                 [-0.0 -0.0 -0.0]] assert num.allclose(min(min(w)), -517.877771593) assert num.allclose(max(max(w)), 0.0) assert num.allclose(num.min(w), -517.877771593) assert num.allclose(num.max(w), 0.0) assert num.allclose(slide.a3D, 518.38797486)
• ## branches/numpy/anuga/utilities/numerical_tools.py

 r6304 n = num.searchsorted(num.sort(a), bins) n = num.concatenate( [n, [len(a)]] ) n = num.concatenate([n, [len(a)]], axis=0)    #??default# hist = n[1:]-n[:-1]
• ## branches/numpy/anuga/utilities/system_tools.py

 r6415 '''Convert 1-D list of strings to 2-D list of chars.''' if not l: return [] if l == ['']: l = [' '] maxlen = reduce(max, map(len, l)) ll = [x.ljust(maxlen) for x in l]
• ## branches/numpy/anuga/utilities/test_system_tools.py

 r6415 import numpy as num import zlib import os from os.path import join, split, sep from Scientific.IO.NetCDF import NetCDFFile from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a from anuga.config import netcdf_float from anuga.config import netcdf_float, netcdf_char, netcdf_int os.remove(tmp_name) os.remove(tmp_name) # Binary NetCDF File X 2 (use mktemp's name) fid.variables['test_array'][:] = test_array fid.close() # Second file filename2 = mktemp(suffix='.nc', dir='.') fid.variables['test_array'][:] = test_array fid.close() checksum1 = compute_checksum(filename1) checksum2 = compute_checksum(filename2) if path == '': path = '.' + sep filename = path + sep +  'crc_test_file.png' MAX_CHARS = 10 MAX_ENTRIES = 100000 MAX_ENTRIES = 10000 A_INT = ord('a') Z_INT = ord('z') self.failUnlessEqual(new_str_list, str_list) # special test - input list is [''] def test_string_to_char2(self): # generate a special list shown bad in load_mesh testing str_list = [''] x = string_to_char(str_list) new_str_list = char_to_string(x) self.failUnlessEqual(new_str_list, str_list) ################################################################################ ################################################################################ ## # @brief Helper function to write a list of strings to a NetCDF file. # @param filename Path to the file to write. # @param l The list of strings to write. def helper_write_msh_file(self, filename, l): # open the NetCDF file fd = NetCDFFile(filename, netcdf_mode_w) fd.description = 'Test file - string arrays' # convert list of strings to num.array al = num.array(string_to_char(l), num.character) # write the list fd.createDimension('num_of_strings', al.shape[0]) fd.createDimension('size_of_strings', al.shape[1]) var = fd.createVariable('strings', netcdf_char, ('num_of_strings', 'size_of_strings')) var[:] = al fd.close() ## # @brief Helper function to read a NetCDF file and return a list of strings. # @param filename Path to the file to read. # @return A list of strings from the file. def helper_read_msh_file(self, filename): fid = NetCDFFile(filename, netcdf_mode_r) mesh = {} # Get the 'strings' variable strings = fid.variables['strings'][:] fid.close() return char_to_string(strings) # test random strings to a NetCDF file def test_string_to_netcdf(self): import random MAX_CHARS = 10 MAX_ENTRIES = 10000 A_INT = ord('a') Z_INT = ord('z') FILENAME = 'test.msh' # generate some random strings in a list, with guaranteed lengths str_list = ['x' * MAX_CHARS]        # make first maximum length for entry in xrange(MAX_ENTRIES): length = random.randint(1, MAX_CHARS) s = '' for c in range(length): s += chr(random.randint(A_INT, Z_INT)) str_list.append(s) self.helper_write_msh_file(FILENAME, str_list) new_str_list = self.helper_read_msh_file(FILENAME) self.failUnlessEqual(new_str_list, str_list) os.remove(FILENAME) # special test - list [''] to a NetCDF file def test_string_to_netcdf2(self): FILENAME = 'test.msh' # generate some random strings in a list, with guaranteed lengths str_list = [''] self.helper_write_msh_file(FILENAME, str_list) new_str_list = self.helper_read_msh_file(FILENAME) self.failUnlessEqual(new_str_list, str_list) os.remove(FILENAME) #        print ('test_string_to_char: str_list=%s, new_str_list=%s' #               % (str(str_list), str(new_str_list))) ################################################################################ if __name__ == "__main__": #suite = unittest.makeSuite(Test_system_tools, 'test') suite = unittest.makeSuite(Test_system_tools, 'test_string_to_char') suite = unittest.makeSuite(Test_system_tools, 'test') runner = unittest.TextTestRunner() runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.