Changeset 8501

Timestamp:

Aug 9, 2012, 11:20:54 AM (13 years ago)

Author:

steve

Message:

Moving more mesh code to numpy

Location:

trunk/anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/general_mesh.py (modified) (2 diffs)
• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (1 diff)
• abstract_2d_finite_volumes/test_neighbour_mesh.py (modified) (1 diff)
• operators/erosion_operators.py (modified) (1 diff)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -512 +512 @@
         vertex_coordinates = num.zeros((3*M, 2), num.float)
 
-#        k0 = self.triangles[:,0]
-#        k1 = self.triangles[:,1]
-#        k2 = self.triangles[:,2]
-#
-#        vertex_coordinates[3*k0,:]   = self.nodes[k0,:]
-#        vertex_coordinates[3*k1+1,:] = self.nodes[k1,:]
-#        vertex_coordinates[3*k2+2,:] = self.nodes[k2,:]
-
-        for i in range(M):
-            for j in range(3):
-                k = self.triangles[i,j] # Index of vertex j in triangle i
-                vertex_coordinates[3*i+j,:] = self.nodes[k]
+        k0 = self.triangles[:,0]
+        k1 = self.triangles[:,1]
+        k2 = self.triangles[:,2]
+
+#        I = num.arange(M,dtype=num.int)
+#
+#        V0 = V[0:3*M:3, :]
+#        V1 = V[1:3*M:3, :]
+#        V2 = V[2:3*M:3, :]
+
+        vertex_coordinates[0:3*M:3,:] = self.nodes[k0,:]
+        vertex_coordinates[1:3*M:3,:] = self.nodes[k1,:]
+        vertex_coordinates[2:3*M:3,:] = self.nodes[k2,:]
+
+#        for i in range(M):
+#            for j in range(3):
+#                k = self.triangles[i,j] # Index of vertex j in triangle i
+#                vertex_coordinates[3*i+j,:] = self.nodes[k]
 
         return vertex_coordinates
@@ -758 +764 @@
 
         # Count number of triangles per node
-        number_of_triangles_per_node = num.zeros(self.number_of_nodes,
-                                                 num.int)       #array default#
-        for volume_id, triangle in enumerate(self.get_triangles()):
-            for vertex_id in triangle:
-                number_of_triangles_per_node[vertex_id] += 1
+#        number_of_triangles_per_node = num.zeros(self.number_of_nodes,
+#                                                 num.int)       #array default#
+#        for volume_id, triangle in enumerate(self.get_triangles()):
+#            for vertex_id in triangle:
+#                number_of_triangles_per_node[vertex_id] += 1
+
+        # Need to pad number_of_triangles_per_node in case lone nodes at end of list
+        #number_of_triangles_per_node = num.zeros(self.number_of_nodes, num.int)
+
+        number_of_triangles_per_node = num.bincount(self.triangles.flatten())
+
+        number_of_lone_nodes = self.number_of_nodes - len(number_of_triangles_per_node)
+
+        #print number_of_lone_nodes
+        if number_of_lone_nodes > 0:
+            number_of_triangles_per_node =  \
+               num.append(number_of_triangles_per_node,num.zeros(number_of_lone_nodes,num.int))
+
+        #assert num.allclose(number_of_triangles_per_node_new, number_of_triangles_per_node)
 
         # Allocate space for inverted structure
         number_of_entries = num.sum(number_of_triangles_per_node)
-        vertex_value_indices = num.zeros(number_of_entries, num.int) #array default#
-
-        # Register (triangle, vertex) indices for each node
-        vertexlist = [None] * self.number_of_nodes
-        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
+
+        assert number_of_entries == 3*self.number_of_triangles
+        
+        #vertex_value_indices = num.zeros(number_of_entries, num.int) #array default#
+
+        # Array of vertex_indices (3*vol_id+vertex_id) sorted into contiguous
+        # order around each node. Use with number_of_triangles_per_node to
+        # find vertices associated with a node.
+        # ie There are  number_of_triangles_per_node[i] vertices
+        vertex_value_indices = num.argsort(self.triangles.flatten())
+
+        node_index = num.zeros((self.number_of_nodes)+1, dtype = num.int)
+
         k = 0
-        for vertices in vertexlist:
-            if vertices is not None:
-                for volume_id, vertex_id in vertices:
-                    vertex_value_indices[k] = 3*volume_id + vertex_id
-                    k += 1
-
-        # Save structure
+        node_index[0] = 0
+        for i in range(self.number_of_nodes):
+            node_index[i+1] = node_index[i] + number_of_triangles_per_node[i]
+
+
+        # Save structures
+        self.node_index = node_index
         self.number_of_triangles_per_node = number_of_triangles_per_node
         self.vertex_value_indices = vertex_value_indices

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -793 +793 @@
                        len(self.vertex_value_indices)
 
+        
         # Check number of triangles per node
         count = [0]*self.number_of_nodes

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -1140 +1140 @@
         c = [1.0, 3.0]
         d = [2.0, 4.0]
-
-        points = [a, b, c, d]
-        vertices = [[0,1,2]]
+        e = [4.0, 3.0]
+
+        points = [a, b, d, c, e]
+        vertices = [[0,1,3]]
 
         mesh = Mesh(points, vertices)
         mesh.check_integrity()
         loners = mesh.get_lone_vertices()
-        self.failUnless(loners==[3],
+
+        #print loners
+        self.failUnless(loners==[2,4],
                         'FAILED!')
 

trunk/anuga_core/source/anuga/operators/erosion_operators.py

@@ -375 +375 @@
 
 
-        node_index = num.zeros((self.domain.number_of_nodes)+1, dtype = num.int)
-
-        k = 0
-        node_index[0] = 0
-        for i in range(self.domain.number_of_nodes):
-            node_index[i+1] = node_index[i] + self.domain.number_of_triangles_per_node[i]
-
-        self.node_index = node_index
+        self.node_index = self.domain.node_index
 
         vertex_ids =[]