Index: /trunk/anuga_core/compile_all.py
===================================================================
--- /trunk/anuga_core/compile_all.py	(revision 8247)
+++ /trunk/anuga_core/compile_all.py	(revision 8248)
@@ -28,4 +28,8 @@
 
 os.chdir('..')
+os.chdir('geometry')
+execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
+
+os.chdir('..')
 os.chdir('structures')
 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
@@ -43,7 +47,4 @@
 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 
-os.chdir('..')
-os.chdir('shallow_water_balanced')
-execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 
 os.chdir('..')
Index: /trunk/anuga_core/source/anuga/__init__.py
===================================================================
--- /trunk/anuga_core/source/anuga/__init__.py	(revision 8247)
+++ /trunk/anuga_core/source/anuga/__init__.py	(revision 8248)
@@ -305,5 +305,5 @@
     """
 
-    from anuga.shallow_water.shallow_water_domain import Domain
+    #from anuga.shallow_water.shallow_water_domain import Domain
     from anuga.pmesh.mesh_interface import create_mesh_from_regions
     
@@ -321,4 +321,5 @@
                              use_cache=False,
                              verbose=verbose)
+
     domain = Domain(mesh_filename, use_cache=False, verbose=verbose)
 
Index: /trunk/anuga_core/source/anuga/compile_all.py
===================================================================
--- /trunk/anuga_core/source/anuga/compile_all.py	(revision 8247)
+++ /trunk/anuga_core/source/anuga/compile_all.py	(revision 8248)
@@ -10,10 +10,15 @@
 execfile('compile.py')
 
-os.chdir('..')
-os.chdir('advection')
-execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
+#os.chdir('..')
+#os.chdir('advection')
+#execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 
 os.chdir('..')
 os.chdir('operators')
+execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
+
+
+os.chdir('..')
+os.chdir('geometry')
 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 
@@ -34,7 +39,7 @@
 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 
-os.chdir('..')
-os.chdir('shallow_water_balanced')
-execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
+#os.chdir('..')
+#os.chdir('shallow_water_balanced')
+#execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 
 os.chdir('..')
Index: /trunk/anuga_core/source/anuga/operators/base_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/base_operator.py	(revision 8247)
+++ /trunk/anuga_core/source/anuga/operators/base_operator.py	(revision 8248)
@@ -28,6 +28,7 @@
 
     def __parallel_safe(self):
-
-        return True
+        """By default an operator is not parallel safe
+        """
+        return False
 
     def statistics(self):
@@ -35,4 +36,5 @@
         message = 'You need to implement operator statistics for your operator'
         return message
+    
 
     def print_statistics(self):
Index: /trunk/anuga_core/source/anuga/operators/mannings_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/mannings_operator.py	(revision 8248)
+++ /trunk/anuga_core/source/anuga/operators/mannings_operator.py	(revision 8248)
@@ -0,0 +1,60 @@
+# To change this template, choose Tools | Templates
+# and open the template in the editor.
+
+__author__="steve"
+__date__ ="$11/11/2011 1:52:17 PM$"
+
+from anuga.operators.base_operator import Operator
+
+
+
+class Mannings_operator(Operator):
+    """
+    Class for setting up a mannings opeator to apply Mannings fricition
+
+    Applying
+
+    d/dt uh =  uh  sqrt( uh^2 + vh^2 ) / h^{7/3}
+
+    d/dt vh =  vh  sqrt( uh^2 + vh^2 ) / h^{7/3}
+
+
+    """
+
+    def __init__(self, domain, verbose=False):
+        if verbose: log.critical('Mannings Operator: Beginning Initialisation')
+
+
+        Operator.__init__(self,domain)
+
+
+
+        if verbose: log.critical('Mannings Operator: Initialisation Done')
+
+
+    def __call__(self):
+
+        #timestep = self.domain.get_timestep()
+        raise Exception('Need to implement __call__ for your operator')
+
+
+    def __parallel_safe(self):
+        """
+        This operator only works with centroid values of quantities
+        and so is parallel safe.
+        """
+        return True
+
+
+    def statistics(self):
+
+        message = 'Manning Operator '
+        return message
+
+
+    def timestepping_statistics(self):
+
+        message = '    Manning Operator: '
+        return message
+
+
Index: /trunk/anuga_core/source/anuga/operators/mannings_operator_ext.c
===================================================================
--- /trunk/anuga_core/source/anuga/operators/mannings_operator_ext.c	(revision 8248)
+++ /trunk/anuga_core/source/anuga/operators/mannings_operator_ext.c	(revision 8248)
@@ -0,0 +1,317 @@
+// Python - C extension module for shallow_water.py
+//
+// To compile (Python2.3):
+//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
+//  gcc -shared domain_ext.o  -o domain_ext.so
+//
+// or use python compile.py
+//
+// See the module shallow_water_domain.py for more documentation on 
+// how to use this module
+//
+//
+// Ole Nielsen, GA 2004
+
+
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+#include <stdio.h>
+#include "numpy_shim.h"
+
+// Shared code snippets
+#include "util_ext.h"
+
+
+
+
+
+const double pi = 3.14159265358979;
+
+
+
+
+void _manning_friction_flat(double g, double eps, int N,
+               double* w, double* zv,
+               double* uh, double* vh,
+               double* eta, double* xmom, double* ymom) {
+
+  int k, k3;
+  double S, h, z, z0, z1, z2;
+
+  for (k=0; k<N; k++) {
+    if (eta[k] > eps) {
+      k3 = 3*k;
+      // Get bathymetry
+      z0 = zv[k3 + 0];
+      z1 = zv[k3 + 1];
+      z2 = zv[k3 + 2];
+      z = (z0+z1+z2)/3.0;
+      h = w[k]-z;
+      if (h >= eps) {
+        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
+        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
+        //S /= exp((7.0/3.0)*log(h));      //seems to save about 15% over manning_friction
+        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
+
+
+        //Update momentum
+        xmom[k] += S*uh[k];
+        ymom[k] += S*vh[k];
+      }
+    }
+  }
+}
+
+
+void _manning_friction_sloped(double g, double eps, int N,
+               double* x, double* w, double* zv,
+               double* uh, double* vh,
+               double* eta, double* xmom_update, double* ymom_update) {
+
+  int k, k3, k6;
+  double S, h, z, z0, z1, z2, zs, zx, zy;
+  double x0,y0,x1,y1,x2,y2;
+
+  for (k=0; k<N; k++) {
+    if (eta[k] > eps) {
+      k3 = 3*k;
+      // Get bathymetry
+      z0 = zv[k3 + 0];
+      z1 = zv[k3 + 1];
+      z2 = zv[k3 + 2];
+
+      // Compute bed slope
+      k6 = 6*k;  // base index
+  
+      x0 = x[k6 + 0];
+      y0 = x[k6 + 1];
+      x1 = x[k6 + 2];
+      y1 = x[k6 + 3];
+      x2 = x[k6 + 4];
+      y2 = x[k6 + 5];
+
+      _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
+
+      zs = sqrt(1.0 + zx*zx + zy*zy);
+      z = (z0+z1+z2)/3.0;
+      h = w[k]-z;
+      if (h >= eps) {
+        S = -g * eta[k]*eta[k] * zs * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
+        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
+        //S /= exp((7.0/3.0)*log(h));      //seems to save about 15% over manning_friction
+        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
+
+
+        //Update momentum
+        xmom_update[k] += S*uh[k];
+        ymom_update[k] += S*vh[k];
+      }
+    }
+  }
+}
+
+
+
+void _chezy_friction(double g, double eps, int N,
+               double* x, double* w, double* zv,
+               double* uh, double* vh,
+               double* chezy, double* xmom_update, double* ymom_update) {
+
+  int k, k3, k6;
+  double S, h, z, z0, z1, z2, zs, zx, zy;
+  double x0,y0,x1,y1,x2,y2;
+
+  for (k=0; k<N; k++) {
+    if (chezy[k] > eps) {
+      k3 = 3*k;
+      // Get bathymetry
+      z0 = zv[k3 + 0];
+      z1 = zv[k3 + 1];
+      z2 = zv[k3 + 2];
+
+      // Compute bed slope
+      k6 = 6*k;  // base index
+
+      x0 = x[k6 + 0];
+      y0 = x[k6 + 1];
+      x1 = x[k6 + 2];
+      y1 = x[k6 + 3];
+      x2 = x[k6 + 4];
+      y2 = x[k6 + 5];
+
+      _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
+
+      zs = sqrt(1.0 + zx*zx + zy*zy);
+      z = (z0+z1+z2)/3.0;
+      h = w[k]-z;
+      if (h >= eps) {
+        S = -g * chezy[k] * zs * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
+        S /= pow(h, 2.0);
+
+        //Update momentum
+        xmom_update[k] += S*uh[k];
+        ymom_update[k] += S*vh[k];
+      }
+    }
+  }
+}
+
+
+PyObject *manning_friction_sloped(PyObject *self, PyObject *args) {
+  //
+  // manning_friction_sloped(g, eps, x, h, uh, vh, z, eta, xmom_update, ymom_update)
+  //
+
+
+  PyArrayObject *x, *w, *z, *uh, *vh, *eta, *xmom, *ymom;
+  int N;
+  double g, eps;
+
+  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
+            &g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
+  }
+
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(x);
+  CHECK_C_CONTIG(w);
+  CHECK_C_CONTIG(z);
+  CHECK_C_CONTIG(uh);
+  CHECK_C_CONTIG(vh);
+  CHECK_C_CONTIG(eta);
+  CHECK_C_CONTIG(xmom);
+  CHECK_C_CONTIG(ymom);
+
+  N = w -> dimensions[0];
+
+  _manning_friction_sloped(g, eps, N,
+            (double*) x -> data,
+            (double*) w -> data,
+            (double*) z -> data,
+            (double*) uh -> data,
+            (double*) vh -> data,
+            (double*) eta -> data,
+            (double*) xmom -> data,
+            (double*) ymom -> data);
+
+  return Py_BuildValue("");
+}
+
+
+PyObject *chezy_friction(PyObject *self, PyObject *args) {
+  //
+  // chezy_friction(g, eps, x, h, uh, vh, z, chezy, xmom_update, ymom_update)
+  //
+
+
+  PyArrayObject *x, *w, *z, *uh, *vh, *chezy, *xmom, *ymom;
+  int N;
+  double g, eps;
+
+  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
+            &g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
+  }
+
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(x);
+  CHECK_C_CONTIG(w);
+  CHECK_C_CONTIG(z);
+  CHECK_C_CONTIG(uh);
+  CHECK_C_CONTIG(vh);
+  CHECK_C_CONTIG(chezy);
+  CHECK_C_CONTIG(xmom);
+  CHECK_C_CONTIG(ymom);
+
+  N = w -> dimensions[0];
+
+  _chezy_friction(g, eps, N,
+            (double*) x -> data,
+            (double*) w -> data,
+            (double*) z -> data,
+            (double*) uh -> data,
+            (double*) vh -> data,
+            (double*) chezy -> data,
+            (double*) xmom -> data,
+            (double*) ymom -> data);
+
+  return Py_BuildValue("");
+}
+
+
+PyObject *manning_friction_flat(PyObject *self, PyObject *args) {
+  //
+  // manning_friction_flat(g, eps, h, uh, vh, z, eta, xmom_update, ymom_update)
+  //
+
+
+  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
+  int N;
+  double g, eps;
+
+  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
+            &g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
+  }
+
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(w);
+  CHECK_C_CONTIG(z);
+  CHECK_C_CONTIG(uh);
+  CHECK_C_CONTIG(vh);
+  CHECK_C_CONTIG(eta);
+  CHECK_C_CONTIG(xmom);
+  CHECK_C_CONTIG(ymom);
+
+  N = w -> dimensions[0];
+
+  _manning_friction_flat(g, eps, N,
+            (double*) w -> data,
+            (double*) z -> data,
+            (double*) uh -> data,
+            (double*) vh -> data,
+            (double*) eta -> data,
+            (double*) xmom -> data,
+            (double*) ymom -> data);
+
+  return Py_BuildValue("");
+}
+
+
+
+//=========================================================================
+// Python Glue
+//=========================================================================
+
+
+
+
+
+
+
+
+//-------------------------------
+// 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.
+   */
+
+  {"manning_friction_flat", manning_friction_flat, METH_VARARGS, "Print out"},
+  {"manning_friction_sloped", manning_friction_sloped, METH_VARARGS, "Print out"},
+  {"chezy_friction", chezy_friction, METH_VARARGS, "Print out"},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initmannings_operator_ext(void){
+  Py_InitModule("mannings_operator_ext", MethodTable);
+
+  import_array(); // Necessary for handling of NumPY structures
+}
Index: /trunk/anuga_validation/okushiri_2005/run_okushiri.py
===================================================================
--- /trunk/anuga_validation/okushiri_2005/run_okushiri.py	(revision 8247)
+++ /trunk/anuga_validation/okushiri_2005/run_okushiri.py	(revision 8248)
@@ -22,10 +22,10 @@
 
 # Module imports
-from anuga.interface import Domain
-from anuga.interface import Reflective_boundary
-from anuga.interface import Time_boundary
-from anuga.interface import Transmissive_momentum_set_stage_boundary
-from anuga.interface import Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
-from anuga.interface import file_function
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import Time_boundary
+from anuga import Transmissive_momentum_set_stage_boundary
+from anuga import Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
+from anuga import file_function
 
 import project
@@ -84,5 +84,7 @@
 triangle_id = domain.get_triangle_containing_point([4.521, 1.196])
 # This should get triangle id 32833 with centroid (4.5244, 1.1972)
+print 'triangle_id = ',triangle_id
 bdry_id = domain.get_triangle_containing_point([0.000, 1.696])
+# This should get triangle id 56474
 print 'bdry_id = ',bdry_id
 
@@ -96,5 +98,5 @@
 for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
     domain.write_time()
-    print  function(domain.get_time())[0],' bdry_d = ',bdry_id,' ',\
+    print '   ',function(domain.get_time())[0],' bdry_d = ',bdry_id,' ',\
           domain.get_conserved_quantities(bdry_id, edge=0)[0],' ',\
           domain.get_conserved_quantities(bdry_id, edge=1)[0],' ',\
Index: /trunk/anuga_work/shallow_water_balanced_steve/run_step.py
===================================================================
--- /trunk/anuga_work/shallow_water_balanced_steve/run_step.py	(revision 8247)
+++ /trunk/anuga_work/shallow_water_balanced_steve/run_step.py	(revision 8248)
@@ -12,4 +12,5 @@
 import sys
 import anuga
+#from anuga import Domain
 from swb_domain import *
 
@@ -46,4 +47,19 @@
 # structured mesh
 points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+
+create_mesh_from_regions(bounding_polygon,
+                             boundary_tags,
+                             maximum_triangle_area=maximum_triangle_area,
+                             interior_regions=interior_regions,
+                             filename=mesh_filename,
+                             interior_holes=interior_holes,
+                             hole_tags=hole_tags,
+                             poly_geo_reference=poly_geo_reference,
+                             mesh_geo_reference=mesh_geo_reference,
+                             minimum_triangle_angle=minimum_triangle_angle,
+                             fail_if_polygons_outside=fail_if_polygons_outside,
+                             use_cache=False,
+                             verbose=verbose)
 
 domain = Domain(points, vertices, boundary) 
Index: /trunk/anuga_work/shallow_water_balanced_steve/run_wave.py
===================================================================
--- /trunk/anuga_work/shallow_water_balanced_steve/run_wave.py	(revision 8247)
+++ /trunk/anuga_work/shallow_water_balanced_steve/run_wave.py	(revision 8248)
@@ -12,5 +12,6 @@
 import sys
 import anuga
-from swb_domain import *
+from anuga import Domain
+#from swb_domain import *
 
 from math import cos
@@ -56,4 +57,5 @@
 domain.set_timestepping_method('rk2')
 domain.set_default_order(2)
+domain.set_beta(2.0)
 
 print domain.get_timestepping_method()
@@ -79,5 +81,5 @@
 Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
 amplitude = 1
-wave_length = 1000.0
+wave_length = 300.0
 Bw = anuga.Time_boundary(domain=domain,     # Time dependent boundary
 ## Sine wave
Index: /trunk/anuga_work/shallow_water_balanced_steve/swb_domain.py
===================================================================
--- /trunk/anuga_work/shallow_water_balanced_steve/swb_domain.py	(revision 8247)
+++ /trunk/anuga_work/shallow_water_balanced_steve/swb_domain.py	(revision 8248)
@@ -77,4 +77,5 @@
         self.forcing_terms[1] = swb_gravity
 
+        print "Swb_domain"
 
     def check_integrity(self):
