Changeset 4471

Timestamp:

May 21, 2007, 5:56:14 PM (17 years ago)

Author:

ole

Message:

Rewrote quantity.get_vertex_values in C and introduced new data structure for keeping track of which triangles contribute to each node. This addresses ticket:154 by making storing orders of magnitudes faster.
There is no longer the option of specifying 'reduction' as averaging is now hard wired.

Location:

anuga_core

Files:

• documentation/user_manual/anuga_user_manual.tex (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/general_mesh.py (modified) (2 diffs)
• source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/quantity.py (modified) (5 diffs)
• source/anuga/abstract_2d_finite_volumes/quantity_ext.c (modified) (4 diffs)
• source/anuga/abstract_2d_finite_volumes/test_general_mesh.py (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/test_util.py (modified) (1 diff)
• source/anuga/shallow_water/test_data_manager.py (modified) (1 diff)

anuga_core/documentation/user_manual/anuga_user_manual.tex

@@ -1773 +1773 @@
 
 
-\begin{methoddesc} {set\_store\_vertices\_uniquely}{flag, reduction=None}
+\begin{methoddesc} {set\_store\_vertices\_uniquely}{flag}
 Decide whether vertex values should be stored uniquely as
 computed in the model or whether they should be reduced to one
-value per vertex using self.reduction.
+value per vertex using averaging.
 \end{methoddesc}
 

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -394 +394 @@
         """
 
+        # FIXME (Ole): Refactor this based on algorithm in test and get
+        # rid of the old vertexlist
+        
         vertexlist = [None]*self.number_of_nodes
         for i in range(self.number_of_triangles):
@@ -410 +413 @@
                 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)
+        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
+
+                    k += 1
+
+        assert k == number_of_entries    
         self.vertexlist = vertexlist
+
+        self.number_of_triangles_per_node = number_of_triangles_per_node
+        self.vertex_value_indices = vertex_value_indices
+        self.vertexlist = vertexlist
+
+        #print
+        #print number_of_triangles_per_node
+        #print vertex_value_indices
+        #print vertexlist
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -701 +701 @@
         #
         #See domain.set_boundary
+
+        # Check integrity of vertex_value_indices
+        current_node = 0
+        k = 0 # Track triangles touching on node
+        for index in self.vertex_value_indices:
+
+            if self.number_of_triangles_per_node[current_node] == 0:
+                # Node is lone - i.e. not part of the mesh
+                continue
+            
+            k += 1
+            
+            volume_id = index / 3
+            vertex_id = index % 3
+            
+            msg = 'Triangle %d, vertex %d points to %d. Should have been %d'\
+                  %(volume_id, vertex_id, self.triangles[volume_id, vertex_id], current_node)
+            assert self.triangles[volume_id, vertex_id] == current_node, msg
+                        
+            if self.number_of_triangles_per_node[current_node] == k:
+                # Move on to next node
+                k = 0
+                current_node += 1
+
 
     def get_lone_vertices(self):

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -943 +943 @@
                 
                 # FIXME (Ole): Should we merge this with get_vertex_values
-                # and use the concept of a reduction operator?
                 sum = 0
                 for triangle_id, vertex_id in triangles:
@@ -1033 +1032 @@
                           xy=True,
                           smooth=None,
-                          precision=None,
-                          reduction=None):
+                          precision=None):
         """Return vertex values like an OBJ format i.e. one value per node.
 
@@ -1046 +1044 @@
 
         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
+        coordinate set will be smoothed by taking the average of vertex values for each node.
+        In this case vertex coordinates will be
         de-duplicated corresponding to the original nodes as obtained from
         the method general_mesh.get_nodes()
@@ -1078 +1076 @@
 
         if smooth is True:
-            # Ensure continuous vertex values by combining
-            # values at each node using the reduction operator
+            # Ensure continuous vertex values by averaging
+            # values at each node
             
-            if reduction is None:
-                # Take default from domain                
-                reduction = self.domain.reduction
-
             V = self.domain.get_triangles()
             N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
-            A = zeros(N, precision)
+            A = zeros(N, Float)
             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
-                contributions = []
-                for volume_id, vertex_id in L:
+            
+            if 1:
+                # Fast C version
+                average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
+                                      ensure_numeric(self.domain.number_of_triangles_per_node),
+                                      ensure_numeric(self.vertex_values),
+                                      A)
+                A = A.astype(precision)
+            else:    
+
+                # Slow Python version
+                
+                current_node = 0
+                k = 0 # Track triangles touching on node
+                total = 0.0
+                for index in self.domain.vertex_value_indices:
+                    k += 1
+                    
+                    volume_id = index / 3
+                    vertex_id = index % 3
+                 
+                    #assert V[volume_id, vertex_id] == current_node
+                
                     v = self.vertex_values[volume_id, vertex_id]
-                    contributions.append(v)
-
-                A[k] = reduction(contributions)
+                    total += v
+
+                    #print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
+                    if self.domain.number_of_triangles_per_node[current_node] == k:
+                        A[current_node] = total/k
+                
+                    
+                        # Move on to next node
+                        total = 0.0
+                        k = 0
+                        current_node += 1
+
+
 
         else:
@@ -1450 +1467 @@
 
     from quantity_ext import compute_gradients, limit,\
-    extrapolate_second_order, interpolate_from_vertices_to_edges, update
+    extrapolate_second_order, interpolate_from_vertices_to_edges, update,\
+    average_vertex_values

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -218 +218 @@
 
 
+int _average_vertex_values(int N,
+                           long* vertex_value_indices,
+                           long* number_of_triangles_per_node,
+                           double* vertex_values, 
+                           double* A) {
+  //Average vertex values to obtain one value per node
+
+  int i, index; 
+  int k = 0; //Track triangles touching each node
+  int current_node = 0;
+  double total = 0.0;
+
+  for (i=0; i<N; i++) {
+    index = vertex_value_indices[i];
+    k += 1;
+            
+    //volume_id = index / 3
+    //vertex_id = index % 3
+    //total += self.vertex_values[volume_id, vertex_id]
+    total += vertex_values[index];
+      
+    //printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);
+    if (number_of_triangles_per_node[current_node] == k) {
+      A[current_node] = total/k;
+                
+      // Move on to next node
+      total = 0.0;
+      k = 0;
+      current_node += 1;
+    }
+  }
+                           
+  return 0;
+}
+
+
 /////////////////////////////////////////////////
 // Gateways to Python
@@ -283 +319 @@
 
         err = _interpolate(N,
-                     (double*) vertex_values -> data,
-                     (double*) edge_values -> data);
+                           (double*) vertex_values -> data,
+                           (double*) edge_values -> data);
 
         if (err != 0) {
@@ -298 +334 @@
         return Py_BuildValue("");
 }
