Index: /trunk/anuga_core/compile_all.py
===================================================================
--- /trunk/anuga_core/compile_all.py	(revision 8128)
+++ /trunk/anuga_core/compile_all.py	(revision 8129)
@@ -21,4 +21,12 @@
 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('structures')
 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 
Index: /trunk/anuga_core/source/anuga/__init__.py
===================================================================
--- /trunk/anuga_core/source/anuga/__init__.py	(revision 8128)
+++ /trunk/anuga_core/source/anuga/__init__.py	(revision 8129)
@@ -320,5 +320,10 @@
         #psyco.background() # attempt to profile code - only compile most used
         
-        
+
+#---------------------------
+# Operators
+#---------------------------
+from anuga.operators.operator import Operator
+
 #---------------------------
 # Structures
Index: /trunk/anuga_core/source/anuga/compile_all.py
===================================================================
--- /trunk/anuga_core/source/anuga/compile_all.py	(revision 8128)
+++ /trunk/anuga_core/source/anuga/compile_all.py	(revision 8129)
@@ -12,4 +12,12 @@
 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('structures')
 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 
Index: /trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 8128)
+++ /trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 8129)
@@ -71,7 +71,7 @@
     if (dq1>=dq2){
       if (dq1>=0.0)
-	*qmax=dq0+dq1;
+    *qmax=dq0+dq1;
       else
-	*qmax=dq0;
+    *qmax=dq0;
       
       *qmin=dq0+dq2;
@@ -80,7 +80,7 @@
     else{// dq1<dq2
       if (dq2>0)
-	*qmax=dq0+dq2;
+    *qmax=dq0+dq2;
       else
-	*qmax=dq0;
+    *qmax=dq0;
       
       *qmin=dq0+dq1;    
@@ -91,7 +91,7 @@
     if (dq1<=dq2){
       if (dq1<0.0)
-	*qmin=dq0+dq1;
+    *qmin=dq0+dq1;
       else
-	*qmin=dq0;
+    *qmin=dq0;
       
       *qmax=dq0+dq2;    
@@ -100,7 +100,7 @@
     else{// dq1>dq2
       if (dq2<0.0)
-	*qmin=dq0+dq2;
+    *qmin=dq0+dq2;
       else
-	*qmin=dq0;
+    *qmin=dq0;
       
       *qmax=dq0+dq1;
@@ -175,8 +175,8 @@
 //static inline double _compute_speed(double *uh, 
 double _compute_speed(double *uh, 
-		      double *h, 
-		      double epsilon, 
-		      double h0,
-		      double limiting_threshold) {
+              double *h,
+              double epsilon,
+              double h0,
+              double limiting_threshold) {
   
   double u;
@@ -263,7 +263,7 @@
                            double n1, double n2,
                            double epsilon, 
-			   double h0,
-			   double limiting_threshold, 
-			   double g,
+               double h0,
+               double limiting_threshold,
+               double g,
                            double *edgeflux, double *max_speed) 
 {
@@ -304,8 +304,8 @@
   z = 0.5*(z_left + z_right); // Average elevation values. 
                               // Even though this will nominally allow 
-			      // for discontinuities in the elevation data, 
-			      // there is currently no numerical support for
-			      // this so results may be strange near 
-			      // jumps in the bed.
+                  // for discontinuities in the elevation data,
+                  // there is currently no numerical support for
+                  // this so results may be strange near
+                  // jumps in the bed.
 
   // Compute speeds in x-direction
@@ -314,5 +314,5 @@
   uh_left = q_left_rotated[1];
   u_left = _compute_speed(&uh_left, &h_left, 
-			  epsilon, h0, limiting_threshold);
+              epsilon, h0, limiting_threshold);
 
   w_right = q_right_rotated[0];
@@ -320,5 +320,5 @@
   uh_right = q_right_rotated[1];
   u_right = _compute_speed(&uh_right, &h_right, 
-			   epsilon, h0, limiting_threshold);
+               epsilon, h0, limiting_threshold);
 
   // Momentum in y-direction
@@ -330,7 +330,7 @@
   // All validation tests pass, so do we really need it anymore? 
   _compute_speed(&vh_left, &h_left, 
-		 epsilon, h0, limiting_threshold);
+         epsilon, h0, limiting_threshold);
   _compute_speed(&vh_right, &h_right, 
-		 epsilon, h0, limiting_threshold);
+         epsilon, h0, limiting_threshold);
 
   // Maximal and minimal wave speeds
@@ -459,5 +459,5 @@
   uh_left = q_left_rotated[1];
   u_left =_compute_speed(&uh_left, &h_left, 
-			 epsilon, h0, limiting_threshold);
+             epsilon, h0, limiting_threshold);
 
   w_right = q_right_rotated[0];
