Index: anuga_core/source/anuga/advection/advection.py
===================================================================
--- anuga_core/source/anuga/advection/advection.py	(revision 4975)
+++ anuga_core/source/anuga/advection/advection.py	(revision 4978)
@@ -66,9 +66,9 @@
                                 numproc=numproc)
 
-
+        import Numeric
         if velocity is None:
-            self.velocity = [1,0]
+            self.velocity = Numeric.array([1,0],'d')
         else:
-            self.velocity = velocity
+            self.velocity = Numeric.array(velocity,'d')
 
         #Only first is implemented for advection
@@ -108,13 +108,77 @@
         return flux, max_speed
 
+
+
     def compute_fluxes(self):
-
-        
-        self.compute_fluxes_ext()
-##         try:
-##             self.compute_fluxes_ext()
-##         except:
-##             self.compute_fluxes_python()
- 
+        """Compute all fluxes and the timestep suitable for all volumes
+        in domain.
+
+        Compute total flux for each conserved quantity using "flux_function"
+
+        Fluxes across each edge are scaled by edgelengths and summed up
+        Resulting flux is then scaled by area and stored in
+        domain.explicit_update
+
+        The maximal allowable speed computed by the flux_function
+        for each volume
+        is converted to a timestep that must not be exceeded. The minimum of
+        those is computed as the next overall timestep.
+
+        Post conditions:
+        domain.explicit_update is reset to computed flux values
+        domain.timestep is set to the largest step satisfying all volumes.
+        """
+
+        import sys
+        from Numeric import zeros, Float
+        from anuga.config import max_timestep
+
+
+        huge_timestep = float(sys.maxint)
+        Stage = self.quantities['stage']
+
+        """
+        print "======================================"
+        print "BEFORE compute_fluxes"
+        print "stage_update",Stage.explicit_update
+        print "stage_edge",Stage.edge_values
+        print "stage_bdry",Stage.boundary_values
+        print "neighbours",self.neighbours
+        print "neighbour_edges",self.neighbour_edges
+        print "normals",self.normals
+        print "areas",self.areas
+        print "radii",self.radii
+        print "edgelengths",self.edgelengths
+        print "tri_full_flag",self.tri_full_flag
+        print "huge_timestep",huge_timestep
+        print "max_timestep",max_timestep
+        print "velocity",self.velocity
+        """
+
+        import advection_ext		
+        self.timestep = advection_ext.compute_fluxes(self, Stage, huge_timestep, max_timestep)
+
+
+
+    def evolve(self,
+               yieldstep = None,
+               finaltime = None,
+               duration = None,
+               skip_initial_step = False):
+
+        """Specialisation of basic evolve method from parent class
+        """
+
+        #Call basic machinery from parent class
+        for t in Generic_domain.evolve(self,
+                                       yieldstep=yieldstep,
+                                       finaltime=finaltime,
+                                       duration=duration,
+                                       skip_initial_step=skip_initial_step):
+
+            #Pass control on to outer loop for more specific actions
+            yield(t)
+
+
 
 
@@ -207,104 +271,2 @@
         self.timestep = timestep
 
