# Changeset 3945

Ignore:
Timestamp:
Nov 8, 2006, 5:35:23 PM (17 years ago)
Message:

One large step towards major cleanup. This has mainly to do with
the way vertex coordinates are handled internally.

Location:
anuga_core/source/anuga
Files:
12 edited

Unmodified
Removed
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

 r3928 class General_mesh: """Collection of triangular elements (purely geometric) """Collection of 2D triangular elements A triangular element is defined in terms of three vertex ids, ordered counter clock-wise, each corresponding to a given coordinate set. Vertices from different elements can point to the same coordinate set. Coordinate sets are implemented as an N x 2 Numeric array containing ordered counter clock-wise, each corresponding to a given node which is represented as a coordinate set (x,y). Vertices from different triangles can point to the same node. The nodes are implemented as an Nx2 Numeric array containing the x and y coordinates. To instantiate: Mesh(coordinates, triangles) Mesh(nodes, triangles) where coordinates is either a list of 2-tuples or an Mx2 Numeric array of nodes is either a list of 2-tuples or an Nx2 Numeric array of floats representing all x, y coordinates in the mesh. triangles is either a list of 3-tuples or an Nx3 Numeric array of triangles is either a list of 3-tuples or an Mx3 Numeric array of integers representing indices of all vertices in the mesh. Each vertex is identified by its index i in [0, M-1]. Each vertex is identified by its index i in [0, N-1]. Example: a = [0.0, 0.0] b = [0.0, 2.0] e = [2.0, 2.0] points = [a, b, c, e] triangles = [ [1,0,2], [1,2,3] ]   #bac, bce mesh = Mesh(points, triangles) #creates two triangles: bac and bce nodes = [a, b, c, e] triangles = [ [1,0,2], [1,2,3] ]   # bac, bce # Create mesh with two triangles: bac and bce mesh = Mesh(nodes, triangles) Other: In addition mesh computes an Nx6 array called vertex_coordinates. In addition mesh computes an Mx6 array called vertex_coordinates. This structure is derived from coordinates and contains for each triangle the three x,y coordinates at the vertices. This is a cut down version of mesh from mesh.py See neighbourmesh.py for a specialisation of the general mesh class which includes information about neighbours and the mesh boundary. The mesh object is purely geometrical and contains no information about quantities defined on the mesh. """ #FIXME: It would be a good idea to use geospatial data as an alternative #input def __init__(self, coordinates, triangles, def __init__(self, nodes, triangles, geo_reference=None, number_of_full_nodes=None, number_of_full_triangles=None, geo_reference=None, verbose=False): """ Build triangles from x,y coordinates (sequence of 2-tuples or Mx2 Numeric array of floats) and triangles (sequence of 3-tuples or Nx3 Numeric array of non-negative integers). origin is a 3-tuple consisting of UTM zone, easting and northing. If specified coordinates are assumed to be relative to this origin. """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' self.triangles = array(triangles,Int) self.coordinates = array(coordinates,Float) self.triangles = array(triangles, Int) self.nodes = array(nodes, Float) # Register number of elements and nodes self.number_of_triangles = N = self.triangles.shape[0] self.number_of_nodes = self.coordinates.shape[0] self.number_of_nodes = self.nodes.shape[0] if number_of_full_nodes is None: self.number_of_full_nodes=self.number_of_nodes self.number_of_full_nodes = self.number_of_nodes 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) self.number_of_full_triangles=number_of_full_triangles self.number_of_full_triangles = number_of_full_triangles # Input checks msg = 'Triangles must an Nx3 Numeric array or a sequence of 3-tuples. ' 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) assert len(self.triangles.shape) == 2, msg msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples' msg = 'Nodes must an Nx2 Numeric array or a sequence of 2-tuples' msg += 'The supplied array has the shape: %s'\ %str(self.coordinates.shape) assert len(self.coordinates.shape) == 2, msg %str(self.nodes.shape) assert len(self.nodes.shape) == 2, msg msg = 'Vertex indices reference non-existing coordinate sets' assert max(max(self.triangles)) <= self.coordinates.shape[0], msg assert max(max(self.triangles)) <= self.nodes.shape[0], msg # FIXME: Maybe move to statistics? # Or use with get_extent xy_extent = [ min(self.coordinates[:,0]), min(self.coordinates[:,1]) , max(self.coordinates[:,0]), max(self.coordinates[:,1]) ] xy_extent = [ min(self.nodes[:,0]), min(self.nodes[:,1]) , max(self.nodes[:,0]), max(self.nodes[:,1]) ] self.xy_extent = array(xy_extent, Float) if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N) x0 = V[i, 0]; y0 = V[i, 1] x1 = V[i, 2]; y1 = V[i, 3] x2 = V[i, 4]; y2 = V[i, 5] x0, y0 = V[3*i, :] x1, y1 = V[3*i+1, :] x2, y2 = V[3*i+2, :] # Area def __repr__(self): return 'Mesh: %d vertices, %d triangles'\ %(self.coordinates.shape[0], len(self)) %(self.nodes.shape[0], len(self)) def get_normals(self): def get_vertex_coordinates(self, unique=False, obj=False, absolute=False): """Return all vertex coordinates. Return all vertex coordinates for all triangles as an Nx6 array (ordered as x0, y0, x1, y1, x2, y2 for each triangle) if obj is True, the x/y pairs are returned in a 3*N x 2 array. FIXME, we might make that the default. FIXME Maybe use keyword: continuous = False for this condition? FIXME - Maybe use something referring to unique vertices? def get_nodes(self, absolute=False): """Return all node coordinates ordered in an Nx2 array. This is the same format they were provided in the constructor i.e. without any duplication. 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. (To see which, switch to default absolute=True and run tests). (To see which, switch to default absolute=True and run tests). """ N = self.number_of_full_nodes if unique is True: V = self.coordinates[:N,:] if absolute is True: if not self.geo_reference.is_absolute(): V = self.geo_reference.get_absolute(V) return V V = self.nodes[:N,:] if absolute is True: if not self.geo_reference.is_absolute(): V = self.geo_reference.get_absolute(V) return V def get_vertex_coordinates(self, unique=False, absolute=False): """Return all vertex coordinates. Return all vertex coordinates for all triangles as a 3*N x 2 array where the jth vertex of the ith triangle is located in row 3*i+j. Boolean keyword unique will cause the points to be returned as they were provided in the constructor i.e. without any duplication in an N x 2 array. """ if unique is True: return self.get_nodes(absolute) V = self.vertex_coordinates if absolute is True: if not self.geo_reference.is_absolute(): V = self.geo_reference.get_absolute(V) V0 = self.geo_reference.get_absolute(V[:N,0:2]) V1 = self.geo_reference.get_absolute(V[:N,2:4]) V2 = self.geo_reference.get_absolute(V[:N,4:6]) # This does double the memory need V = concatenate( (V0, V1, V2), axis=1 ) if obj is True: N = V.shape[0] return reshape(V, (3*N, 2)) else: return V return V V = self.get_vertex_coordinates(absolute=absolute) return V[i, 2*j:2*j+2] ##return self.vertex_coordinates[i, 2*j:2*j+2] return V[3*i+j, :] def compute_vertex_coordinates(self): """Return vertex coordinates for all triangles as an Nx6 array (ordered as x0, y0, x1, y1, x2, y2 for each triangle) """ #FIXME (Ole): Perhaps they should be ordered as in obj files?? #See quantity.get_vertex_values #FIXME (Ole) - oh yes they should """Return all vertex coordinates for all triangles as a 3*N x 2 array where the jth vertex of the ith triangle is located in row 3*i+j. This function is used to precompute this important structure. Use get_vertex coordinates to retrieve the points. """ N = self.number_of_triangles vertex_coordinates = zeros((N, 6), Float) vertex_coordinates = zeros((3*N, 2), Float) for i in range(N): for j in range(3): k = self.triangles[i,j]  #Index of vertex 0 v_k = self.coordinates[k] vertex_coordinates[i, 2*j+0] = v_k[0] vertex_coordinates[i, 2*j+1] = v_k[1] v_k = self.nodes[k] vertex_coordinates[3*i+j,:] = v_k return vertex_coordinates def get_vertices(self, indices=None): Preconditions: self.coordinates and self.triangles are defined self.nodes and self.triangles are defined Postcondition: """ vertexlist = [None]*len(self.coordinates) vertexlist = [None]*len(self.nodes) for i in range(self.number_of_triangles):
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

 r3931 for i, (vol_id, edge_id) in enumerate(boundary_keys): x0 = V[vol_id, 0]; y0 = V[vol_id, 1] x1 = V[vol_id, 2]; y1 = V[vol_id, 3] x2 = V[vol_id, 4]; y2 = V[vol_id, 5] base_index = 3*vol_id x0, y0 = V[base_index, :] x1, y1 = V[base_index+1, :] x2, y2 = V[base_index+2, :] #Compute midpoints if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

 r3928 if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N) x0 = V[i, 0]; y0 = V[i, 1] x1 = V[i, 2]; y1 = V[i, 3] x2 = V[i, 4]; y2 = V[i, 5] x0, y0 = V[3*i, :] x1, y1 = V[3*i+1, :] x2, y2 = V[3*i+2, :] #x0 = V[i, 0]; y0 = V[i, 1] #x1 = V[i, 2]; y1 = V[i, 3] #x2 = V[i, 4]; y2 = V[i, 5] #Compute centroid #Check each triangle for i in range(N): x0, y0 = V[3*i, :] x1, y1 = V[3*i+1, :] x2, y2 = V[3*i+2, :] x0 = V[i, 0]; y0 = V[i, 1] x1 = V[i, 2]; y1 = V[i, 3] x2 = V[i, 4]; y2 = V[i, 5] #Check that area hasn't been compromised area = self.areas[i]
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

 r3928 from Numeric import array, zeros, Float, less, concatenate, NewAxis, argmax, allclose from anuga.utilities.numerical_tools import ensure_numeric from anuga.utilities.numerical_tools import ensure_numeric, is_scalar class Quantity: is_subset = True # FIXME (Ole): Now we can compute the arrays once and for all # for both centroids and vertices and the use set_values_from_array if location == 'centroids': P = take(self.domain.centroid_coordinates, indices) V = take(self.domain.centroid_coordinates, indices) if is_subset: self.set_values(f(P[:,0], P[:,1]), self.set_values(f(V[:,0], V[:,1]), location = location, indices = indices) else: self.set_values(f(P[:,0], P[:,1]), location = location) self.set_values(f(V[:,0], V[:,1]), location = location) elif location == 'vertices': P = self.domain.vertex_coordinates V = self.domain.get_vertex_coordinates() x = V[:,0]; y = V[:,1]; values = f(x, y) if is_scalar(values): # Function returned a constant value self.set_values_from_constant(values, location, indices, verbose) return if is_subset: #Brute force for e in indices: for i in range(3): self.vertex_values[e,i] = f(P[e,2*i], P[e,2*i+1]) for i in indices: for j in range(3): self.vertex_values[i,j] = values[3*i+j] else: for i in range(3): self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1]) for j in range(3): self.vertex_values[:,j] = values[j::3] else: raise 'Not implemented: %s' %location raise ms coordinates = self.domain.coordinates triangles = self.domain.triangles coordinates = self.domain.get_nodes() triangles = self.domain.triangles      #FIXME elif location == 'unique vertices': if (indices ==  None): indices=range(self.domain.coordinates.shape[0]) indices=range(self.domain.number_of_nodes) vert_values = [] #Go through list of unique vertices if indices == None: assert A.shape[0] == self.domain.coordinates.shape[0] assert A.shape[0] == self.domain.get_nodes().shape[0] vertex_list = range(A.shape[0]) else: if xy is True: X = self.domain.coordinates[:,0].astype(precision) Y = self.domain.coordinates[:,1].astype(precision) X = self.domain.get_nodes()[:,0].astype(precision) Y = self.domain.get_nodes()[:,1].astype(precision) return X, Y, A, V from Numeric import zeros, Float N = len(quantity.domain) N = quantity.domain.number_of_nodes beta_w = quantity.domain.beta_w
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

 r3928 assert allclose(domain.get_vertex_coordinates(unique=True), domain.coordinates) assert allclose(domain.get_vertex_coordinates(unique=True), domain.nodes) #assert allclose(domain.get_vertex_coordinates(), ...TODO assert domain.get_vertices([0,4]) == [domain.triangles[0], domain.triangles[4]] def test_areas(self): from mesh_factory import rectangular
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

 r3688 #Vertex coordinates #V = mesh.get_vertex_coordinates() #assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0]) V = mesh.get_vertex_coordinates() assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0]) V = mesh.get_vertex_coordinates(obj=True) assert allclose(V, [ [0.0, 0.0], [4.0, 0.0], V = mesh.get_vertex_coordinates() x0 = V[0,0] y0 = V[0,1] x1 = V[0,2] y1 = V[0,3] x2 = V[0,4] y2 = V[0,5] x0 = V[0, 0]; y0 = V[0, 1] x1 = V[1, 0]; y1 = V[1, 1] x2 = V[2, 0]; y2 = V[2, 1] #x0 = V[0,0] #y0 = V[0,1] #x1 = V[0,2] #y1 = V[0,3] #x2 = V[0,4] #y2 = V[0,5] m0 = [(x1 + x2)/2, (y1 + y2)/2] V = mesh.get_vertex_coordinates() x0 = V[0,0] y0 = V[0,1] x1 = V[0,2] y1 = V[0,3] x2 = V[0,4] y2 = V[0,5] x0 = V[0, 0]; y0 = V[0, 1] x1 = V[1, 0]; y1 = V[1, 1] x2 = V[2, 0]; y2 = V[2, 1] #x0 = V[0,0] #y0 = V[0,1] #x1 = V[0,2] #y1 = V[0,3] #x2 = V[0,4] #y2 = V[0,5] m0 = [(x1 + x2)/2, (y1 + y2)/2]