@@ -465,5 +465,5 @@
   uh_right = q_right_rotated[1];
   u_right =_compute_speed(&uh_right, &h_right,
-			  epsilon, h0, limiting_threshold);
+              epsilon, h0, limiting_threshold);
 
 
@@ -553,7 +553,7 @@
 
 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) {
+               double* x, double* w, double* zv,
+               double* uh, double* vh,
+               double* eta, double* xmom_update, double* ymom_update) {
 
   int k, k3, k6;
@@ -602,7 +602,7 @@
 
 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) {
+               double* x, double* w, double* zv,
+               double* uh, double* vh,
+               double* chezy, double* xmom_update, double* ymom_update) {
 
   int k, k3, k6;
@@ -785,5 +785,5 @@
       // FIXME: Try with this one precomputed
       for (i=0; i<3; i++) {
-	dz = max(dz, fabs(zv[k3+i]-zc[k]));
+    dz = max(dz, fabs(zv[k3+i]-zc[k]));
       }
     }
@@ -813,7 +813,7 @@
       
       if (dz > 0.0) {
-	alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );      
+    alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );
       } else {
-	alpha = 1.0;  // Flat bed
+    alpha = 1.0;  // Flat bed
       }
       //printf("Using old style limiter\n");
@@ -828,27 +828,27 @@
     
       if (hmin < H0) {
-	alpha = 1.0;
-	for (i=0; i<3; i++) {
-	  
-	  h_diff = hc_k - hv[i];      
-	  if (h_diff <= 0) {
-	    // Deep water triangle is further away from bed than 
-	    // shallow water (hbar < h). Any alpha will do
-      
-	  } else {  
-	    // Denominator is positive which means that we need some of the 
-	    // h-limited stage.
-	    
-	    alpha = min(alpha, (hc_k - H0)/h_diff);     
-	  }
-	}
-
-	// Ensure alpha in [0,1]
-	if (alpha>1.0) alpha=1.0;
-	if (alpha<0.0) alpha=0.0;
-	
+    alpha = 1.0;
+    for (i=0; i<3; i++) {
+
+      h_diff = hc_k - hv[i];
+      if (h_diff <= 0) {
+        // Deep water triangle is further away from bed than
+        // shallow water (hbar < h). Any alpha will do
+      
       } else {
-	// Use w-limited stage exclusively in deeper water.
-	alpha = 1.0;       
+        // Denominator is positive which means that we need some of the
+        // h-limited stage.
+
+        alpha = min(alpha, (hc_k - H0)/h_diff);
+      }
+    }
+
+    // Ensure alpha in [0,1]
+    if (alpha>1.0) alpha=1.0;
+    if (alpha<0.0) alpha=0.0;
+
+      } else {
+    // Use w-limited stage exclusively in deeper water.
+    alpha = 1.0;
       }
     }
@@ -875,44 +875,44 @@
     if (alpha < 1) {      
       for (i=0; i<3; i++) {
-	
-	wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i]; 
-
-	// Update momentum at vertices
-	if (use_centroid_velocities == 1) {
-	  // This is a simple, efficient and robust option
-	  // It uses first order approximation of velocities, but retains
-	  // the order used by stage.
+
+    wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
+
+    // Update momentum at vertices
+    if (use_centroid_velocities == 1) {
+      // This is a simple, efficient and robust option
+      // It uses first order approximation of velocities, but retains
+      // the order used by stage.
     
-	  // Speeds at centroids
-	  if (hc_k > epsilon) {
-	    uc = xmomc[k]/hc_k;
-	    vc = ymomc[k]/hc_k;
-	  } else {
-	    uc = 0.0;
-	    vc = 0.0;
-	  }
-	  
-	  // Vertex momenta guaranteed to be consistent with depth guaranteeing
-	  // controlled speed
-	  hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
-	  xmomv[k3+i] = uc*hv[i];
-	  ymomv[k3+i] = vc*hv[i];
-      
-	} else {
-	  // Update momentum as a linear combination of
-	  // xmomc and ymomc (shallow) and momentum
-	  // from extrapolator xmomv and ymomv (deep).
-	  // This assumes that values from xmomv and ymomv have
-	  // been established e.g. by the gradient limiter.
-	  
-	  // FIXME (Ole): I think this should be used with vertex momenta
-	  // computed above using centroid_velocities instead of xmomc 
-	  // and ymomc as they'll be more representative first order
-	  // values.
-	  
-	  xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
-	  ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
-	  
-	}
+      // Speeds at centroids
+      if (hc_k > epsilon) {
+        uc = xmomc[k]/hc_k;
+        vc = ymomc[k]/hc_k;
+      } else {
+        uc = 0.0;
+        vc = 0.0;
+      }
+
+      // Vertex momenta guaranteed to be consistent with depth guaranteeing
+      // controlled speed
+      hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
+      xmomv[k3+i] = uc*hv[i];
+      ymomv[k3+i] = vc*hv[i];
+      
+    } else {
+      // Update momentum as a linear combination of
+      // xmomc and ymomc (shallow) and momentum
+      // from extrapolator xmomv and ymomv (deep).
+      // This assumes that values from xmomv and ymomv have
+      // been established e.g. by the gradient limiter.
+
+      // FIXME (Ole): I think this should be used with vertex momenta
+      // computed above using centroid_velocities instead of xmomc
+      // and ymomc as they'll be more representative first order
+      // values.
+
+      xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
+      ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
+
+    }
       }
     }