-
-    def compute_fluxes_ext(self):
-        """Compute all fluxes and the timestep suitable for all volumes
-        in domain.
-
-        Compute total flux for each conserved quantity using "flux_function"
-
-        Fluxes across each edge are scaled by edgelengths and summed up
-        Resulting flux is then scaled by area and stored in
-        domain.explicit_update
-
-        The maximal allowable speed computed by the flux_function
-        for each volume
-        is converted to a timestep that must not be exceeded. The minimum of
-        those is computed as the next overall timestep.
-
-        Post conditions:
-        domain.explicit_update is reset to computed flux values
-        domain.timestep is set to the largest step satisfying all volumes.
-        """
-
-        import sys
-        from Numeric import zeros, Float
-        from anuga.config import max_timestep
-
-        ntri = len(self)
-
-        timestep = max_timestep
-
-        #Shortcuts
-        Stage = self.quantities['stage']
-
-        #Arrays
-        neighbours      = self.neighbours
-        neighbour_edges = self.neighbour_edges
-        normals         = self.normals
-        areas           = self.areas
-        radii           = self.radii
-        edgelengths     = self.edgelengths
-        tri_full_flag   = self.tri_full_flag
-
-        stage_edge      = Stage.edge_values
-        stage_bdry      = Stage.boundary_values
-        stage_update    = Stage.explicit_update
-
-        huge_timestep = float(sys.maxint)
-
-        v = self.velocity
-
-        nbdry = len(stage_bdry)
-
-        from advection_ext import compute_fluxes
-
-        print "stage_update",stage_update
-        print "stage_edge",stage_edge
-        print "stage_bdry",stage_bdry
-        print "neighbours",neighbours
-        print "neighbour_edges",neighbour_edges
-        print "normals",normals
-        print "areas",areas
-        print "radii",radii
-        print "edgelengths",edgelengths
-        print "tri_full_flag",tri_full_flag
-        print "huge_timestep",huge_timestep
-        print "max_timestep",max_timestep
-        print "v",v
-        print "ntri",ntri
-        print "nbdry",nbdry
-
-
-		
-        timestep = compute_fluxes(stage_update,stage_edge,stage_bdry,
-                                  neighbours,neighbour_edges,normals,
-                                  areas,radii,edgelengths,
-                                  tri_full_flag,
-                                  huge_timestep,max_timestep,v,ntri,nbdry)
-
-        print "stage_update",stage_update        
-        
-        #print 'timestep out2 =',timestep
-
-        self.timestep = timestep
-
-
-    def evolve(self,
-               yieldstep = None,
-               finaltime = None,
-               duration = None,
-               skip_initial_step = False):
-
-        """Specialisation of basic evolve method from parent class
-        """
-
-        #Call basic machinery from parent class
-        for t in Generic_domain.evolve(self,
-                                       yieldstep=yieldstep,
-                                       finaltime=finaltime,
-                                       duration=duration,
-                                       skip_initial_step=skip_initial_step):
-
-            #Pass control on to outer loop for more specific actions
-            yield(t)
Index: anuga_core/source/anuga/advection/advection_ext.c
===================================================================
--- anuga_core/source/anuga/advection/advection_ext.c	(revision 4975)
+++ anuga_core/source/anuga/advection/advection_ext.c	(revision 4978)
@@ -1,60 +1,41 @@
+// Python - C extension module for advection.py
+//
+// To compile (Python2.X):
+// use python ../utilities/compile.py
+//
+// See the module advection.py
+//
+//
+// Steve Roberts, ANU 2008
+
+
+#include "Python.h"
+#include "Numeric/arrayobject.h"
 #include "math.h"
 #include "stdio.h"
 
