Changeset 8713

Timestamp:

Feb 22, 2013, 4:25:41 PM (12 years ago)

Author:

steve

Message:

Looked to speed up check_integrity in neighbour_mesh.

Location:

trunk/anuga_core

Files:

• documentation/user_manual/demos/cairns/project.py (modified) (1 diff)
• documentation/user_manual/demos/cairns/runcairns.py (modified) (2 diffs)
• source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/neighbour_mesh_ext.c (modified) (2 diffs)

trunk/anuga_core/documentation/user_manual/demos/cairns/project.py

@@ -85 +85 @@
 slide_origin = [451871, 8128376]   # Assume to be on continental shelf
 slide_depth = 500.
+
+
+#------------------------------------------------------------------------------
+# Data for Tides
+#------------------------------------------------------------------------------
+tide = 0.0

trunk/anuga_core/documentation/user_manual/demos/cairns/runcairns.py

@@ -68 +68 @@
 # Setup initial conditions
 #------------------------------------------------------------------------------
-tide = 0.0
+tide = project.tide
 domain.set_quantity('stage', tide)
 domain.set_quantity('friction', 0.0)
@@ -146 +146 @@
 if project.scenario == 'fixed_wave':
     # Save every two mins leading up to wave approaching land
-    for t in domain.evolve(yieldstep=120, finaltime=5000): 
+    for t in domain.evolve(yieldstep=2*60, finaltime=5000): 
         print domain.timestepping_statistics()
         print domain.boundary_statistics(tags='ocean_east')    
 
     # Save every 30 secs as wave starts inundating ashore