@@ -945,8 +945,8 @@
       if (hc < minimum_allowed_height) {
         
-	// Set momentum to zero and ensure h is non negative
-	xmomc[k] = 0.0;
-	ymomc[k] = 0.0;
-	if (hc <= 0.0) wc[k] = zc[k];
+    // Set momentum to zero and ensure h is non negative
+    xmomc[k] = 0.0;
+    ymomc[k] = 0.0;
+    if (hc <= 0.0) wc[k] = zc[k];
       }
     }
@@ -968,22 +968,22 @@
         } else {
           //Reduce excessive speeds derived from division by small hc
-	  //FIXME (Ole): This may be unnecessary with new slope limiters 
-	  //in effect.
+      //FIXME (Ole): This may be unnecessary with new slope limiters
+      //in effect.
           
           u = xmomc[k]/hc;
-	  if (fabs(u) > maximum_allowed_speed) {
-	    reduced_speed = maximum_allowed_speed * u/fabs(u);
-	    //printf("Speed (u) has been reduced from %.3f to %.3f\n",
-	    //   u, reduced_speed);
-	    xmomc[k] = reduced_speed * hc;
-	  }
-	  
+      if (fabs(u) > maximum_allowed_speed) {
+        reduced_speed = maximum_allowed_speed * u/fabs(u);
+        //printf("Speed (u) has been reduced from %.3f to %.3f\n",
+        //   u, reduced_speed);
+        xmomc[k] = reduced_speed * hc;
+      }
+
           v = ymomc[k]/hc;
-	  if (fabs(v) > maximum_allowed_speed) {
-	    reduced_speed = maximum_allowed_speed * v/fabs(v);
-	    //printf("Speed (v) has been reduced from %.3f to %.3f\n",
-	    //   v, reduced_speed);
-	    ymomc[k] = reduced_speed * hc;
-	  }
+      if (fabs(v) > maximum_allowed_speed) {
+        reduced_speed = maximum_allowed_speed * v/fabs(v);
+        //printf("Speed (v) has been reduced from %.3f to %.3f\n",
+        //   v, reduced_speed);
+        ymomc[k] = reduced_speed * hc;
+      }
         }
       }
@@ -1056,5 +1056,5 @@
               // But evidence suggests that h0 can be as little as
               // epsilon!
-	      
+
   limiting_threshold = 10*H0; // Avoid applying limiter below this
                               // threshold for performance reasons.
@@ -1062,13 +1062,13 @@
   
   _flux_function_central((double*) ql -> data, 
-			 (double*) qr -> data, 
-			 zl, 
-			 zr,                         
-			 ((double*) normal -> data)[0],
-			 ((double*) normal -> data)[1],          
-			 epsilon, h0, limiting_threshold,
-			 g,
-			 (double*) edgeflux -> data, 
-			 &max_speed);
+             (double*) qr -> data,
+             zl,
+             zr,
+             ((double*) normal -> data)[0],
+             ((double*) normal -> data)[1],
+             epsilon, h0, limiting_threshold,
+             g,
+             (double*) edgeflux -> data,
+             &max_speed);
   
   return Py_BuildValue("d", max_speed);  
@@ -1156,5 +1156,5 @@
 
   if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
-			&g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
+            &g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
       report_python_error(AT, "could not parse input arguments");
       return NULL;
@@ -1174,12 +1174,12 @@
 
   _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);
+            (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("");
@@ -1198,5 +1198,5 @@
 
   if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
-			&g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
+            &g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
       report_python_error(AT, "could not parse input arguments");
       return NULL;
@@ -1216,12 +1216,12 @@
 
   _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);