-void  print_double_array(char* name, double* array, int n, int m){
-
-    int k,i,km;
-
-    printf("%s = [",name);
-    for (k=0; k<n; k++){
-	km = m*k;
-	printf("[");
-	for (i=0; i<m ; i++){
-	    printf("%g ",array[km+i]);
-	}
-	if (k==(n-1))
-	    printf("]");
-	else
-	    printf("]\n");
-    }
-    printf("]\n");
-}
-
-void  print_int_array(char* name, int* array, int n, int m){
-
-    int k,i,km;
-
-    printf("%s = [",name);
-    for (k=0; k<n; k++){
-	km = m*k;
-	printf("[");
-	for (i=0; i<m ; i++){
-	    printf("%i ",array[km+i]);
-	}
-	if (k==(n-1))
-	    printf("]");
-	else
-	    printf("]\n");
-    }
-    printf("]\n");
-}
-
-
-double compute_fluxes(
-		    double* stage_update,
-		    double* stage_edge,
-		    double* stage_bdry,
-                    int* neighbours,
-		    int* neighbour_edges,
-		    double* normals,
-                    double* areas,
-		    double* radii,
-		    double* edgelengths,
-		    int*   tri_full_flag,
+// Shared code snippets
+#include "util_ext.h"
+
+
+
+double _compute_fluxes(
+		    double* quantity_update,
+		    double* quantity_edge,
+		    double* quantity_bdry,
+                    long*   domain_neighbours,
+		    long*   domain_neighbour_edges,
+		    double* domain_normals,
+                    double* domain_areas,
+		    double* domain_radii,
+		    double* domain_edgelengths,
+		    long*   domain_tri_full_flag,
+		    double* domain_velocity,
                     double  huge_timestep,
                     double  max_timestep,
-		    double* v,
 		    int ntri,
 		    int nbdry){
+
+ 
         //Loop
 
@@ -72,26 +53,4 @@
 	timestep = max_timestep;
 
-
-	printf("======================================================\n");
-	printf("INSIDE compute_fluxes\n");
-
-
-        print_double_array( "stage_update",stage_update, ntri, 1);
-        print_double_array( "stage_edge",stage_edge, ntri, 3);
-        print_double_array( "stage_bdry",stage_bdry, nbdry, 1);
-        print_int_array( "neighbours",neighbours, ntri, 3);
-        print_int_array( "neighbour_edges",neighbour_edges, ntri, 3);
-        print_double_array( "normals",normals, ntri, 6);
-        print_double_array( "areas",areas, ntri, 1);
-        print_double_array( "radii",radii, ntri, 1);
-        print_double_array( "edgelengths",edgelengths, ntri, 3);
-        print_int_array( "tri_full_flag",tri_full_flag, ntri, 1);
-	printf("huge_timestep = %g\n",huge_timestep);
-	printf("max_timestep = %g\n",max_timestep);
-        print_double_array( "v",v, 2, 1);
-	printf("ntri = %i\n",ntri);
-	printf("nbdry = %i\n",nbdry);
-
-
         for (k=0; k<ntri; k++){
             optimal_timestep = huge_timestep;
@@ -101,16 +60,16 @@
 	        k_i = k3+i;
                 //Quantities inside volume facing neighbour i
-                ql = stage_edge[k_i];
+                ql = quantity_edge[k_i];
 
 
                 //Quantities at neighbour on nearest face
-                n = neighbours[k_i];
+                n = domain_neighbours[k_i];
                 if (n < 0) {
                     m = -n-1; //Convert neg flag to index
-                    qr = stage_bdry[m];
+                    qr = quantity_bdry[m];
                 } else {
-                    m = neighbour_edges[k_i];
+                    m = domain_neighbour_edges[k_i];
 		    n_m = 3*n+m;
-                    qr = stage_edge[n_m];
+                    qr = quantity_edge[n_m];
                 }
 
@@ -119,10 +78,10 @@
                 for (j=0; j<2; j++){
 		    k_i_j = 6*k+2*i+j;
-                    normal[j] = normals[k_i_j];
+                    normal[j] = domain_normals[k_i_j];
                 }
 
 
                 //Flux computation using provided function
-                normal_velocity = v[0]*normal[0] + v[1]*normal[1];
+                normal_velocity = domain_velocity[0]*normal[0] + domain_velocity[1]*normal[1];
 
                 if (normal_velocity < 0) {
@@ -133,10 +92,11 @@
 
                 max_speed = fabs(normal_velocity);
-                flux = flux - edgeflux * edgelengths[k_i];
+                flux = flux - edgeflux * domain_edgelengths[k_i];
 
                 //Update optimal_timestep
-                if (tri_full_flag[k] == 1) {
+                if (domain_tri_full_flag[k] == 1) {
                     if (max_speed != 0.0) {
-                        optimal_timestep = (optimal_timestep>radii[k]/max_speed) ? radii[k]/max_speed : optimal_timestep;
+                        optimal_timestep = (optimal_timestep>domain_radii[k]/max_speed) ? 
+			  domain_radii[k]/max_speed : optimal_timestep;
                     }
                 }
@@ -146,23 +106,138 @@
             //Normalise by area and store for when all conserved
             //quantities get updated
-            stage_update[k] = flux/areas[k];
+            quantity_update[k] = flux/domain_areas[k];
 
             timestep = (timestep>optimal_timestep) ? optimal_timestep : timestep;
 
         }
-
-	//printf("timestep out = %g\n",timestep);
-
-
-
-	printf("INSIDE compute_fluxes, end \n");
-
-        print_double_array( "stage_update",stage_update, ntri, 1);
-
-	printf("FINISHED compute_fluxes\n");
-	printf("======================================================\n");
-
 
 	return timestep;
 }
 
+
+//-----------------------------------------------------
+// Python method Wrappers 
+//-----------------------------------------------------
+
+
+
+PyObject *compute_fluxes(PyObject *self, PyObject *args) {
+  /*Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for the conserved quantity stage.
+
+    The maximal allowable speed computed for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Python call:
+    timestep = advection_ext.compute_fluxes(domain,quantity,huge_timestep,max_timestep)
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+
+  */
+
+  PyObject *domain, *quantity;
+  PyArrayObject 
+    * quantity_update,
+    * quantity_edge,
+    * quantity_bdry,
+    * domain_neighbours,
+    * domain_neighbour_edges,
+    * domain_normals,
+    * domain_areas,
+    * domain_radii,
+    * domain_edgelengths,
+    * domain_tri_full_flag,
+    * domain_velocity;
+                  
+
+  // Local variables
+  int ntri, nbdry;
+  double huge_timestep, max_timestep;
+  double timestep;
+
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOdd",
+			&domain,
+			&quantity,
+			&huge_timestep,
+			&max_timestep)) {
+    PyErr_SetString(PyExc_RuntimeError, "advection_ext.c: compute_fluxes could not parse input");
+    return NULL;
+  }
+
+  quantity_edge          = get_consecutive_array(quantity, "edge_values");
+  quantity_bdry          = get_consecutive_array(quantity, "boundary_values");
+  quantity_update        = get_consecutive_array(quantity, "explicit_update");
+  domain_neighbours      = get_consecutive_array(domain,   "neighbours");
+  domain_neighbour_edges = get_consecutive_array(domain,   "neighbour_edges");
+  domain_normals         = get_consecutive_array(domain,   "normals");
+  domain_areas           = get_consecutive_array(domain,   "areas");
+  domain_radii           = get_consecutive_array(domain,   "radii");
+  domain_edgelengths     = get_consecutive_array(domain,   "edgelengths");
+  domain_tri_full_flag   = get_consecutive_array(domain,   "tri_full_flag");
+  domain_velocity        = get_consecutive_array(domain,   "velocity");  
+    
+  ntri  = quantity_edge -> dimensions[0];
+  nbdry = quantity_bdry -> dimensions[0]; 
+
+  timestep = _compute_fluxes((double*) quantity_update -> data,
+			     (double*) quantity_edge -> data,
+			     (double*) quantity_bdry -> data,
+			     (long*)   domain_neighbours -> data,
+			     (long*)   domain_neighbour_edges -> data,
+			     (double*) domain_normals -> data,
+			     (double*) domain_areas ->data,
+			     (double*) domain_radii -> data,
+			     (double*) domain_edgelengths -> data,
+			     (long*)   domain_tri_full_flag -> data,
+			     (double*) domain_velocity -> data,
+			     huge_timestep,
+			     max_timestep,
+			     ntri,
+			     nbdry);
+
+  // Release and return
+  Py_DECREF(quantity_update);
+  Py_DECREF(quantity_edge);
+  Py_DECREF(quantity_bdry);
+  Py_DECREF(domain_neighbours);
+  Py_DECREF(domain_neighbour_edges);
+  Py_DECREF(domain_normals);
+  Py_DECREF(domain_areas);
+  Py_DECREF(domain_radii);
+  Py_DECREF(domain_edgelengths);
+  Py_DECREF(domain_tri_full_flag);
+  Py_DECREF(domain_velocity);
+
+
+  return Py_BuildValue("d", timestep);
+}
+
+
+
+//-------------------------------
+// Method table for python module
+//-------------------------------
+static struct PyMethodDef MethodTable[] = {
+  /* The cast of the function is necessary since PyCFunction values
+   * only take two PyObject* parameters, and rotate() takes
+   * three.
+   */
+
+  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initadvection_ext(void){
+  Py_InitModule("advection_ext", MethodTable);
+  import_array(); // Necessary for handling of NumPY structures
+}
Index: anuga_core/source/anuga/advection/advection_ext.pyf
===================================================================
--- anuga_core/source/anuga/advection/advection_ext.pyf	(revision 4975)
+++ 	(revision )
@@ -1,36 +1,0 @@
-
-! File advection_ext.pyf
-python module advection_ext
-interface
-	function compute_fluxes(stage_update,stage_edge,stage_bdry,&
-			neighbours,neighbour_edges,&
-			normals,areas,radii,edgelengths,&
-			tri_full_flag,huge_timestep,max_timestep,v,n,nbdry)
-		intent(c) compute_fluxes
-		double precision compute_fluxes
-
-		double precision stage_update(n)	
-		intent(in,out)    stage_update
-		intent(c)         stage_update	
-		double precision stage_edge(n,3)
-		double precision stage_bdry(nbdry)
-                integer           neighbours(n,3)
-		integer           neighbour_edges(n,3)
-		double precision normals(n,6)
-                double precision areas(n)
-		double precision radii(n)
-		double precision edgelengths(n,3)
-		integer          tri_full_flag(n)
-                double precision huge_timestep, max_timestep
-                intent(c)        huge_timestep, max_timestep
-		double precision v(2)
-                integer          n, nbdry
-		intent(c)        n, nbdry
-
-
-	end function compute_fluxes
-end interface
-end python module
-
-
-
Index: anuga_core/source/anuga/advection/compile_all.py
===================================================================
--- anuga_core/source/anuga/advection/compile_all.py	(revision 4975)
+++ 	(revision )
@@ -1,14 +1,0 @@
-"""f2py script -
-
-"""
-
-import os
-
-separation_line = '---------------------------------------'
-print separation_line
-
-command = 'python setup.py build_src build_ext --inplace'
-print 'Trying to run %s in directory %s' % (command, os.getcwd())
-os.system(command)
-
-
Index: anuga_core/source/anuga/advection/setup.py
===================================================================
--- anuga_core/source/anuga/advection/setup.py	(revision 4975)
+++ 	(revision )
@@ -1,18 +1,0 @@
-"""f2py script -
-
-"""
-
-def configuration(parent_package='',top_path=None):
-    from numpy.distutils.misc_util import Configuration
-    config = Configuration('',parent_package,top_path)
-
-    config.add_extension('advection_ext',
-                         sources = ['advection_ext.pyf','advection_ext.c'])
-    return config
-
-#-------------------------------------------------------------
-if __name__ == "__main__":
-    from numpy.distutils.core import setup
-    setup(configuration=configuration)
-
-
Index: anuga_core/source/anuga/advection/test_advection.py
===================================================================
--- anuga_core/source/anuga/advection/test_advection.py	(revision 4975)
+++ anuga_core/source/anuga/advection/test_advection.py	(revision 4978)
@@ -60,8 +60,4 @@
         U = -domain.quantities['stage'].explicit_update
         R = -0.5/domain.areas[0]
-
-        print 'U ', U
-        print 'R ', R
-        
 
         assert U==R, '%s %s' %(U, R)
Index: anuga_core/source/anuga/advection/test_all.py
===================================================================
--- anuga_core/source/anuga/advection/test_all.py	(revision 4978)
+++ anuga_core/source/anuga/advection/test_all.py	(revision 4978)
@@ -0,0 +1,102 @@
+"""Regression testing framework
+This module will search for scripts in the same directory named
+test_*.py.  Each such script should be a test suite that tests a
+module through PyUnit. This script will aggregate all
+found test suites into one big test suite and run them all at once.
+"""
+
+# Author: Mark Pilgrim
+# Modified by Ole Nielsen
+
+import unittest
+import os
+
+
+#List files that should be excluded from the testing process.
+#E.g. if they are known to fail and under development
+exclude = [] 
+
+
+def get_test_files(path):
+
+    import sys
+
+    files = os.listdir(path)
+
+    #Check sub directories
+    test_files = []
+    for file in files:
+        if os.path.isdir(file):
+            sys.path.append(file)
+            #print 'Recursing into', file
+            test_files += get_test_files(path + os.sep + file)
+        elif file[:5] == 'test_' and file[-2:] == 'py':
+            #print 'Appending', file
+            test_files.append(file)
+        else:
+            pass
+    return test_files
+
+
+
+def regressionTest():
+    import sys, os, re, unittest
+    path = os.path.split(sys.argv[0])[0] or os.getcwd()
+
+
+    files = get_test_files(path)
+
+
+
+    #test = re.compile('^test_[\w]*.py$', re.IGNORECASE)
+    #files = filter(test.search, files)
+
+
+    try:
+        files.remove(__file__)  #Remove self from list (Ver 2.3. or later)
+    except:
+        files.remove('test_all.py')
+
+    print 'Testing:'
+    for file in files:
+        print '  ' + file
+
+    if globals().has_key('exclude'):
+        for file in exclude:
+            files.remove(file)
+            print 'WARNING: File '+ file + ' excluded from testing'
+
+
+    filenameToModuleName = lambda f: os.path.splitext(f)[0]
+    #print "files",files 
+    moduleNames = map(filenameToModuleName, files)
+    #print "moduleNames",moduleNames
+    modules = map(__import__, moduleNames)
+    load = unittest.defaultTestLoader.loadTestsFromModule
+    return unittest.TestSuite(map(load, modules))
+
+if __name__ == '__main__':
+
+    from os import sep
+
+    #Attempt to compile all extensions
+    #execfile('..' + sep + 'utilities' + sep + 'compile.py')
+
+    #FIXME: Temporary measure
+    #os.chdir('..' + sep + 'utilities')
+    #execfile('compile.py')
+    #os.chdir('..' + sep + 'pyvolution')    
+    
+    #FIXME: Temporary measure
+    #os.chdir('..' + sep + 'triangle')
+    #execfile('compile.py')
+    #os.chdir('..' + sep + 'pyvolution')   
+    
+    #os.system('python compile.py') 
+
+    #print regressionTest()
+    #unittest.main(defaultTest='regressionTest')
+    
+    suite = regressionTest()
+    runner = unittest.TextTestRunner() #(verbosity=2)
+    runner.run(suite)
