Changeset 4478

Timestamp:

May 23, 2007, 12:13:06 PM (17 years ago)

Author:

ole

Message:

Retired old structure vertexlist completely and replaced all references to it with the new inverted triangle structured defined in general_mesh.py. This closes ticket:154

All tests and autovalidation passed.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/general_mesh.py (modified) (3 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (4 diffs)
• abstract_2d_finite_volumes/test_general_mesh.py (modified) (3 diffs)
• fit_interpolate/search_functions.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -210 +210 @@
 
         
-        # Build vertex list
-        if verbose: print 'Building vertex list'         
-        self.build_vertexlist()
+        # Build structure listing which trianglse belong to which nodet.
+        if verbose: print 'Building inverted triangle structure'         
+        self.build_inverted_triangle_structure()
 
             
@@ -379 +379 @@
 
 
-    def build_vertexlist(self):
-        """Build information about which triangles belong to each node
+    def get_triangles_and_vertices_per_node(self, node=None):
+        """Get triangles associated with given node.
+
+        Return list of triangle_ids, vertex_ids for specified node.
+        If node in None or absent, this information will be returned
+        for all (full) nodes.
+        """
+
+        triangle_list = []
+        if node is not None:
+            # Get index for this node
+            first = sum(self.number_of_triangles_per_node[:node])
+            
+            # Get number of triangles for this node
+            count = self.number_of_triangles_per_node[node]
+
+            for i in range(count):
+                index = self.vertex_value_indices[first+i]
+
+                volume_id = index / 3
+                vertex_id = index % 3
+
+                triangle_list.append( (volume_id, vertex_id) )
+        else:
+            # Get info for all nodes recursively.
+            # If need be, we can speed this up by
+            # 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):
+        """Build structure listing triangles belonging to each node
 
         Two arrays are created and store as mesh attributes
@@ -441 +475 @@
         """
 
-        # FIXME (Ole): Refactor this based on algorithm in test and get
-        # rid of the old vertexlist
-        
+        # Count number of triangles per node
+        number_of_triangles_per_node = zeros(self.number_of_nodes)
+        for volume_id, triangle in enumerate(self.triangles):
+            for vertex_id in triangle:
+                number_of_triangles_per_node[vertex_id] += 1
+
+        # Allocate space for inverted structure
+        number_of_entries = sum(number_of_triangles_per_node)
+        vertex_value_indices = zeros(number_of_entries)
+
+        # Register (triangle, vertex) indices for each node
         vertexlist = [None]*self.number_of_nodes
-        for i in range(self.number_of_triangles):
-
-            a = self.triangles[i, 0]
-            b = self.triangles[i, 1]
-            c = self.triangles[i, 2]
-
-            #Register the vertices v as lists of
-            #(triangle_id, vertex_id) tuples associated with them
-            #This is used for averaging multiple vertex values.
-            for vertex_id, v in enumerate([a,b,c]):
-                if vertexlist[v] is None:
-                    vertexlist[v] = []
-
-                vertexlist[v].append( (i, vertex_id) )
-
-
-        # Register number of triangles touching each nodes
-        number_of_triangles_per_node = zeros(self.number_of_nodes)        
-        number_of_entries = 0        
-        for i, vertices in enumerate(vertexlist):
-            try:
-                v = len(vertices)
-            except:
-                v = 0 # Lone node - e.g. not part of the mesh
-                
-            number_of_triangles_per_node[i] = v
-            number_of_entries += v
-                
-
-        # Build vertex value index array
-        vertex_value_indices = zeros(number_of_entries)
+        for volume_id in range(self.number_of_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]):
+                if vertexlist[node_id] is None:
+                    vertexlist[node_id] = []
+        
+                vertexlist[node_id].append( (volume_id, vertex_id) )
+
+
+        # Build inverted triangle index array
         k = 0
         for vertices in vertexlist:
             if vertices is not None:
-                for i, v_id in vertices:
-                    vertex_value_indices[k] = 3*i + v_id
-
+                for volume_id, vertex_id in vertices:
+                    vertex_value_indices[k] = 3*volume_id + vertex_id
+                                        
                     k += 1
 
-        assert k == number_of_entries    
-        self.vertexlist = vertexlist
-
+        # Save structure
         self.number_of_triangles_per_node = number_of_triangles_per_node
         self.vertex_value_indices = vertex_value_indices
-        
-
-        #print
-        #print number_of_triangles_per_node
-        #print vertex_value_indices
-        #print vertexlist
-
+
+
+        
 
     def get_extent(self, absolute=False):

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -407 +407 @@
         elif location == 'unique vertices':
             if indices is None:
-                self.edge_values[:] = X
+                self.edge_values[:] = X  #FIXME (Ole): Shouldn't this be vertex_values?
             else:
 
                 #Go through list of unique vertices
                 for unique_vert_id in indices:
-                    triangles = self.domain.vertexlist[unique_vert_id]
-
+
+                    triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
+                    
                     #In case there are unused points
-                    if triangles is None: continue
-
+                    if len(triangles) == 0:
+                        continue
+                    
                     #Go through all triangle, vertex pairs
                     #and set corresponding vertex value
@@ -932 +934 @@
             #Go through list of unique vertices
             for unique_vert_id in indices:
-                triangles = self.domain.vertexlist[unique_vert_id]
-
+                triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
+                    
                 #In case there are unused points
-                if triangles is None:
+                if len(triangles) == 0:
                     msg = 'Unique vertex not associated with triangles'
                     raise msg
@@ -958 +960 @@
         """Set vertex values for all unique vertices based on input array A
         which has one entry per unique vertex, i.e.
-        one value for each row in array self.domain.coordinates or
-        one value for each row in vertexlist.
+        one value for each row in array self.domain.nodes.
 
         indices is the list of vertex_id's that will be set.
@@ -982 +983 @@
 
         #Go through list of unique vertices
+        
         for i_index, unique_vert_id in enumerate(vertex_list):
-            triangles = self.domain.vertexlist[unique_vert_id]
-
-            if triangles is None: continue #In case there are unused points
+
+
+            triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
+                    
+            #In case there are unused points
+            if len(triangles) == 0: continue
 
             #Go through all triangle, vertex pairs

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -96 +96 @@
                     count[i] += 1
 
+            #print count
+            #
             assert allclose(count, domain.number_of_triangles_per_node)
-
+            
             # Check indices
             current_node = 0
@@ -115 +117 @@
                 
 
-
+    def test_get_triangles_and_vertices_per_node(self):
+        """test_get_triangles_and_vertices_per_node -
+
+        Test that tuples of triangle, vertex can be extracted
+        from inverted triangles structure
+
+        """
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        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]])
+
+        domain = General_mesh(nodes, triangles)
+
+        # One node
+        L = domain.get_triangles_and_vertices_per_node(node=2)
+
+        assert allclose(L[0], [0, 2])
+        assert allclose(L[1], [1, 1])
+        assert allclose(L[2], [2, 1])
+
+        # All nodes
+        ALL = domain.get_triangles_and_vertices_per_node()
+
+        assert len(ALL) == 6
+        for i, Lref in enumerate(ALL):
+            L = domain.get_triangles_and_vertices_per_node(node=i)
+            assert allclose(L, Lref)
+            
+
+        
+
+        
 
 
@@ -156 +197 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_General_Mesh,'test')
+    suite = unittest.makeSuite(Test_General_Mesh,'test_vertex_value_indices')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/fit_interpolate/search_functions.py

@@ -64 +64 @@
         #FIXME (DSG-DSG): this catches verts with no triangle.
         #Currently pmesh is producing these.
-        #this should be stopped, 
-        if mesh.vertexlist[v] is None:
+        #this should be stopped,
+
+        if mesh.number_of_triangles_per_node[v] == 0:
             continue
+        
         #for each triangle id (k) which has v as a vertex
-        for k, _ in mesh.vertexlist[v]:
+
+        vertexlist = mesh.get_triangles_and_vertices_per_node(node=v)
+        for k, _ in vertexlist:
             #Get the three vertex_points of candidate triangle
             xi0 = mesh.get_vertex_coordinate(k, 0)