+
+
+PyObject *average_vertex_values(PyObject *self, PyObject *args) {
+
+        PyArrayObject 
+          *vertex_value_indices, 
+          *number_of_triangles_per_node,
+          *vertex_values,
+          *A;
+        
+
+        int N, err;
+
+        // Convert Python arguments to C
+        if (!PyArg_ParseTuple(args, "OOOO",
+                              &vertex_value_indices, 
+                              &number_of_triangles_per_node,
+                              &vertex_values, 
+                              &A)) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: average_vertex_values could not parse input");
+          return NULL;
+        }
+        
+        N = vertex_value_indices -> dimensions[0];
+        // printf("Got parameters, N=%d\n", N);
+        err = _average_vertex_values(N,
+                                     (long*) vertex_value_indices -> data,
+                                     (long*) number_of_triangles_per_node -> data,
+                                     (double*) vertex_values -> data, 
+                                     (double*) A -> data);
+
+        if (err != 0) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "average_vertex_values could not be computed");
+          return NULL;
+        }
+
+        return Py_BuildValue("");
+}
+
 
 
@@ -567 +644 @@
                 interpolate_from_vertices_to_edges,
                 METH_VARARGS, "Print out"},
+        {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"},            
         {NULL, NULL, 0, NULL}   // sentinel
 };

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -61 +61 @@
                         [triangles[0], triangles[4]])
         
+
+    def test_vertex_value_indices(self):
+        """Check that structures are correct.
+        """
+        from mesh_factory import rectangular
+        from Numeric import zeros, Float
+
+        #Create basic mesh
+        nodes, triangles, _ = rectangular(3, 6)
+        domain = General_mesh(nodes, triangles)
+
+        assert sum(domain.number_of_triangles_per_node) ==\
+               len(domain.vertex_value_indices)
+
+        # Check number of triangles per node
+        count = [0]*domain.number_of_nodes
+        for triangle in triangles:
+            for i in triangle:
+                count[i] += 1
+
+        assert allclose(count, domain.number_of_triangles_per_node)
+
+        #print nodes
+        #print triangles
+        #print domain.number_of_triangles_per_node
+        #print domain.vertex_value_indices        
+
+
+        # Check indices
+        current_node = 0
+        k = 0 # Track triangles touching on node
+        for index in domain.vertex_value_indices:
+            k += 1
+            
+            triangle = index / 3
+            vertex = index % 3
+
+            assert triangles[triangle, vertex] == current_node
+
+            if domain.number_of_triangles_per_node[current_node] == k:
+                # Move on to next node
+                k = 0
+                current_node += 1
+                
+
+
+
+
     def test_areas(self):
         from mesh_factory import rectangular

anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

@@ -1331 +1331 @@
     suite = unittest.makeSuite(Test_Util,'test')
 #    suite = unittest.makeSuite(Test_Util,'test_store_parameters')
-    runner = unittest.TextTestRunner()
+    runner = unittest.TextTestRunner(verbosity=0)
     runner.run(suite)
 

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -393 +393 @@
         os.remove(sww.filename)
 
-    def test_sww_variable3(self):
+    def no_test_sww_variable3(self):
         """Test that sww information can be written correctly
         multiple timesteps using a different reduction operator (min)
         """
 
+        # Different reduction in sww files has been made obsolete.
+        
         import time, os
         from Numeric import array, zeros, allclose, Float, concatenate