-    for t in domain.evolve(yieldstep=10, finaltime=10000, 
+    for t in domain.evolve(yieldstep=60*0.5, finaltime=10000, 
                            skip_initial_step=True):
         print domain.timestepping_statistics()

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

@@ -926 +926 @@
         
         # Check number of triangles per node
-        count = [0]*self.number_of_nodes
-        for triangle in self.triangles:
-            for i in triangle:
-                count[i] += 1
-
-        assert num.allclose(count, self.number_of_triangles_per_node)
-
-
-        # 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
+#        count = [0]*self.number_of_nodes
+#        for triangle in self.triangles:
+#            for i in triangle:
+#                count[i] += 1
+
+        count  = num.bincount(self.triangles.flat)
+
+
+        ncount = len(count)
+        #print len(count)
+        #print len(self.number_of_triangles_per_node)
+
+
+        number_of_lone_nodes = self.number_of_nodes - len(self.number_of_triangles_per_node)
+
+
+        assert num.allclose(count, self.number_of_triangles_per_node[:ncount])
+
+
+        from neighbour_mesh_ext import check_integrity_c
+
+
+        #print self.vertex_value_indices.shape
+        #print self.triangles.shape
+        #print self.node_index.shape
+        #print self.number_of_triangles_per_node.shape
+
+        check_integrity_c(self.vertex_value_indices,
+                          self.triangles,
+                          self.node_index,
+                          self.number_of_triangles_per_node)
+
+
+
+#        # 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
 
 

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh_ext.c

@@ -77 +77 @@
 
 
+static PyObject *CheckIntegrity( PyObject *self, PyObject *args )
+{
+        int nt, nt3, tri, n_node, n_node_1;
+        int ok;
+        int current_node;
+        int k, i;
+        int index;
+   
+        //char errorMsg[50];
+
+        double cumsum;
+
+        long *vertex_value_indices;
+        long *triangles;
+        long *node_index;
+        long *number_of_triangles_per_node;
+
+
+        PyObject *pyobj_vertex_value_indices;
+        PyObject *pyobj_triangles;
+        PyObject *pyobj_node_index;
+        PyObject *pyobj_number_of_triangles_per_node;
+
+        //Py_ssize_t pos;
+
+        ok = PyArg_ParseTuple( args, "OOOO",
+                &pyobj_vertex_value_indices,
+                &pyobj_triangles,
+                &pyobj_node_index,
+                &pyobj_number_of_triangles_per_node
+                );
+        if( !ok ){
+                report_python_error(AT, "could not parse input arguments");
+                return NULL;
+        }
+
+        vertex_value_indices = IDATA( pyobj_vertex_value_indices );
+        nt3 = DIMENSHAPE( pyobj_vertex_value_indices, 0 );
+
+
+        //printf(" n3 %d \n",nt3);
+
+        triangles = IDATA( pyobj_triangles );
+        nt = DIMENSHAPE( pyobj_triangles, 0 );
+        tri = DIMENSHAPE( pyobj_triangles, 1 );
+
+        //printf(" nt %d tri %d \n",nt, tri);
+
+        node_index = IDATA( pyobj_node_index );
+        n_node_1 = DIMENSHAPE( pyobj_node_index, 0 );
+
+
+        //printf(" n_node_1 %d  \n",n_node_1);
+
+        number_of_triangles_per_node = IDATA( pyobj_number_of_triangles_per_node );
+        n_node = DIMENSHAPE( pyobj_number_of_triangles_per_node, 0 );
+
+
+        //printf(" n_node %d \n",n_node);
+
+        //----------------------------------------------------------------------
+        // Check that the arrays are consistent
+        //----------------------------------------------------------------------
+        if ( nt3 != 3*nt) {
+                report_python_error(AT, "Mismatch in size of triangles and vertex_value_indices");
+                return NULL;
+        }
+
+        if ( n_node_1 != n_node + 1 ) {
+                report_python_error(AT, "Mismatch in size of node_index and number_of_triangles_per_node");
+                return NULL;
+        }
+
+        if ( tri != 3 ) {
+                report_python_error(AT, "Triangle array should be nt by 3");
+                return NULL;
+        }
+
+        //----------------------------------------------------------------------
+        // Check consistence between triangles, and node_index and
+        // number_of_triangles_per_node
+        //----------------------------------------------------------------------
+        current_node = 0;
+        k = 0; //Track triangles touching on node
+        for (i=0; i<nt3 ; i++) {
+            index = vertex_value_indices[i];
+
+            if (number_of_triangles_per_node[current_node] == 0) {
+                // Node is lone - i.e. not part of the mesh
+                continue;
+            }
+
+            k += 1;
+
+            if ( triangles[index] != current_node) {
+                report_python_error(AT, "Inconsistency between triangles and vertex_values_indices");
+                return NULL;
+            }
+
+            if ( number_of_triangles_per_node[current_node] == k) {
+                // Move on to next node
+                k = 0;
+                current_node += 1;
+            }
+        }
+
+        cumsum = 0.0;
+        for (i=0; i<n_node; i++) {
+            cumsum += number_of_triangles_per_node[i];
+            if ( cumsum != node_index[i+1]) {
+                report_python_error(AT, "Inconsistency between node_index and number_of_triangles_per_node");
+                return NULL;
+            }
+        }
+
+
+/*
+        pos = 0;
+        if( PyDict_Size( pyobj_boundary ) ){
+                // iterate through dictionary entries
+                while( PyDict_Next( pyobj_boundary, &pos, &pyobj_dictKey, &pyobj_dictVal ) ){
+                        volID  = PyLong_AsLong( PyTuple_GetItem( pyobj_dictKey, 0 ) );
+                        edgeID = PyLong_AsLong( PyTuple_GetItem( pyobj_dictKey, 1 ) );
+
+                        if (!(volID<aDimen && edgeID<bDimen)) {
+                            sprintf(errorMsg, "Segment (%d, %d) does not exist", volID, edgeID);
+                            report_python_error(AT, errorMsg);
+                            return NULL;
+                        }
+                }
+        }
+
+        for(volID=0; volID<numTriangle; volID++){
+                for(edgeID=0; edgeID<3; edgeID++){
+                        if( neighbours[volID*3+edgeID] < 0 ){
+                                pyobj_dictKey = PyTuple_New( 2 );
+                                PyTuple_SetItem( pyobj_dictKey, 0, PyInt_FromLong( volID ) );
+                                PyTuple_SetItem( pyobj_dictKey, 1, PyInt_FromLong( edgeID ) );
+
+                                if( !PyDict_Contains( pyobj_boundary, pyobj_dictKey ) )
+                                        PyDict_SetItem( pyobj_boundary, pyobj_dictKey, PyString_FromString( defaultTag ) );
+                        }
+                }
+        }
+*/
+
+        return Py_BuildValue("");
+}
+
+
+
+
+
+
+
 /*
 static PyObject *BoundaryDictionaryConstruct( PyObject *self, PyObject *args )
@@ -114 +269 @@
 static PyMethodDef MF_module_methods[] = {
         {"boundary_dictionary_construct", BoundaryDictionaryConstruct, METH_VARARGS},
+        {"check_integrity_c", CheckIntegrity, METH_VARARGS},
         {NULL, NULL}    //sentinel
 };