• ## anuga_core/source/anuga/fit_interpolate/fit.py

 r3633 max_vertices_per_cell) m = self.mesh.coordinates.shape[0] #Nbr of basis functions (vertices) m = self.mesh.number_of_nodes # Nbr of basis functions (vertices) self.AtA = None #"element stiffness matrices: m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex) m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex) self.D = Sparse(m,m) if self.AtA == None: # AtA and Atz need ot be initialised. m = self.mesh.coordinates.shape[0] #Nbr of vertices m = self.mesh.number_of_nodes if len(z.shape) > 1: att_num = z.shape[1] #Check sanity m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex) m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex) n = self.point_count if n
• ## anuga_core/source/anuga/fit_interpolate/interpolate.py

 r3941 print '\n WARNING: No points within the mesh! \n' m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex) n = point_coordinates.shape[0]     #Nbr of data points m = self.mesh.number_of_nodes  # Nbr of basis functions (1/vertex) n = point_coordinates.shape[0] # Nbr of data points if verbose: print 'Number of datapoints: %d' %n
• ## anuga_core/source/anuga/utilities/numerical_tools.py

 r3850 if x < 0: return -1 if x == 0: return 0 def is_scalar(x): """True if x is a scalar (constant numeric value) """ from types import IntType, FloatType if type(x) in [IntType, FloatType]: return True else: return False def angle(v1, v2=None):