+            (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("");
@@ -1240,5 +1240,5 @@
 
   if (!PyArg_ParseTuple(args, "ddOOOOOOO",
-			&g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
+            &g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
       report_python_error(AT, "could not parse input arguments");
       return NULL;
@@ -1257,11 +1257,11 @@
 
   _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);
+            (double*) w -> data,
+            (double*) z -> data,
+            (double*) uh -> data,
+            (double*) vh -> data,
+            (double*) eta -> data,
+            (double*) xmom -> data,
+            (double*) ymom -> data);
 
   return Py_BuildValue("");
@@ -1484,5 +1484,5 @@
       dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
       
-	  inv_area2 = 1.0/area2;
+      inv_area2 = 1.0/area2;
       // Calculate the gradient of stage on the auxiliary triangle
       a = dy2*dq1 - dy1*dq2;
@@ -1517,6 +1517,6 @@
       //for (i=0;i<3;i++)
       stage_vertex_values[k3+0] = stage_centroid_values[k] + dqv[0];
-	  stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1];
-	  stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2];
+      stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1];
+      stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2];
       
       
Index: /trunk/anuga_core/source/anuga/structures/structure_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/structure_operator.py	(revision 8128)
+++ /trunk/anuga_core/source/anuga/structures/structure_operator.py	(revision 8129)
@@ -6,8 +6,9 @@
 from anuga.utilities.system_tools import log_to_file
 from anuga.utilities.numerical_tools import ensure_numeric
-
-
-
-class Structure_operator:
+from anuga.operators.operator import Operator
+
+
+
+class Structure_operator(anuga.Operator):
     """Structure Operator - transfer water from one rectangular box to another.
     Sets up the geometry of problem
@@ -37,15 +38,10 @@
                  logging,
                  verbose):
-        
-        self.domain = domain
-        self.domain.set_fractional_step_operator(self)
+
+        anuga.Operator.__init__(self,domain)
+        
         self.end_points = ensure_numeric(end_points)
         self.exchange_lines = ensure_numeric(exchange_lines)
         self.enquiry_points = ensure_numeric(enquiry_points)
-
-
-        if domain.numproc > 1:
-            msg = 'Not implemented to run in parallel'
-            assert self.__parallel_safe(), msg
 
         
@@ -259,8 +255,4 @@
             self.enquiry_points.append(centre_point1 + gap)
             
-
-    def __parallel_safe(self):
-
-        return False
 
     def discharge_routine(self):
Index: /trunk/anuga_core/source/anuga/structures/test_boyd_box_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_boyd_box_operator.py	(revision 8128)
+++ /trunk/anuga_core/source/anuga/structures/test_boyd_box_operator.py	(revision 8129)
@@ -3,12 +3,9 @@
 
 import unittest
-import os.path
-import sys
-
-from anuga.utilities.system_tools import get_pathname_from_package
+
+
 from anuga.structures.boyd_box_operator import Boyd_box_operator
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
 from anuga.shallow_water.shallow_water_domain import Domain
-from anuga.shallow_water.forcing import Rainfall, Inflow
 import numpy
 
Index: /trunk/anuga_core/source/anuga/structures/test_inlet_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/test_inlet_operator.py	(revision 8128)
+++ /trunk/anuga_core/source/anuga/structures/test_inlet_operator.py	(revision 8129)
@@ -155,8 +155,20 @@
         finaltime = 3.0
         line1 = [[95.0, 10.0], [105.0, 10.0]]
-        Q1 = file_function('inlet_operator_test1.tms', quantities=['hydrograph'])
+
+        
+        import os
+        baseDir = os.getcwd()
+
+        try:
+            os.chdir('structures')
+        except:
+            pass
+
+        Q1 = file_function(filename='inlet_operator_test1.tms', quantities=['hydrograph'])
         
         line2 = [[10.0, 90.0], [20.0, 90.0]]
-        Q2 = file_function('inlet_operator_test2.tms', quantities=['hydrograph'])
+        Q2 = file_function(filename='inlet_operator_test2.tms', quantities=['hydrograph'])
+
+        os.chdir(baseDir)
         
         Inlet_operator(domain, line1, Q1)
Index: /trunk/anuga_core/source/anuga/test_all.py
===================================================================
--- /trunk/anuga_core/source/anuga/test_all.py	(revision 8128)
+++ /trunk/anuga_core/source/anuga/test_all.py	(revision 8129)
@@ -23,5 +23,5 @@
 
 # Directories that should not be searched for test files.
-exclude_dirs = ['pypar_dist',  'shallow_water_balanced',  'structures', # Special requirements
+exclude_dirs = ['shallow_water_balanced', 'operators' , # Special requirements
                 '.svn',          # subversion
                 'props', 'wcprops', 'prop-base', 'text-base', 'tmp']
Index: /trunk/anuga_core/test_all.py
===================================================================
--- /trunk/anuga_core/test_all.py	(revision 8128)
+++ /trunk/anuga_core/test_all.py	(revision 8129)
