# Changeset 3954

Ignore:
Timestamp:
Nov 9, 2006, 10:41:21 AM (17 years ago)
Message:

More cleanup of triangle and node formats + better comments

Location:
anuga_core/source
Files:
6 edited

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

 r3945 def get_nodes(self, absolute=False): """Return all node coordinates ordered in an Nx2 array. """Return all nodes in mesh. The nodes are ordered in an Nx2 array where N is the number of nodes. This is the same format they were provided in the constructor i.e. without any duplication. 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 def get_vertex_coordinates(self, 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 M the number of triangles in the mesh. 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 #[:3*M,:] if absolute is True: if not self.geo_reference.is_absolute(): def compute_vertex_coordinates(self): """Return all vertex coordinates for all triangles as a 3*N x 2 array """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. """ N = self.number_of_triangles vertex_coordinates = zeros((3*N, 2), Float) for i in range(N): M = self.number_of_triangles vertex_coordinates = zeros((3*M, 2), Float) for i in range(M): for j in range(3): k = self.triangles[i,j]  #Index of vertex 0 v_k = self.nodes[k] vertex_coordinates[3*i+j,:] = v_k k = self.triangles[i,j] # Index of vertex j in triangle i vertex_coordinates[3*i+j,:] = self.nodes[k] return vertex_coordinates def get_vertices(self, indices=None): """Get connectivity indices is the set of element ids of interest """ N = self.number_of_full_triangles def get_triangles(self, indices=None): """Get mesh triangles. Return Mx3 integer array where M is the number of triangles. Each row corresponds to one triangle and the three entries are indices into the mesh nodes which can be obtained using the method get_nodes() Optional argument, indices is the set of triangle ids of interest. """ M = self.number_of_full_triangles if indices is None: #indices = range(len(self))  #len(self)=number of elements indices = range(N) return  take(self.triangles, indices) indices = range(M) return take(self.triangles, indices) #FIXME - merge these two (get_vertices and get_triangles) def get_triangles(self, obj=False): """Get connetivity Return triangles (triplets of indices into point coordinates) If obj is True return structure commensurate with replicated points, allowing for discontinuities (FIXME: Need good name for this concept) """ if obj is True: m = len(self)  #Number of triangles M = 3*m        #Total number of unique vertices T = reshape(array(range(M)).astype(Int), (m,3)) else: T = self.triangles def get_disconnected_triangles(self): """Get mesh based on nodes obtained from get_vertex_coordinates. Return array Mx3 array of integers where each row corresponds to a triangle. A triangle is a triplet of indices into point coordinates obtained from get_vertex_coordinates and each index appears only once This provides a mesh where no triangles share nodes (hence the name disconnected triangles) and different nodes may have the same coordinates. This version of the mesh is useful for storing meshes with discontinuities at each node and is e.g. used for storing data in sww files. 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]] """ M = len(self) # Number of triangles K = 3*M       # Total number of unique vertices T = reshape(array(range(K)).astype(Int), (M,3)) return T def get_unique_vertices(self,  indices=None): triangles = self.get_vertices(indices=indices) """FIXME(Ole): This function needs a docstring """ triangles = self.get_triangles(indices=indices) unique_verts = {} for triangle in triangles: def build_vertexlist(self): """Build vertexlist index by vertex ids and for each entry (point id) """Build vertexlist indexed by vertex ids and for each entry (point id) build a list of (triangles, vertex_id) pairs that use the point as vertex. The vertex list will have length N, where N is the number of nodes in the mesh. Preconditions: self.nodes and self.triangles are defined """ vertexlist = [None]*len(self.nodes) vertexlist = [None]*self.number_of_nodes for i in range(self.number_of_triangles):
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

 r3945 """ from Numeric import array, zeros, Float, less, concatenate, NewAxis, argmax, allclose from Numeric import array, zeros, Float, less, concatenate, NewAxis,\ argmax, allclose from anuga.utilities.numerical_tools import ensure_numeric, is_scalar #Method for outputting model results #FIXME: Split up into geometric and numeric stuff. #FIXME: Geometric (X,Y,V) should live in mesh.py #FIXME: STill remember to move XY to mesh # Methods for outputting model results def get_vertex_values(self, xy=True, smooth = None, precision = None, reduction = None): """Return vertex values like an OBJ format smooth=None, precision=None, reduction=None): """Return vertex values like an OBJ format i.e. one value per node. The vertex values are returned as one sequence in the 1D float array A. The connectivity is represented as an integer array, V, of dimension M x 3, where M is the number of volumes. Each row has three indices into the X, Y, A arrays defining the triangle. Mx3, where M is the number of triangles. Each row has three indices defining the triangle and they correspond to elements in the arrays X, Y and A. if smooth is True, vertex values corresponding to one common coordinate set will be smoothed according to the given reduction operator. In this case vertex coordinates will be de-duplicated. de-duplicated corresponding to the original nodes as obtained from the method general_mesh.get_nodes() If no smoothings is required, vertex coordinates and values will be aggregated as a concatenation of values at vertices 0, vertices 1 and vertices 2 vertices 0, vertices 1 and vertices 2. This corresponds to the node coordinates obtained from the method general_mesh.get_vertex_coordinates() if smooth is None: # Take default from domain smooth = self.domain.smooth if precision is None: precision = Float #Create connectivity if smooth == True: if smooth is True: # Ensure continuous vertex values by combining # values at each node using the reduction operator if reduction is None: # Take default from domain reduction = self.domain.reduction V = self.domain.get_vertices() N = len(self.domain.vertexlist) V = self.domain.get_triangles() N = self.domain.number_of_full_nodes # Ignore ghost nodes if any A = zeros(N, precision) #Smoothing loop points = self.domain.get_nodes() # Reduction loop for k in range(N): L = self.domain.vertexlist[k] #Go through all triangle, vertex pairs #contributing to vertex k and register vertex value if L is None: continue #In case there are unused points # Go through all triangle, vertex pairs # contributing to vertex k and register vertex value if L is None: continue # In case there are unused points contributions = [] for volume_id, vertex_id in L: A[k] = reduction(contributions) if xy is True: X = self.domain.get_nodes()[:,0].astype(precision) Y = self.domain.get_nodes()[:,1].astype(precision) return X, Y, A, V else: return A, V else: #Don't smooth #obj machinery moved to general_mesh # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]] # These vert_id's will relate to the verts created below #m = len(self.domain)  #Number of volumes #M = 3*m        #Total number of unique vertices #V = reshape(array(range(M)).astype(Int), (m,3)) V = self.domain.get_triangles(obj=True) #FIXME use get_vertices, when ready A = self.vertex_values.flat #Do vertex coordinates if xy is True: C = self.domain.get_vertex_coordinates() X = C[:,0:6:2].copy() Y = C[:,1:6:2].copy() return X.flat, Y.flat, A, V else: return A, V else: # Allow discontinuous vertex values V = self.domain.get_disconnected_triangles() points = self.domain.get_vertex_coordinates() A = self.vertex_values.flat.astype(precision) # Return if xy is True: X = points[:,0].astype(precision) Y = points[:,1].astype(precision) return X, Y, A, V else: return A, V
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

 r3945 #Create basic mesh points, vertices, boundary = rectangular(1, 3) domain = General_mesh(points, vertices) nodes, triangles, _ = rectangular(1, 3) domain = General_mesh(nodes, triangles) assert allclose(domain.get_vertex_coordinates(unique=True), domain.nodes) assert allclose(domain.get_nodes(), nodes) #assert allclose(domain.get_vertex_coordinates(), ...TODO #assert allclose(domain.get_vertex_coordinates(absolute=True), ...TODO M = domain.number_of_triangles V = domain.get_vertex_coordinates() assert V.shape[0] == 3*M for i in range(M): for j in range(3): k = triangles[i,j]  #Index of vertex j in triangle i assert allclose(V[3*i+j,:], nodes[k]) #Create basic mesh points, vertices, boundary = rectangular(1, 3) domain = General_mesh(points, vertices) nodes, triangles, _ = rectangular(1, 3) domain = General_mesh(nodes, triangles) value = [7] #indexes = [1]  #FIXME (Ole): Should this be used assert  domain.get_vertices() == domain.triangles assert domain.get_vertices([0,4]) == [domain.triangles[0], domain.triangles[4]] assert allclose(domain.get_triangles(), triangles) assert allclose(domain.get_triangles([0,4]), [triangles[0], triangles[4]]) def test_areas(self):
• ## anuga_core/source/anuga/shallow_water/data_manager.py

 r3953 precision=self.precision) x[:] = X.astype(self.precision) y[:] = Y.astype(self.precision) # FIXME (HACK) if len(z) <> len(Z): truncation = self.domain.number_of_full_nodes Z = Z[:truncation] print len(z), len(Z), truncation, len(self.domain.get_nodes()) z[:] = Z.astype(self.precision) volumes[:] = V.astype(volumes.typecode()) x[:] = X y[:] = Y z[:] = Z #FIXME: Backwards compatibility z = fid.variables['z'] z[:] = Z.astype(self.precision) z[:] = Z ################################ volumes[:] = V.astype(volumes.typecode()) #Close precision = self.precision) # HACK truncation = self.domain.number_of_full_nodes # HACK if stage.shape[1] <> len(A): A = A[:truncation] #FIXME: Make this general (see below) outfile.variables['x'][:] = x - xllcorner outfile.variables['y'][:] = y - yllcorner outfile.variables['z'][:] = z             #FIXME HACK outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat. outfile.variables['elevation'][:] = z outfile.variables['time'][:] = times   #Store time relative outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64
• ## anuga_core/source/anuga/shallow_water/shallow_water_domain.py

 r3928 value per vertex using self.reduction. """ # FIXME (Ole): how about using the word continuous vertex values? self.smooth = not flag
• ## anuga_core/source/anuga_parallel/print_stats.py

 r3514 def build_full_flag(domain, ghost_recv_dict): tri_full_flag = ones(len(domain.get_vertices()), Int8) tri_full_flag = ones(len(domain.get_triangles()), Int8) for i in ghost_recv_dict.keys(): for id in ghost_recv_dict[i][0]:
Note: See TracChangeset for help on using the changeset viewer.