Changeset 8563

Timestamp:

Sep 9, 2012, 8:30:25 AM (13 years ago)

Author:

steve

Message:

looking at pmesh2domain and trying to speed up calc_sides

Location:

trunk/anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/neighbour_mesh_ext.c (modified) (3 diffs)
• abstract_2d_finite_volumes/pmesh2domain.py (modified) (1 diff)
• abstract_2d_finite_volumes/test_neighbour_mesh.py (modified) (1 diff)
• file_conversion/sww2dem.py (modified) (1 diff)

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

@@ -3 +3 @@
 #include <stdio.h>
 
+#include "util_ext.h"
+
 #define DDATA(p) ((double*)(((PyArrayObject *)p)->data))
 #define IDATA(p) ((long*)(((PyArrayObject *)p)->data))
+#define DIMENSHAPE(p, d) (((PyArrayObject *)p)->dimensions[d])
+
+
+#ifndef NDEBUG
+#define ASSERT(condition, message) \
+        do {\
+                if(!(condition)) { \
+                        printf("Assertion `%s` failed in %s line %d: %s\n", \
+                                          #condition, __FILE__, __LINE__, message); \
+                        exit(EXIT_FAILURE); \
+                } \
+        } while (0)
+#else
+#define ASSERT(condition, message) do { } while (0)
+#endif
+
+static PyObject *BoundaryDictionaryConstruct( PyObject *self, PyObject *args )
+{
+        int aDimen, bDimen, ok, numTriangle, volID, edgeID;
+        char *defaultTag;
+        char errorMsg[50];
+        long *neighbours;
+        PyObject *pyobj_dictKey, *pyobj_dictVal;
+        PyObject *pyobj_neighbours;
+        PyObject *pyobj_boundary;
+        Py_ssize_t pos;
+
+        ok = PyArg_ParseTuple( args, "isOO", &numTriangle, &defaultTag, &pyobj_neighbours, &pyobj_boundary );
+        if( !ok ){
+                report_python_error(AT, "could not parse input arguments");
+                return NULL;
+        }
+
+        neighbours = IDATA( pyobj_neighbours );
+        aDimen = DIMENSHAPE( pyobj_neighbours, 0 );
+        bDimen = DIMENSHAPE( pyobj_neighbours, 1 );
+
+        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("O", pyobj_boundary);
+}
 
 
 
-
+/*
 static PyObject *BoundaryDictionaryConstruct( PyObject *self, PyObject *args )
 {
-        int i, j, initSize, ok, numTriangle;
+        int i, j, ok, numTriangle;
         char *defaultTag;
         long *neighbours;
@@ -20 +89 @@
         ok = PyArg_ParseTuple( args, "isOO", &numTriangle, &defaultTag, &pyobj_neighbours, &pyobj_boundary );
         if( !ok ){
-                fprintf( stderr, "BuildBoundaryFromScratch func: argument parsing error\n" );
-                exit(1);
+                report_python_error(AT, "could not parse input arguments");
+                return NULL;
         }
 
@@ -41 +110 @@
         return Py_BuildValue("O", pyobj_boundary);
 }
+*/
 
 static PyMethodDef MF_module_methods[] = {

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

@@ -199 +199 @@
     return sides
 
+def calc_sides_zip(triangles):
+    '''Build dictionary mapping from sides (2-tuple of points)
+    to left hand side neighbouring triangle
+    '''
+
+    sides = {}
+
+
+    triangles = num.array(triangles,num.int)
+
+
+    a = triangles[:,0]
+    b = triangles[:,1]
+    c = triangles[:,2]
+
+    id = num.arange(len(triangles))
+    face = num.ones(len(triangles),num.int)
+
+    sides.update(dict(zip(zip(a,b),zip(id,2*face))))
+    sides.update(dict(zip(zip(b,c),zip(id,0*face))))
+    sides.update(dict(zip(zip(c,a),zip(id,1*face))))
+
+    return sides
+
+def calc_sides_c(triangles):
+    '''Build dictionary mapping from sides (2-tuple of points)
+    to left hand side neighbouring triangle
+    '''
+
+    sides = {}
+
+
+    triangles = num.array(triangles,num.int)
+    ntriangles = len(triangles)
+
+#    print 'calc_sides'
+#    print type(triangles)
+
+    print ntriangles
+
+    from pmesh2domain_ext import sides_dictionary_construct
+    sides = sides_dictionary_construct(triangles, sides)
+
+
+#    old_sides = {}
+#
+#    for id, triangle in enumerate(triangles):
+#        a = int(triangle[0])
+#        b = int(triangle[1])
+#        c = int(triangle[2])
+#
+#        old_sides[a,b] = (id, 2) #(id, face)
+#        old_sides[b,c] = (id, 0) #(id, face)
+#        old_sides[c,a] = (id, 1) #(id, face)
+
+
+    #print sides
+    #print old_sides
+
+
+    return sides

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

@@ -693 +693 @@
             mesh = Mesh(points, vertices, {(5,0): 'x'})
         except AssertionError:
+            pass
+        except RuntimeError:
             pass
         else:

trunk/anuga_core/source/anuga/file_conversion/sww2dem.py

@@ -329 +329 @@
 
 
-
-        #
-
         interp = Interpolate(vertex_points, volumes, verbose = verbose)