Changeset 9306

Timestamp:

Aug 15, 2014, 4:48:51 PM (11 years ago)

Author:

davies

Message:

Cleaning up C interface for DE functions

Location:

trunk/anuga_core/source/anuga/shallow_water

Files:

• shallow_water_domain.py (modified) (3 diffs)
• swDE1_domain_ext.c (modified) (72 diffs)
• sw_domain.h (modified) (20 diffs)

trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -1697 +1697 @@
             # Flux calculation and gravity incorporated in same
             # procedure
-            # Also, fluxes are limited to prevent cells exporting more water
-            # than they contain
-
-            # FIXME SR: This needs cleaning up, should just be passing through
-            # the domain as in other compute flux calls
+
             from swDE1_domain_ext import compute_fluxes_ext_central \
                                       as compute_fluxes_ext
-            #if self.timestepping_method=='euler':
-            #    raise Exception, "DE1 doesn't seems to work with euler at present \
-            #                     (boundary conditions?), and is mostly designed for rk2"
-            Stage = self.quantities['stage']
-            Xmom = self.quantities['xmomentum']
-            Ymom = self.quantities['ymomentum']
-            Bed = self.quantities['elevation']
-            height = self.quantities['height']
 
             timestep = self.evolve_max_timestep 
 
-            flux_timestep = compute_fluxes_ext(timestep,
-                                               self.epsilon,
-                                               self.minimum_allowed_height,
-                                               self.g,
-                                               self.boundary_flux_sum,
-                                               self.neighbours,
-                                               self.neighbour_edges,
-                                               self.normals,
-                                               self.edgelengths,
-                                               self.radii,
-                                               self.areas,
-                                               self.tri_full_flag,
-                                               Stage.edge_values,
-                                               Xmom.edge_values,
-                                               Ymom.edge_values,
-                                               Bed.edge_values,
-                                               height.edge_values,
-                                               Stage.boundary_values,
-                                               Xmom.boundary_values,
-                                               Ymom.boundary_values,
-                                               self.boundary_flux_type,
-                                               Stage.explicit_update,
-                                               Xmom.explicit_update,
-                                               Ymom.explicit_update,
-                                               self.already_computed_flux,
-                                               self.max_speed,
-                                               int(self.optimise_dry_cells),
-                                               self.timestep_fluxcalls,
-                                               Stage.centroid_values,
-                                               Xmom.centroid_values,
-                                               Ymom.centroid_values,
-                                               Bed.centroid_values, 
-                                               height.centroid_values,
-                                               self.edge_flux_type,
-                                               self.riverwallData.riverwall_elevation,
-                                               self.riverwallData.hydraulic_properties_rowIndex,
-                                               int(self.riverwallData.ncol_hydraulic_properties),
-                                               self.riverwallData.hydraulic_properties,
-                                               self.edge_flux_work,
-                                               self.pressuregrad_work,
-                                               self.edge_timestep,
-                                               self.update_next_flux,
-                                               self.allow_timestep_increase)
+            flux_timestep = compute_fluxes_ext(self, timestep)
 
             self.flux_timestep = flux_timestep
@@ -1834 +1780 @@
             # Do extrapolation step
             from swDE1_domain_ext import extrapolate_second_order_edge_sw as extrapol2
-
-            Stage = self.quantities['stage']
-            Xmom = self.quantities['xmomentum']
-            Ymom = self.quantities['ymomentum']
-            Elevation = self.quantities['elevation']
-            height=self.quantities['height']
-
-            extrapol2(self,
-                      self.surrogate_neighbours,
-                      self.neighbour_edges,
-                      self.number_of_boundaries,
-                      self.centroid_coordinates,
-                      Stage.centroid_values,
-                      Xmom.centroid_values,
-                      Ymom.centroid_values,
-                      Elevation.centroid_values,
-                      height.centroid_values,
-                      self.edge_coordinates,
-                      Stage.edge_values,
-                      Xmom.edge_values,
-                      Ymom.edge_values,
-                      Elevation.edge_values,
-                      height.edge_values,
-                      Stage.vertex_values,
-                      Xmom.vertex_values,
-                      Ymom.vertex_values,
-                      Elevation.vertex_values,
-                      height.vertex_values,
-                      int(self.optimise_dry_cells),
-                      int(self.extrapolate_velocity_second_order),
-                      self.x_centroid_work,
-                      self.y_centroid_work,
-                      self.update_extrapolation)
+            extrapol2(self)
+
         else:
             # Code for original method
@@ -2633 +2548 @@
                                   as compute_flux_update_frequency_ext
 
-        compute_flux_update_frequency_ext(
-                self.timestep,
-                self.neighbours,
-                self.neighbour_edges,
-                self.already_computed_flux,
-                self.edge_timestep,
-                self.flux_update_frequency,
-                self.update_extrapolation,
-                self.max_flux_update_frequency,
-                self.update_next_flux,
-                self.allow_timestep_increase)
+        compute_flux_update_frequency_ext(self, self.timestep)
         
     def report_water_volume_statistics(self, verbose=True, returnStats=False):

trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c

@@ -25 +25 @@
 // Shared code snippets
 #include "util_ext.h"
+#include "sw_domain.h"
 
 
@@ -451 +452 @@
 ////////////////////////////////////////////////////////////////
 
-int _compute_flux_update_frequency(int number_of_elements, 
-                                  long *neighbours, long *neighbour_edges,
-                                  long *already_computed_flux,
-                                  int max_flux_update_frequency, double *edge_timestep,
-                                  double timestep, double notSoFast,
-                                  long *flux_update_frequency,
-                                  long *update_extrapolation, 
-                                  long *update_next_flux,
-                                  long *allow_timestep_increase){
+int _compute_flux_update_frequency(struct domain *D, double timestep){
     // Compute the 'flux_update_frequency' for each edge.
     //
@@ -474 +467 @@
     int k, i, k3, ki, m, n, nm, ii, j, ii2;
     long fuf;
+    double notSoFast=1.0;
     static int cyclic_number_of_steps=-1;
 
     // QUICK EXIT
-    if(max_flux_update_frequency==1){
+    if(D->max_flux_update_frequency==1){
         return 0;
     }
@@ -483 +477 @@
     // Count the steps
     cyclic_number_of_steps++;
-    if(cyclic_number_of_steps==max_flux_update_frequency){
+    if(cyclic_number_of_steps==D->max_flux_update_frequency){
         // The flux was just updated in every cell
         cyclic_number_of_steps=0;
@@ -491 +485 @@
     // PART 1: ONLY OCCURS FOLLOWING FLUX UPDATE
 
-    for ( k = 0; k < number_of_elements; k++){
+    for ( k = 0; k < D->number_of_elements; k++){
         for ( i = 0; i < 3; i++){
             ki = k*3 + i;
-            if((Mod_of_power_2(cyclic_number_of_steps, flux_update_frequency[ki])==0)){
+            if((Mod_of_power_2(cyclic_number_of_steps, D->flux_update_frequency[ki])==0)){
                 // The flux was just updated, along with the edge_timestep
                 // So we better recompute the flux_update_frequency
-                n=neighbours[ki];
+                n=D->neighbours[ki];
                 if(n>=0){
-                    m = neighbour_edges[ki];
+                    m = D->neighbour_edges[ki];
                     nm = n * 3 + m; // Linear index (triangle n, edge m)
                 }
@@ -506 +500 @@
                 // (Multiply already_computed_flux by -1 on the first update,
                 // and again on the 2nd)
-                if(already_computed_flux[ki] > 0 ){
+                if(D->already_computed_flux[ki] > 0 ){
                     // We have not fixed this flux value yet
-                    already_computed_flux[ki] *=-1;
+                    D->already_computed_flux[ki] *=-1;
                     if(n>=0){
-                        already_computed_flux[nm] *=-1;
+                        D->already_computed_flux[nm] *=-1;
                     }
                 }else{
                     // We have fixed this flux value already 
-                    already_computed_flux[ki] *=-1;
+                    D->already_computed_flux[ki] *=-1;
                     if(n>=0){
-                        already_computed_flux[nm] *=-1;
+                        D->already_computed_flux[nm] *=-1;
                     }
                     continue;
@@ -525 +519 @@
                 // NOTE: edge_timestep[ki]/timestep can be very large [so int overflows].
                 //       Do not pull the (int) inside the min term
-                fuf = (int)min((edge_timestep[ki]/timestep)*notSoFast,max_flux_update_frequency*1.);
+                fuf = (int)min((D->edge_timestep[ki]/timestep)*notSoFast,D->max_flux_update_frequency*1.);
                 // Account for neighbour
                 if(n>=0){
-                    fuf = min( (int)min(edge_timestep[nm]/timestep*notSoFast, max_flux_update_frequency*1.), fuf);
+                    fuf = min( (int)min(D->edge_timestep[nm]/timestep*notSoFast, D->max_flux_update_frequency*1.), fuf);
                 }
 
@@ -536 +530 @@
                 }
                 // Deal with large fuf 
-                if(fuf>max_flux_update_frequency){
-                    fuf = max_flux_update_frequency;
+                if(fuf> D->max_flux_update_frequency){
+                    fuf = D->max_flux_update_frequency;
                 }
                 //// Deal with intermediate cases
                 ii=2;
-                while(ii<max_flux_update_frequency){
+                while(ii< D->max_flux_update_frequency){
                     // Set it to 1,2,4, 8, ...
                     ii2=2*ii;
@@ -552 +546 @@
 
                 // Set the values
-                flux_update_frequency[ki]=fuf;
+                D->flux_update_frequency[ki]=fuf;
                 if(n>=0){
-                    flux_update_frequency[nm]= fuf;
+                    D->flux_update_frequency[nm]= fuf;
                 }
                 
@@ -567 +561 @@
     // Experiences suggests this is numerically important
     // (But, it can result in the same edge having different flux_update_freq)
-    for( k=0; k<number_of_elements; k++){
+    for( k=0; k< D->number_of_elements; k++){
         k3=3*k;
-        ii = 1*min(flux_update_frequency[k3],
-                 min(flux_update_frequency[k3+1],
-                     flux_update_frequency[k3+2]));
+        ii = 1*min(D->flux_update_frequency[k3],
+                 min(D->flux_update_frequency[k3+1],
+                     D->flux_update_frequency[k3+2]));
         
-        flux_update_frequency[k3]=min(ii, flux_update_frequency[k3]);
-        flux_update_frequency[k3+1]=min(ii, flux_update_frequency[k3+1]);
-        flux_update_frequency[k3+2]=min(ii,flux_update_frequency[k3+2]);
+        D->flux_update_frequency[k3]=min(ii, D->flux_update_frequency[k3]);
+        D->flux_update_frequency[k3+1]=min(ii, D->flux_update_frequency[k3+1]);
+        D->flux_update_frequency[k3+2]=min(ii,D->flux_update_frequency[k3+2]);
             
     }
@@ -585 +579 @@
     // will have its flux_update_frequency updated after a few timesteps (i.e. before max_flux_update_frequency timesteps)
     // OTOH, could this cause oscillations in flux_update_frequency?
-    for( k = 0; k < number_of_elements; k++){
-        update_extrapolation[k]=0;
+    for( k = 0; k < D->number_of_elements; k++){
+        D->update_extrapolation[k]=0;
         for( i = 0; i< 3; i++){
             ki=3*k+i;
             // Account for neighbour
-            n=neighbours[ki];
+            n=D->neighbours[ki];
             if(n>=0){
-                m = neighbour_edges[ki];
+                m = D->neighbour_edges[ki];
                 nm = n * 3 + m; // Linear index (triangle n, edge m)
-                flux_update_frequency[ki]=min(flux_update_frequency[ki], flux_update_frequency[nm]);
+                D->flux_update_frequency[ki]=min(D->flux_update_frequency[ki], D->flux_update_frequency[nm]);
             }
             // Do we need to update the extrapolation?
             // (We do if the next flux computation will actually compute a flux!)
-            if(Mod_of_power_2((cyclic_number_of_steps+1),flux_update_frequency[ki])==0){
-                update_next_flux[ki]=1;
-                update_extrapolation[k]=1;
+            if(Mod_of_power_2((cyclic_number_of_steps+1),D->flux_update_frequency[ki])==0){
+                D->update_next_flux[ki]=1;
+                D->update_extrapolation[k]=1;
             }else{
-                update_next_flux[ki]=0;
+                D->update_next_flux[ki]=0;
             }
         }
@@ -608 +602 @@
 
     // Check whether the timestep can be increased in the next compute_fluxes call
-    if(cyclic_number_of_steps+1==max_flux_update_frequency){
+    if(cyclic_number_of_steps+1==D->max_flux_update_frequency){
         // All fluxes will be updated on the next timestep
         // We also allow the timestep to increase then
-        allow_timestep_increase[0]=1;
+        D->allow_timestep_increase[0]=1;
     }else{
-        allow_timestep_increase[0]=0;
+        D->allow_timestep_increase[0]=0;
     }
 
@@ -712 +706 @@
 
 // Computational function for flux computation
-double _compute_fluxes_central(int number_of_elements,
-        double timestep,
-        double epsilon,
-        double H0,
-        double g,
-        double* boundary_flux_sum,
-        long* neighbours,
-        long* neighbour_edges,
-        double* normals,
-        double* edgelengths,
-        double* radii,
-        double* areas,
-        long* tri_full_flag,
-        double* stage_edge_values,
-        double* xmom_edge_values,
-        double* ymom_edge_values,
-        double* bed_edge_values,
-        double* height_edge_values,
-        double* stage_boundary_values,
-        double* xmom_boundary_values,
-        double* ymom_boundary_values,
-        long*   boundary_flux_type,
-        double* stage_explicit_update,
-        double* xmom_explicit_update,
-        double* ymom_explicit_update,
-        long* already_computed_flux,
-        double* max_speed_array,
-        int optimise_dry_cells, 
-        int timestep_fluxcalls,
-        double* stage_centroid_values,
-        double* xmom_centroid_values,
-        double* ymom_centroid_values,
-        double* bed_centroid_values,
-        double* height_centroid_values,
-        long* edge_flux_type,
-        double* riverwall_elevation,
-        long* riverwall_rowIndex,
-        int ncol_riverwall_hydraulic_properties,
-        double* riverwall_hydraulic_properties,
-        double* edge_flux_work,
-        double* pressuregrad_work,
-        double* edge_timestep,
-        long* update_next_flux,
-        long* allow_timestep_increase) {
+double _compute_fluxes_central(struct domain *D, double timestep){
+
     // Local variables
-    double max_speed, length, inv_area, zl, zr;
+    double max_speed_local, length, inv_area, zl, zr;
     double h_left, h_right, z_half ;  // For andusse scheme
     // FIXME: limiting_threshold is not used for DE1
-    double limiting_threshold = 10*H0;
+    double limiting_threshold = 10*D->H0;
     //
     int k, i, m, n,j, ii;
@@ -781 +733 @@
     // Set explicit_update to zero for all conserved_quantities.
     // This assumes compute_fluxes called before forcing terms
-    memset((char*) stage_explicit_update, 0, number_of_elements * sizeof (double));
-    memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double));
-    memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
+    memset((char*) D->stage_explicit_update, 0, D->number_of_elements * sizeof (double));
+    memset((char*) D->xmom_explicit_update, 0, D->number_of_elements * sizeof (double));
+    memset((char*) D->ymom_explicit_update, 0, D->number_of_elements * sizeof (double));
 
 
@@ -789 +741 @@
     RiverWall_count=0;
     // Which substep of the timestepping method are we on?
-    substep_count=(call-2)%timestep_fluxcalls;
+    substep_count=(call-2)%D->timestep_fluxcalls;
     
     // Fluxes are not updated every timestep,
     // but all fluxes ARE updated when the following condition holds
-    if(allow_timestep_increase[0]==1){
+    if(D->allow_timestep_increase[0]==1){
         // We can only increase the timestep if all fluxes are allowed to be updated
         // If this is not done the timestep can't increase (since local_timestep is static)
@@ -800 +752 @@
 
     // For all triangles
-    for (k = 0; k < number_of_elements; k++) {
+    for (k = 0; k < D->number_of_elements; k++) {
         speed_max_last=0.0;
 
@@ -809 +761 @@
             ki3 = 3*ki; 
 
-            if ((already_computed_flux[ki] == call) || (update_next_flux[ki]!=1)) {
+            if ((D->already_computed_flux[ki] == call) || (D->update_next_flux[ki]!=1)) {
                 // We've already computed the flux across this edge
                 // Check if it is a riverwall
-                if(edge_flux_type[ki]==1){
+                if(D->edge_flux_type[ki]==1){
                     // Update counter of riverwall edges == index of
                     // riverwall_elevation + riverwall_rowIndex
@@ -821 +773 @@
 
             // Get left hand side values from triangle k, edge i
-            ql[0] = stage_edge_values[ki];
-            ql[1] = xmom_edge_values[ki];
-            ql[2] = ymom_edge_values[ki];
-            zl = bed_edge_values[ki];
-            hc = height_centroid_values[k];
-            zc = bed_centroid_values[k];
-            hle= height_edge_values[ki];
+            ql[0] = D->stage_edge_values[ki];
+            ql[1] = D->xmom_edge_values[ki];
+            ql[2] = D->ymom_edge_values[ki];
+            zl = D->bed_edge_values[ki];
+            hc = D->height_centroid_values[k];
+            zc = D->bed_centroid_values[k];
+            hle= D->height_edge_values[ki];
 
             // Get right hand side values either from neighbouring triangle
             // or from boundary array (Quantities at neighbour on nearest face).
-            n = neighbours[ki];
+            n = D->neighbours[ki];
             hc_n = hc;
-            zc_n = bed_centroid_values[k];
+            zc_n = D->bed_centroid_values[k];
             if (n < 0) {
                 // Neighbour is a boundary condition
                 m = -n - 1; // Convert negative flag to boundary index
 
-                qr[0] = stage_boundary_values[m];
-                qr[1] = xmom_boundary_values[m];
-                qr[2] = ymom_boundary_values[m];
+                qr[0] = D->stage_boundary_values[m];
+                qr[1] = D->xmom_boundary_values[m];
+                qr[2] = D->ymom_boundary_values[m];
                 zr = zl; // Extend bed elevation to boundary
                 hre= max(qr[0]-zr,0.);//hle; 
             } else {
                 // Neighbour is a real triangle
-                hc_n = height_centroid_values[n];
-                zc_n = bed_centroid_values[n];
-                m = neighbour_edges[ki];
+                hc_n = D->height_centroid_values[n];
+                zc_n = D->bed_centroid_values[n];
+                m = D->neighbour_edges[ki];
                 nm = n * 3 + m; // Linear index (triangle n, edge m)
                 nm3 = nm*3;
 
-                qr[0] = stage_edge_values[nm];
-                qr[1] = xmom_edge_values[nm];
-                qr[2] = ymom_edge_values[nm];
-                zr = bed_edge_values[nm];
-                hre = height_edge_values[nm];
+                qr[0] = D->stage_edge_values[nm];
+                qr[1] = D->xmom_edge_values[nm];
+                qr[2] = D->ymom_edge_values[nm];
+                zr = D->bed_edge_values[nm];
+                hre = D->height_edge_values[nm];
             }
           
@@ -862 +814 @@
 
             //// Account for riverwalls
-            if(edge_flux_type[ki]==1){
+            if(D->edge_flux_type[ki]==1){
                 // Update counter of riverwall edges == index of
                 // riverwall_elevation + riverwall_rowIndex
@@ -868 +820 @@
                 
                 // Set central bed to riverwall elevation
-                z_half=max(riverwall_elevation[RiverWall_count-1], z_half) ;
+                z_half=max(D->riverwall_elevation[RiverWall_count-1], z_half) ;
 
             }
@@ -881 +833 @@
                     h_left, h_right,
                     hle, hre,
-                    normals[ki2], normals[ki2 + 1],
-                    epsilon, z_half, limiting_threshold, g,
-                    edgeflux, &max_speed, &pressure_flux, hc, hc_n);
+                    D->normals[ki2],D->normals[ki2 + 1],
+                    D->epsilon, z_half, limiting_threshold, D->g,
+                    edgeflux, &max_speed_local, &pressure_flux, hc, hc_n);
 
             // Force weir discharge to match weir theory
-            if(edge_flux_type[ki]==1){
-                weir_height=max(riverwall_elevation[RiverWall_count-1]-min(zl,zr), 0.); // Reference weir height  
+            if(D->edge_flux_type[ki]==1){
+                weir_height=max(D->riverwall_elevation[RiverWall_count-1]-min(zl,zr), 0.); // Reference weir height  
                 // If the weir height is zero, avoid the weir computation entirely
                 if(weir_height>0.){
@@ -893 +845 @@
                     // Use first-order h's for weir -- as the 'upstream/downstream' heads are
                     //  measured away from the weir itself
-                    h_left_tmp= max(stage_centroid_values[k]-z_half,0.);
+                    h_left_tmp= max(D->stage_centroid_values[k]-z_half,0.);
                     if(n>=0){
-                        h_right_tmp= max(stage_centroid_values[n]-z_half,0.);
+                        h_right_tmp= max(D->stage_centroid_values[n]-z_half,0.);
                     }else{
                         h_right_tmp= max(hc_n+zr-z_half,0.);
@@ -902 +854 @@
                     //////////////////////////////////////////////////////////////////////////////////
                     // Get Qfactor index - multiply the idealised weir discharge by this constant factor
-                    ii=riverwall_rowIndex[RiverWall_count-1]*ncol_riverwall_hydraulic_properties;
-                    Qfactor=riverwall_hydraulic_properties[ii];
+                    ii = D->riverwall_rowIndex[RiverWall_count-1] * D->ncol_riverwall_hydraulic_properties;
+                    Qfactor = D->riverwall_hydraulic_properties[ii];
                     // Get s1, submergence ratio at which we start blending with the shallow water solution 
                     ii+=1;
-                    s1=riverwall_hydraulic_properties[ii];
+                    s1= D->riverwall_hydraulic_properties[ii];
                     // Get s2, submergence ratio at which we entirely use the shallow water solution 
                     ii+=1;
-                    s2=riverwall_hydraulic_properties[ii];
+                    s2= D->riverwall_hydraulic_properties[ii];
                     // Get h1, tailwater head / weir height at which we start blending with the shallow water solution 
                     ii+=1;
-                    h1=riverwall_hydraulic_properties[ii];
+                    h1= D->riverwall_hydraulic_properties[ii];
                     // Get h2, tailwater head / weir height at which we entirely use the shallow water solution 
                     ii+=1;
-                    h2=riverwall_hydraulic_properties[ii];
+                    h2= D->riverwall_hydraulic_properties[ii];
                     
                     // Weir flux adjustment 
-                    adjust_edgeflux_with_weir(edgeflux, h_left_tmp, h_right_tmp, g, 
+                    adjust_edgeflux_with_weir(edgeflux, h_left_tmp, h_right_tmp, D->g, 
                                               weir_height, Qfactor, 
                                               s1, s2, h1, h2);
@@ -926 +878 @@
             
             // Multiply edgeflux by edgelength
-            length = edgelengths[ki];
+            length = D->edgelengths[ki];
             edgeflux[0] *= length;
             edgeflux[1] *= length;
@@ -938 +890 @@
             //    edgeflux[1] = 0.;
             //    edgeflux[2] = 0.;
-            //    //max_speed=0.;
+            //    //max_speed_local=0.;
             //    //pressure_flux=0.;
             //}
@@ -946 +898 @@
             //    edgeflux[1] = 0.;
             //    edgeflux[2] = 0.;
-            //    //max_speed=0.;
+            //    //max_speed_local=0.;
             //    //pressure_flux=0.;
             //}
 
-            edge_flux_work[ki3 + 0 ] = -edgeflux[0];
-            edge_flux_work[ki3 + 1 ] = -edgeflux[1];
-            edge_flux_work[ki3 + 2 ] = -edgeflux[2];
+            D->edge_flux_work[ki3 + 0 ] = -edgeflux[0];
+            D->edge_flux_work[ki3 + 1 ] = -edgeflux[1];
+            D->edge_flux_work[ki3 + 2 ] = -edgeflux[2];
 
             // bedslope_work contains all gravity related terms
-            bedslope_work=length*(-g*0.5*(h_left*h_left - hle*hle -(hle+hc)*(zl-zc))+pressure_flux);
-
-            pressuregrad_work[ki]=bedslope_work;
+            bedslope_work=length*(- D->g *0.5*(h_left*h_left - hle*hle -(hle+hc)*(zl-zc))+pressure_flux);
+
+            D->pressuregrad_work[ki]=bedslope_work;
             
-            already_computed_flux[ki] = call; // #k Done
+            D->already_computed_flux[ki] = call; // #k Done
 
             // Update neighbour n with same flux but reversed sign
             if (n >= 0) {
 
-                edge_flux_work[nm3 + 0 ] = edgeflux[0];
-                edge_flux_work[nm3 + 1 ] = edgeflux[1];
-                edge_flux_work[nm3 + 2 ] = edgeflux[2];
-                bedslope_work=length*(-g*0.5*(h_right*h_right-hre*hre-(hre+hc_n)*(zr-zc_n))+pressure_flux);
-                pressuregrad_work[nm]=bedslope_work;
-
-                already_computed_flux[nm] = call; // #n Done
+                D->edge_flux_work[nm3 + 0 ] = edgeflux[0];
+                D->edge_flux_work[nm3 + 1 ] = edgeflux[1];
+                D->edge_flux_work[nm3 + 2 ] = edgeflux[2];
+                bedslope_work=length*(- D->g *0.5*(h_right*h_right-hre*hre-(hre+hc_n)*(zr-zc_n))+pressure_flux);
+                D->pressuregrad_work[nm]=bedslope_work;
+
+                D->already_computed_flux[nm] = call; // #n Done
             }
 
@@ -979 +931 @@
 
                 // Compute the 'edge-timesteps' (useful for setting flux_update_frequency)
-                edge_timestep[ki] = radii[k] / max(max_speed, epsilon);
+                D->edge_timestep[ki] = D->radii[k] / max(max_speed_local, D->epsilon);
                 if (n >= 0) {
-                    edge_timestep[nm] = radii[n] / max(max_speed, epsilon);
+                    D->edge_timestep[nm] = D->radii[n] / max(max_speed_local, D->epsilon);
                 }
 
                 // Update the timestep
-                if ((tri_full_flag[k] == 1)) {
-
-                    speed_max_last=max(speed_max_last, max_speed);
-
-                    if (max_speed > epsilon) {
+                if ((D->tri_full_flag[k] == 1)) {
+
+                    speed_max_last=max(speed_max_last, max_speed_local);
+
+                    if (max_speed_local > D->epsilon) {
                         // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
 
                         // CFL for triangle k
-                        local_timestep = min(local_timestep, edge_timestep[ki]);
+                        local_timestep = min(local_timestep, D->edge_timestep[ki]);
 
                         if (n >= 0) {
                             // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
-                            local_timestep = min(local_timestep, edge_timestep[nm]);
+                            local_timestep = min(local_timestep, D->edge_timestep[nm]);
                         }
                     }
@@ -1005 +957 @@
         } // End edge i (and neighbour n)
         // Keep track of maximal speeds
-        if(substep_count==0) max_speed_array[k] = speed_max_last; //max_speed;
+        if(substep_count==0) D->max_speed[k] = speed_max_last; //max_speed;
 
 
@@ -1065 +1017 @@
 
     // Now add up stage, xmom, ymom explicit updates
-    for(k=0; k<number_of_elements; k++){
-        hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.);
+    for(k=0; k<D->number_of_elements; k++){
+        hc = max(D->stage_centroid_values[k] - D->bed_centroid_values[k],0.);
 
         for(i=0;i<3;i++){
@@ -1073 +1025 @@
             ki2=ki*2;
             ki3 = ki*3;
-            n=neighbours[ki];
+            n=D->neighbours[ki];
 
             // GD HACK
             // Option to limit advective fluxes
             //if(hc > H0){
-                stage_explicit_update[k] += edge_flux_work[ki3+0];
-                xmom_explicit_update[k] += edge_flux_work[ki3+1];
-                ymom_explicit_update[k] += edge_flux_work[ki3+2];
+                D->stage_explicit_update[k] += D->edge_flux_work[ki3+0];
+                D->xmom_explicit_update[k] += D->edge_flux_work[ki3+1];
+                D->ymom_explicit_update[k] += D->edge_flux_work[ki3+2];
             //}else{
             //    stage_explicit_update[k] += edge_flux_work[ki3+0];
@@ -1089 +1041 @@
             // condition OR a ghost cell, then add the flux to the
             // boundary_flux_integral
-            if( (n<0 & tri_full_flag[k]==1) | ( n>=0 && (tri_full_flag[k]==1 & tri_full_flag[n]==0)) ){ 
+            if( (n<0 & D->tri_full_flag[k]==1) | ( n>=0 && (D->tri_full_flag[k]==1 & D->tri_full_flag[n]==0)) ){ 
                 // boundary_flux_sum is an array with length = timestep_fluxcalls
                 // For each sub-step, we put the boundary flux sum in.
-                boundary_flux_sum[substep_count] += edge_flux_work[ki3];
+                D->boundary_flux_sum[substep_count] += D->edge_flux_work[ki3];
             }
     
@@ -1098 +1050 @@
             // Compute bed slope term
             //if(hc > H0){
-                xmom_explicit_update[k] -= normals[ki2]*pressuregrad_work[ki];
-                ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_work[ki];
+                D->xmom_explicit_update[k] -= D->normals[ki2]*D->pressuregrad_work[ki];
+                D->ymom_explicit_update[k] -= D->normals[ki2+1]*D->pressuregrad_work[ki];
             //}else{
             //    xmom_explicit_update[k] *= 0.;
@@ -1110 +1062 @@
         // Normalise triangle k by area and store for when all conserved
         // quantities get updated
-        inv_area = 1.0 / areas[k];
-        stage_explicit_update[k] *= inv_area;
-        xmom_explicit_update[k] *= inv_area;
-        ymom_explicit_update[k] *= inv_area;
+        inv_area = 1.0 / D->areas[k];
+        D->stage_explicit_update[k] *= inv_area;
+        D->xmom_explicit_update[k] *= inv_area;
+        D->ymom_explicit_update[k] *= inv_area;
    
     }  // end cell k
@@ -1242 +1194 @@
 
 // Computational routine
-int _extrapolate_second_order_edge_sw(int number_of_elements,
-                                 double epsilon,
-                                 double minimum_allowed_height,
-                                 double beta_w,
-                                 double beta_w_dry,
-                                 double beta_uh,
-                                 double beta_uh_dry,
-                                 double beta_vh,
-                                 double beta_vh_dry,
-                                 long* surrogate_neighbours,
-                                 long* neighbour_edges,
-                                 long* number_of_boundaries,
-                                 double* centroid_coordinates,
-                                 double* stage_centroid_values,
-                                 double* xmom_centroid_values,
-                                 double* ymom_centroid_values,
-                                 double* elevation_centroid_values,
-                                 double* height_centroid_values,
-                                 double* edge_coordinates,
-                                 double* stage_edge_values,
-                                 double* xmom_edge_values,
-                                 double* ymom_edge_values,
-                                 double* elevation_edge_values,
-                                 double* height_edge_values,
-                                 double* stage_vertex_values,
-                                 double* xmom_vertex_values,
-                                 double* ymom_vertex_values,
-                                 double* elevation_vertex_values,
-                                 double* height_vertex_values,
-                                 int optimise_dry_cells, 
-                                 int extrapolate_velocity_second_order,
-                                 double* x_centroid_work,
-                                 double* y_centroid_work,
-                                 long* update_extrapolation) {
+//int _extrapolate_second_order_edge_sw(int number_of_elements,
+//                                 double epsilon,
+//                                 double minimum_allowed_height,
+//                                 double beta_w,
+//                                 double beta_w_dry,
+//                                 double beta_uh,
+//                                 double beta_uh_dry,
+//                                 double beta_vh,
+//                                 double beta_vh_dry,
+//                                 long* surrogate_neighbours,
+//                                 long* neighbour_edges,
+//                                 long* number_of_boundaries,
+//                                 double* centroid_coordinates,
+//                                 double* stage_centroid_values,
+//                                 double* xmom_centroid_values,
+//                                 double* ymom_centroid_values,
+//                                 double* bed_centroid_values,
+//                                 double* height_centroid_values,
+//                                 double* edge_coordinates,
+//                                 double* stage_edge_values,
+//                                 double* xmom_edge_values,
+//                                 double* ymom_edge_values,
+//                                 double* bed_edge_values,
+//                                 double* height_edge_values,
+//                                 double* stage_vertex_values,
+//                                 double* xmom_vertex_values,
+//                                 double* ymom_vertex_values,
+//                                 double* bed_vertex_values,
+//                                 double* height_vertex_values,
+//                                 int optimise_dry_cells, 
+//                                 int extrapolate_velocity_second_order,
+//                                 double* x_centroid_work,
+//                                 double* y_centroid_work,
+//                                 long* update_extrapolation) {
+int _extrapolate_second_order_edge_sw(struct domain *D){
                   
   // Local variables
@@ -1287 +1240 @@
   
 
-  memset((char*) x_centroid_work, 0, number_of_elements * sizeof (double));
-  memset((char*) y_centroid_work, 0, number_of_elements * sizeof (double));
-  
-  
-  if(extrapolate_velocity_second_order==1){
+  memset((char*) D->x_centroid_work, 0, D->number_of_elements * sizeof (double));
+  memset((char*) D->y_centroid_work, 0, D->number_of_elements * sizeof (double));
+  
+  
+  if(D->extrapolate_velocity_second_order==1){
 
       // Replace momentum centroid with velocity centroid to allow velocity
       // extrapolation This will be changed back at the end of the routine
-      for (k=0; k<number_of_elements; k++){
+      for (k=0; k< D->number_of_elements; k++){
           
-          height_centroid_values[k] = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);
-
-          dk = height_centroid_values[k]; 
-          if(dk>minimum_allowed_height){
+          D->height_centroid_values[k] = max(D->stage_centroid_values[k] - D->bed_centroid_values[k], 0.);
+
+          dk = D->height_centroid_values[k]; 
+          if(dk> D->minimum_allowed_height){
               dk_inv=1.0/dk;
-              x_centroid_work[k] = xmom_centroid_values[k];
-              xmom_centroid_values[k] = xmom_centroid_values[k]*dk_inv;
-
-              y_centroid_work[k] = ymom_centroid_values[k];
-              ymom_centroid_values[k] = ymom_centroid_values[k]*dk_inv;
+              D->x_centroid_work[k] = D->xmom_centroid_values[k];
+              D->xmom_centroid_values[k] = D->xmom_centroid_values[k]*dk_inv;
+
+              D->y_centroid_work[k] = D->ymom_centroid_values[k];
+              D->ymom_centroid_values[k] = D->ymom_centroid_values[k]*dk_inv;
           }else{
-              x_centroid_work[k] = 0.;
-              xmom_centroid_values[k] = 0.;
-              y_centroid_work[k] = 0.;
-              ymom_centroid_values[k] = 0.;
+              D->x_centroid_work[k] = 0.;
+              D->xmom_centroid_values[k] = 0.;
+              D->y_centroid_work[k] = 0.;
+              D->ymom_centroid_values[k] = 0.;
 
          }
@@ -1321 +1274 @@
   // of water being trapped and unable to lose momentum, which can occur in
   // some situations
-  for (k=0; k<number_of_elements;k++){
+  for (k=0; k< D->number_of_elements;k++){
       
       k3=k*3;
-      k0 = surrogate_neighbours[k3];
-      k1 = surrogate_neighbours[k3 + 1];
-      k2 = surrogate_neighbours[k3 + 2];
-
-      if((height_centroid_values[k0] < minimum_allowed_height | k0==k) &
-         (height_centroid_values[k1] < minimum_allowed_height | k1==k) &
-         (height_centroid_values[k2] < minimum_allowed_height | k2==k)){
+      k0 = D->surrogate_neighbours[k3];
+      k1 = D->surrogate_neighbours[k3 + 1];
+      k2 = D->surrogate_neighbours[k3 + 2];
+
+      if((D->height_centroid_values[k0] < D->minimum_allowed_height | k0==k) &
+         (D->height_centroid_values[k1] < D->minimum_allowed_height | k1==k) &
+         (D->height_centroid_values[k2] < D->minimum_allowed_height | k2==k)){
                   //printf("Surrounded by dry cells\n");
-              x_centroid_work[k] = 0.;
-              xmom_centroid_values[k] = 0.;
-              y_centroid_work[k] = 0.;
-              ymom_centroid_values[k] = 0.;
+              D->x_centroid_work[k] = 0.;
+              D->xmom_centroid_values[k] = 0.;
+              D->y_centroid_work[k] = 0.;
+              D->ymom_centroid_values[k] = 0.;
 
       }
@@ -1343 +1296 @@
 
   // Begin extrapolation routine
-  for (k = 0; k < number_of_elements; k++) 
+  for (k = 0; k < D->number_of_elements; k++) 
   {
 
     // Don't update the extrapolation if the flux will not be computed on the
     // next timestep
-    if(update_extrapolation[k]==0){
+    if(D->update_extrapolation[k]==0){
        continue;
     }
@@ -1357 +1310 @@
     k6=k*6;
 
-    if (number_of_boundaries[k]==3)
+    if (D->number_of_boundaries[k]==3)
     {
       // No neighbours, set gradient on the triangle to zero
       
-      stage_edge_values[k3]   = stage_centroid_values[k];
-      stage_edge_values[k3+1] = stage_centroid_values[k];
-      stage_edge_values[k3+2] = stage_centroid_values[k];
+      D->stage_edge_values[k3]   = D->stage_centroid_values[k];
+      D->stage_edge_values[k3+1] = D->stage_centroid_values[k];
+      D->stage_edge_values[k3+2] = D->stage_centroid_values[k];
 
       //xmom_centroid_values[k] = 0.;
       //ymom_centroid_values[k] = 0.;
       
-      xmom_edge_values[k3]    = xmom_centroid_values[k];
-      xmom_edge_values[k3+1]  = xmom_centroid_values[k];
-      xmom_edge_values[k3+2]  = xmom_centroid_values[k];
-      ymom_edge_values[k3]    = ymom_centroid_values[k];
-      ymom_edge_values[k3+1]  = ymom_centroid_values[k];
-      ymom_edge_values[k3+2]  = ymom_centroid_values[k];
-
-      dk = height_centroid_values[k];
-      height_edge_values[k3] = dk;
-      height_edge_values[k3+1] = dk;
-      height_edge_values[k3+2] = dk;
+      D->xmom_edge_values[k3]    = D->xmom_centroid_values[k];
+      D->xmom_edge_values[k3+1]  = D->xmom_centroid_values[k];
+      D->xmom_edge_values[k3+2]  = D->xmom_centroid_values[k];
+      D->ymom_edge_values[k3]    = D->ymom_centroid_values[k];
+      D->ymom_edge_values[k3+1]  = D->ymom_centroid_values[k];
+      D->ymom_edge_values[k3+2]  = D->ymom_centroid_values[k];
+
+      dk = D->height_centroid_values[k];
+      D->height_edge_values[k3] = dk;
+      D->height_edge_values[k3+1] = dk;
+      D->height_edge_values[k3+2] = dk;
       
       continue;
@@ -1388 +1341 @@
       
       // Get the edge coordinates
-      xv0 = edge_coordinates[k6];   
-      yv0 = edge_coordinates[k6+1];
-      xv1 = edge_coordinates[k6+2]; 
-      yv1 = edge_coordinates[k6+3];
-      xv2 = edge_coordinates[k6+4]; 
-      yv2 = edge_coordinates[k6+5];
+      xv0 = D->edge_coordinates[k6];   
+      yv0 = D->edge_coordinates[k6+1];
+      xv1 = D->edge_coordinates[k6+2]; 
+      yv1 = D->edge_coordinates[k6+3];
+      xv2 = D->edge_coordinates[k6+4]; 
+      yv2 = D->edge_coordinates[k6+5];
       
       // Get the centroid coordinates
       coord_index = 2*k;
-      x = centroid_coordinates[coord_index];
-      y = centroid_coordinates[coord_index+1];
+      x = D->centroid_coordinates[coord_index];
+      y = D->centroid_coordinates[coord_index+1];
       
       // Store x- and y- differentials for the edges of 
@@ -1413 +1366 @@
 
             
-    if (number_of_boundaries[k]<=1)
+    if (D->number_of_boundaries[k]<=1)
     {
       //==============================================
@@ -1426 +1379 @@
       // from this centroid and its two neighbours
       
-      k0 = surrogate_neighbours[k3];
-      k1 = surrogate_neighbours[k3 + 1];
-      k2 = surrogate_neighbours[k3 + 2];
+      k0 = D->surrogate_neighbours[k3];
+      k1 = D->surrogate_neighbours[k3 + 1];
+      k2 = D->surrogate_neighbours[k3 + 2];
 
 
@@ -1434 +1387 @@
       // (really the centroids of neighbouring triangles)
       coord_index = 2*k0;
-      x0 = centroid_coordinates[coord_index];
-      y0 = centroid_coordinates[coord_index+1];
+      x0 = D->centroid_coordinates[coord_index];
+      y0 = D->centroid_coordinates[coord_index+1];
       
       coord_index = 2*k1;
-      x1 = centroid_coordinates[coord_index];
-      y1 = centroid_coordinates[coord_index+1];
+      x1 = D->centroid_coordinates[coord_index];
+      y1 = D->centroid_coordinates[coord_index+1];
       
       coord_index = 2*k2;
-      x2 = centroid_coordinates[coord_index];
-      y2 = centroid_coordinates[coord_index+1];
+      x2 = D->centroid_coordinates[coord_index];
+      y2 = D->centroid_coordinates[coord_index+1];
 
       // Store x- and y- differentials for the vertices 
@@ -1463 +1416 @@
 
           // Isolated wet cell -- constant stage/depth extrapolation 
-          stage_edge_values[k3]   = stage_centroid_values[k];
-          stage_edge_values[k3+1] = stage_centroid_values[k];
-          stage_edge_values[k3+2] = stage_centroid_values[k];
-
-          dk=height_centroid_values[k]; //max(stage_centroid_values[k]-elevation_centroid_values[k],0.);
-          height_edge_values[k3] = dk;
-          height_edge_values[k3+1] = dk;
-          height_edge_values[k3+2] = dk;
+          D->stage_edge_values[k3]   = D->stage_centroid_values[k];
+          D->stage_edge_values[k3+1] = D->stage_centroid_values[k];
+          D->stage_edge_values[k3+2] = D->stage_centroid_values[k];
+
+          dk= D->height_centroid_values[k]; //max(stage_centroid_values[k]-bed_centroid_values[k],0.);
+          D->height_edge_values[k3] = dk;
+          D->height_edge_values[k3+1] = dk;
+          D->height_edge_values[k3+2] = dk;
           
-          xmom_edge_values[k3]    = xmom_centroid_values[k];
-          xmom_edge_values[k3+1]  = xmom_centroid_values[k];
-          xmom_edge_values[k3+2]  = xmom_centroid_values[k];
-          ymom_edge_values[k3]    = ymom_centroid_values[k];
-          ymom_edge_values[k3+1]  = ymom_centroid_values[k];
-          ymom_edge_values[k3+2]  = ymom_centroid_values[k];
+          D->xmom_edge_values[k3]    = D->xmom_centroid_values[k];
+          D->xmom_edge_values[k3+1]  = D->xmom_centroid_values[k];
+          D->xmom_edge_values[k3+2]  = D->xmom_centroid_values[k];
+          D->ymom_edge_values[k3]    = D->ymom_centroid_values[k];
+          D->ymom_edge_values[k3+1]  = D->ymom_centroid_values[k];
+          D->ymom_edge_values[k3+2]  = D->ymom_centroid_values[k];
 
           continue;
@@ -1483 +1436 @@
       
       // Calculate heights of neighbouring cells
-      hc = height_centroid_values[k];
-      h0 = height_centroid_values[k0];
-      h1 = height_centroid_values[k1];
-      h2 = height_centroid_values[k2];
+      hc = D->height_centroid_values[k];
+      h0 = D->height_centroid_values[k0];
+      h1 = D->height_centroid_values[k1];
+      h2 = D->height_centroid_values[k2];
       
       hmin = min(min(h0, min(h1, h2)), hc);
@@ -1507 +1460 @@
       // Set hfactor to zero smothly as hmin--> minimum_allowed_height. This
       // avoids some 'chatter' for very shallow flows 
-      hfactor=min( 1.2*max(hmin-minimum_allowed_height,0.)/(max(hmin,0.)+1.*minimum_allowed_height), hfactor);
+      hfactor=min( 1.2*max(hmin- D->minimum_allowed_height,0.)/(max(hmin,0.)+1.* D->minimum_allowed_height), hfactor);
 
       inv_area2 = 1.0/area2; 
@@ -1514 +1467 @@
       //-----------------------------------      
       
-      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
+      beta_tmp = D->beta_w_dry + (D->beta_w - D->beta_w_dry) * hfactor;
      
       if(beta_tmp>0.){ 
           // Calculate the difference between vertex 0 of the auxiliary 
           // triangle and the centroid of triangle k
-          dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
+          dq0 = D->stage_centroid_values[k0] - D->stage_centroid_values[k];
           
           // Calculate differentials between the vertices 
           // of the auxiliary triangle (centroids of neighbouring triangles)
-          dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
-          dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
+          dq1 = D->stage_centroid_values[k1] - D->stage_centroid_values[k0];
+          dq2 = D->stage_centroid_values[k2] - D->stage_centroid_values[k0];
          
           // Calculate the gradient of stage on the auxiliary triangle
@@ -1545 +1498 @@
           limit_gradient(dqv, qmin, qmax, beta_tmp);
 
-          stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
-          stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
-          stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
+          D->stage_edge_values[k3+0] = D->stage_centroid_values[k] + dqv[0];
+          D->stage_edge_values[k3+1] = D->stage_centroid_values[k] + dqv[1];
+          D->stage_edge_values[k3+2] = D->stage_centroid_values[k] + dqv[2];
       }else{
           // Fast alternative when beta_tmp==0
-          stage_edge_values[k3+0] = stage_centroid_values[k];
-          stage_edge_values[k3+1] = stage_centroid_values[k];
-          stage_edge_values[k3+2] = stage_centroid_values[k];
+          D->stage_edge_values[k3+0] = D->stage_centroid_values[k];
+          D->stage_edge_values[k3+1] = D->stage_centroid_values[k];
+          D->stage_edge_values[k3+2] = D->stage_centroid_values[k];
       }
 
@@ -1563 +1516 @@
           // Calculate the difference between vertex 0 of the auxiliary 
           // triangle and the centroid of triangle k
-          dq0 = height_centroid_values[k0] - height_centroid_values[k];
+          dq0 = D->height_centroid_values[k0] - D->height_centroid_values[k];
           
           // Calculate differentials between the vertices 
           // of the auxiliary triangle (centroids of neighbouring triangles)
-          dq1 = height_centroid_values[k1] - height_centroid_values[k0];
-          dq2 = height_centroid_values[k2] - height_centroid_values[k0];
+          dq1 = D->height_centroid_values[k1] - D->height_centroid_values[k0];
+          dq2 = D->height_centroid_values[k2] - D->height_centroid_values[k0];
          
           // Calculate the gradient of height on the auxiliary triangle
@@ -1593 +1546 @@
           //beta_tmp = 0. + (beta_w - 0.) * hfactor;
 
-          height_edge_values[k3+0] = height_centroid_values[k] + dqv[0];
-          height_edge_values[k3+1] = height_centroid_values[k] + dqv[1];
-          height_edge_values[k3+2] = height_centroid_values[k] + dqv[2];
+          D->height_edge_values[k3+0] = D->height_centroid_values[k] + dqv[0];
+          D->height_edge_values[k3+1] = D->height_centroid_values[k] + dqv[1];
+          D->height_edge_values[k3+2] = D->height_centroid_values[k] + dqv[2];
       }else{
           // Fast alternative when beta_tmp==0
-          height_edge_values[k3+0] = height_centroid_values[k];
-          height_edge_values[k3+1] = height_centroid_values[k];
-          height_edge_values[k3+2] = height_centroid_values[k];
+          D->height_edge_values[k3+0] = D->height_centroid_values[k];
+          D->height_edge_values[k3+1] = D->height_centroid_values[k];
+          D->height_edge_values[k3+2] = D->height_centroid_values[k];
       }
       //-----------------------------------
@@ -1606 +1559 @@
       //-----------------------------------            
 
-      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
+      beta_tmp = D->beta_uh_dry + (D->beta_uh - D->beta_uh_dry) * hfactor;
       if(beta_tmp>0.){
           // Calculate the difference between vertex 0 of the auxiliary 
           // triangle and the centroid of triangle k      
-          dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
+          dq0 = D->xmom_centroid_values[k0] - D->xmom_centroid_values[k];
           
           // Calculate differentials between the vertices 
           // of the auxiliary triangle
-          dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
-          dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
+          dq1 = D->xmom_centroid_values[k1] - D->xmom_centroid_values[k0];
+          dq2 = D->xmom_centroid_values[k2] - D->xmom_centroid_values[k0];
           
           // Calculate the gradient of xmom on the auxiliary triangle
@@ -1641 +1594 @@
           for (i=0; i < 3; i++)
           {
-            xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
+            D->xmom_edge_values[k3+i] = D->xmom_centroid_values[k] + dqv[i];
           }
       }else{
@@ -1647 +1600 @@
           for (i=0; i < 3; i++)
           {
-            xmom_edge_values[k3+i] = xmom_centroid_values[k];
+            D->xmom_edge_values[k3+i] = D->xmom_centroid_values[k];
           }
       }
@@ -1655 +1608 @@
       //-----------------------------------                  
       
-      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
+      beta_tmp = D->beta_vh_dry + (D->beta_vh - D->beta_vh_dry) * hfactor;   
 
       if(beta_tmp>0.){
           // Calculate the difference between vertex 0 of the auxiliary 
           // triangle and the centroid of triangle k
-          dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
+          dq0 = D->ymom_centroid_values[k0] - D->ymom_centroid_values[k];
           
           // Calculate differentials between the vertices 
           // of the auxiliary triangle
-          dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
-          dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
+          dq1 = D->ymom_centroid_values[k1] - D->ymom_centroid_values[k0];
+          dq2 = D->ymom_centroid_values[k2] - D->ymom_centroid_values[k0];
           
           // Calculate the gradient of xmom on the auxiliary triangle
@@ -1691 +1644 @@
           for (i=0;i<3;i++)
           {
-            ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+            D->ymom_edge_values[k3 + i] = D->ymom_centroid_values[k] + dqv[i];
           }
       }else{
@@ -1697 +1650 @@
           for (i=0;i<3;i++)
           {
-            ymom_edge_values[k3 + i] = ymom_centroid_values[k];
+            D->ymom_edge_values[k3 + i] = D->ymom_centroid_values[k];
           }
 
@@ -1718 +1671 @@
       // k2 indexes the edges of triangle k   
       
-          if (surrogate_neighbours[k2] != k)
+          if (D->surrogate_neighbours[k2] != k)
           {
              break;
@@ -1731 +1684 @@
       }
       
-      k1 = surrogate_neighbours[k2];
+      k1 = D->surrogate_neighbours[k2];
       
       // The coordinates of the triangle are already (x,y). 
       // Get centroid of the neighbour (x1,y1)
       coord_index = 2*k1;
-      x1 = centroid_coordinates[coord_index];
-      y1 = centroid_coordinates[coord_index + 1];
+      x1 = D->centroid_coordinates[coord_index];
+      y1 = D->centroid_coordinates[coord_index + 1];
       
       // Compute x- and y- distances between the centroid of 
@@ -1761 +1714 @@
 
       // Compute differentials
-      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
+      dq1 = D->stage_centroid_values[k1] - D->stage_centroid_values[k];
       
       // Calculate the gradient between the centroid of triangle k 
@@ -1786 +1739 @@
       
       // Limit the gradient
-      limit_gradient(dqv, qmin, qmax, beta_w);
-      
-      stage_edge_values[k3] = stage_centroid_values[k] + dqv[0];
-      stage_edge_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
-      stage_edge_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
+      limit_gradient(dqv, qmin, qmax, D->beta_w);
+      
+      D->stage_edge_values[k3] = D->stage_centroid_values[k] + dqv[0];
+      D->stage_edge_values[k3 + 1] = D->stage_centroid_values[k] + dqv[1];
+      D->stage_edge_values[k3 + 2] = D->stage_centroid_values[k] + dqv[2];
 
       //-----------------------------------
@@ -1797 +1750 @@
       
       // Compute differentials
-      dq1 = height_centroid_values[k1] - height_centroid_values[k];
+      dq1 = D->height_centroid_values[k1] - D->height_centroid_values[k];
       
       // Calculate the gradient between the centroid of triangle k 
@@ -1822 +1775 @@
       
       // Limit the gradient
-      limit_gradient(dqv, qmin, qmax, beta_w);
-      
-      height_edge_values[k3] = height_centroid_values[k] + dqv[0];
-      height_edge_values[k3 + 1] = height_centroid_values[k] + dqv[1];
-      height_edge_values[k3 + 2] = height_centroid_values[k] + dqv[2];
+      limit_gradient(dqv, qmin, qmax, D->beta_w);
+      
+      D->height_edge_values[k3] = D->height_centroid_values[k] + dqv[0];
+      D->height_edge_values[k3 + 1] = D->height_centroid_values[k] + dqv[1];
+      D->height_edge_values[k3 + 2] = D->height_centroid_values[k] + dqv[2];
 
       //-----------------------------------
@@ -1833 +1786 @@
       
       // Compute differentials
-      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
+      dq1 = D->xmom_centroid_values[k1] - D->xmom_centroid_values[k];
       
       // Calculate the gradient between the centroid of triangle k 
@@ -1858 +1811 @@
       
       // Limit the gradient      
-      limit_gradient(dqv, qmin, qmax, beta_w);
+      limit_gradient(dqv, qmin, qmax, D->beta_w);
       
       for (i = 0; i < 3;i++)
       {
-          xmom_edge_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
+          D->xmom_edge_values[k3 + i] = D->xmom_centroid_values[k] + dqv[i];
       }
       
@@ -1870 +1823 @@
 
       // Compute differentials
-      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
+      dq1 = D->ymom_centroid_values[k1] - D->ymom_centroid_values[k];
       
       // Calculate the gradient between the centroid of triangle k 
@@ -1895 +1848 @@
       
       // Limit the gradient
-      limit_gradient(dqv, qmin, qmax, beta_w);
+      limit_gradient(dqv, qmin, qmax, D->beta_w);
       
       for (i=0;i<3;i++)
               {
-              ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+              D->ymom_edge_values[k3 + i] = D->ymom_centroid_values[k] + dqv[i];
               }
     } // else [number_of_boundaries==2]
@@ -1906 +1859 @@
 
   // Compute vertex values of quantities
-  for (k=0; k<number_of_elements; k++){
-      if(extrapolate_velocity_second_order==1){
+  for (k=0; k< D->number_of_elements; k++){
+      if(D->extrapolate_velocity_second_order==1){
           //Convert velocity back to momenta at centroids
-          xmom_centroid_values[k] = x_centroid_work[k];
-          ymom_centroid_values[k] = y_centroid_work[k];
+          D->xmom_centroid_values[k] = D->x_centroid_work[k];
+          D->ymom_centroid_values[k] = D->y_centroid_work[k];
       }
      
       // Don't proceed if we didn't update the edge/vertex values
-      if(update_extrapolation[k]==0){
+      if(D->update_extrapolation[k]==0){
          continue;
       }
@@ -1921 +1874 @@
       
       // Compute stage vertex values 
-      stage_vertex_values[k3] = stage_edge_values[k3+1] + stage_edge_values[k3+2] -stage_edge_values[k3] ; 
-      stage_vertex_values[k3+1] =  stage_edge_values[k3] + stage_edge_values[k3+2]-stage_edge_values[k3+1]; 
-      stage_vertex_values[k3+2] =  stage_edge_values[k3] + stage_edge_values[k3+1]-stage_edge_values[k3+2]; 
+      D->stage_vertex_values[k3] = D->stage_edge_values[k3+1] + D->stage_edge_values[k3+2] - D->stage_edge_values[k3] ; 
+      D->stage_vertex_values[k3+1] =  D->stage_edge_values[k3] + D->stage_edge_values[k3+2]- D->stage_edge_values[k3+1]; 
+      D->stage_vertex_values[k3+2] =  D->stage_edge_values[k3] + D->stage_edge_values[k3+1]- D->stage_edge_values[k3+2]; 
       
       // Compute height vertex values 
-      height_vertex_values[k3] = height_edge_values[k3+1] + height_edge_values[k3+2] -height_edge_values[k3] ; 
-      height_vertex_values[k3+1] =  height_edge_values[k3] + height_edge_values[k3+2]-height_edge_values[k3+1]; 
-      height_vertex_values[k3+2] =  height_edge_values[k3] + height_edge_values[k3+1]-height_edge_values[k3+2]; 
+      D->height_vertex_values[k3] = D->height_edge_values[k3+1] + D->height_edge_values[k3+2] - D->height_edge_values[k3] ; 
+      D->height_vertex_values[k3+1] =  D->height_edge_values[k3] + D->height_edge_values[k3+2]- D->height_edge_values[k3+1]; 
+      D->height_vertex_values[k3+2] =  D->height_edge_values[k3] + D->height_edge_values[k3+1]- D->height_edge_values[k3+2]; 
 
       // If needed, convert from velocity to momenta
-      if(extrapolate_velocity_second_order==1){
+      if(D->extrapolate_velocity_second_order==1){
           // Re-compute momenta at edges
           for (i=0; i<3; i++){
-              dk= height_edge_values[k3+i]; 
-              xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*dk;
-              ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*dk;
+              dk= D->height_edge_values[k3+i]; 
+              D->xmom_edge_values[k3+i] = D->xmom_edge_values[k3+i]*dk;
+              D->ymom_edge_values[k3+i] = D->ymom_edge_values[k3+i]*dk;
           }
       }
       // Compute momenta at vertices
-      xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ; 
-      xmom_vertex_values[k3+1] =  xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1]; 
-      xmom_vertex_values[k3+2] =  xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2]; 
-      ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ; 
-      ymom_vertex_values[k3+1] =  ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1]; 
-      ymom_vertex_values[k3+2] =  ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2]; 
+      D->xmom_vertex_values[k3]   =  D->xmom_edge_values[k3+1] + D->xmom_edge_values[k3+2] - D->xmom_edge_values[k3] ; 
+      D->xmom_vertex_values[k3+1] =  D->xmom_edge_values[k3] + D->xmom_edge_values[k3+2]- D->xmom_edge_values[k3+1]; 
+      D->xmom_vertex_values[k3+2] =  D->xmom_edge_values[k3] + D->xmom_edge_values[k3+1]- D->xmom_edge_values[k3+2]; 
+      D->ymom_vertex_values[k3]   =  D->ymom_edge_values[k3+1] + D->ymom_edge_values[k3+2] - D->ymom_edge_values[k3] ; 
+      D->ymom_vertex_values[k3+1] =  D->ymom_edge_values[k3] + D->ymom_edge_values[k3+2]- D->ymom_edge_values[k3+1]; 
+      D->ymom_vertex_values[k3+2] =  D->ymom_edge_values[k3] + D->ymom_edge_values[k3+1]- D->ymom_edge_values[k3+2]; 
 
       // Compute new bed elevation
-      elevation_edge_values[k3]=stage_edge_values[k3]-height_edge_values[k3];
-      elevation_edge_values[k3+1]=stage_edge_values[k3+1]-height_edge_values[k3+1];
-      elevation_edge_values[k3+2]=stage_edge_values[k3+2]-height_edge_values[k3+2];
-      elevation_vertex_values[k3] = elevation_edge_values[k3+1] + elevation_edge_values[k3+2] -elevation_edge_values[k3] ; 
-      elevation_vertex_values[k3+1] =  elevation_edge_values[k3] + elevation_edge_values[k3+2]-elevation_edge_values[k3+1]; 
-      elevation_vertex_values[k3+2] =  elevation_edge_values[k3] + elevation_edge_values[k3+1]-elevation_edge_values[k3+2]; 
+      D->bed_edge_values[k3]= D->stage_edge_values[k3]- D->height_edge_values[k3];
+      D->bed_edge_values[k3+1]= D->stage_edge_values[k3+1]- D->height_edge_values[k3+1];
+      D->bed_edge_values[k3+2]= D->stage_edge_values[k3+2]- D->height_edge_values[k3+2];
+      D->bed_vertex_values[k3] = D->bed_edge_values[k3+1] + D->bed_edge_values[k3+2] - D->bed_edge_values[k3] ; 
+      D->bed_vertex_values[k3+1] =  D->bed_edge_values[k3] + D->bed_edge_values[k3+2] - D->bed_edge_values[k3+1]; 
+      D->bed_vertex_values[k3+2] =  D->bed_edge_values[k3] + D->bed_edge_values[k3+1] - D->bed_edge_values[k3+2]; 
   } 
 
@@ -1984 +1937 @@
     is converted to a timestep that must not be exceeded. The minimum of
     those is computed as the next overall timestep.
-
-    Python call:
-    domain.timestep = compute_fluxes(timestep,
-                                     domain.epsilon,
-                                     domain.H0,
-                                     domain.g,
-                                     domain.neighbours,
-                                     domain.neighbour_edges,
-                                     domain.normals,
-                                     domain.edgelengths,
-                                     domain.radii,
-                                     domain.areas,
-                                     tri_full_flag,
-                                     Stage.edge_values,
-                                     Xmom.edge_values,
-                                     Ymom.edge_values,
-                                     Bed.edge_values,
-                                     Stage.boundary_values,
-                                     Xmom.boundary_values,
-                                     Ymom.boundary_values,
-                                     Stage.explicit_update,
-                                     Xmom.explicit_update,
-                                     Ymom.explicit_update,
-                                     already_computed_flux,
-                                     optimise_dry_cells,
-                                     stage.centroid_values,
-                                     bed.centroid_values,
-                                     domain.riverwall_elevation)
-
-
-    Post conditions:
-      domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
-
-
   */
-
-
-  PyArrayObject *boundary_flux_sum, *neighbours, *neighbour_edges,
-    *normals, *edgelengths, *radii, *areas,
-    *tri_full_flag,
-    *stage_edge_values,
-    *xmom_edge_values,
-    *ymom_edge_values,
-    *bed_edge_values,
-    *height_edge_values,
-    *stage_boundary_values,
-    *xmom_boundary_values,
-    *ymom_boundary_values,
-    *boundary_flux_type,
-    *stage_explicit_update,
-    *xmom_explicit_update,
-    *ymom_explicit_update,
-    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
-    *max_speed_array, //Keeps track of max speeds for each triangle
-    *stage_centroid_values,
-    *xmom_centroid_values,
-    *ymom_centroid_values,
-    *bed_centroid_values,
-    *height_centroid_values,
-    *bed_vertex_values,
-    *edge_flux_type,
-    *riverwall_elevation,
-    *riverwall_rowIndex,
-    *riverwall_hydraulic_properties,
-    *edge_flux_work,
-    *pressuregrad_work,
-    *edge_timestep,
-    *update_next_flux,
-    *allow_timestep_increase;
+  struct domain D;
+  PyObject *domain;
+
+   
+  double timestep;
+  
+  if (!PyArg_ParseTuple(args, "Od", &domain, &timestep)) {
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
+  }
     
-  double timestep, epsilon, H0, g;
-  int optimise_dry_cells, timestep_fluxcalls, ncol_riverwall_hydraulic_properties;
-  
-    
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOOiiOOOOOOOOiOOOOOO",
-            &timestep,
-            &epsilon,
-            &H0,
-            &g,
-            &boundary_flux_sum,
-            &neighbours,
-            &neighbour_edges,
-            &normals,
-            &edgelengths, &radii, &areas,
-            &tri_full_flag,
-            &stage_edge_values,
-            &xmom_edge_values,
-            &ymom_edge_values,
-            &bed_edge_values,
-            &height_edge_values,
-            &stage_boundary_values,
-            &xmom_boundary_values,
-            &ymom_boundary_values,
-            &boundary_flux_type,
-            &stage_explicit_update,
-            &xmom_explicit_update,
-            &ymom_explicit_update,
-            &already_computed_flux,
-            &max_speed_array,
-            &optimise_dry_cells,
-            &timestep_fluxcalls,
-            &stage_centroid_values,
-            &xmom_centroid_values,
-            &ymom_centroid_values,
-            &bed_centroid_values,
-            &height_centroid_values,
-            //&bed_vertex_values,
-            &edge_flux_type,
-            &riverwall_elevation,
-            &riverwall_rowIndex,
-            &ncol_riverwall_hydraulic_properties,
-            &riverwall_hydraulic_properties,
-            &edge_flux_work,
-            &pressuregrad_work,
-            &edge_timestep,
-            &update_next_flux,
-            &allow_timestep_increase
-            )) {
-    report_python_error(AT, "could not parse input arguments");
-    return NULL;
-  }
-
-  // check that numpy array objects arrays are C contiguous memory
-  CHECK_C_CONTIG(neighbours);
-  CHECK_C_CONTIG(neighbour_edges);
-  CHECK_C_CONTIG(normals);
-  CHECK_C_CONTIG(edgelengths);
-  CHECK_C_CONTIG(radii);
-  CHECK_C_CONTIG(areas);
-  CHECK_C_CONTIG(tri_full_flag);
-  CHECK_C_CONTIG(stage_edge_values);
-  CHECK_C_CONTIG(xmom_edge_values);
-  CHECK_C_CONTIG(ymom_edge_values);
-  CHECK_C_CONTIG(bed_edge_values);
-  CHECK_C_CONTIG(height_edge_values);
-  CHECK_C_CONTIG(stage_boundary_values);
-  CHECK_C_CONTIG(xmom_boundary_values);
-  CHECK_C_CONTIG(ymom_boundary_values);
-  CHECK_C_CONTIG(boundary_flux_type);
-  CHECK_C_CONTIG(stage_explicit_update);
-  CHECK_C_CONTIG(xmom_explicit_update);
-  CHECK_C_CONTIG(ymom_explicit_update);
-  CHECK_C_CONTIG(already_computed_flux);
-  CHECK_C_CONTIG(max_speed_array);
-  CHECK_C_CONTIG(stage_centroid_values);
-  CHECK_C_CONTIG(xmom_centroid_values);
-  CHECK_C_CONTIG(ymom_centroid_values);
-  CHECK_C_CONTIG(bed_centroid_values);
-  CHECK_C_CONTIG(height_centroid_values);
-  //CHECK_C_CONTIG(bed_vertex_values);
-  CHECK_C_CONTIG(edge_flux_type);
-  CHECK_C_CONTIG(riverwall_elevation);
-  CHECK_C_CONTIG(riverwall_rowIndex);
-  CHECK_C_CONTIG(riverwall_hydraulic_properties);
-  CHECK_C_CONTIG(edge_flux_work);
-  CHECK_C_CONTIG(pressuregrad_work);
-  CHECK_C_CONTIG(edge_timestep);
-  CHECK_C_CONTIG(update_next_flux);
-  CHECK_C_CONTIG(allow_timestep_increase);
-
-  int number_of_elements = stage_edge_values -> dimensions[0];
-
-  // Call underlying flux computation routine and update 
-  // the explicit update arrays 
-  timestep = _compute_fluxes_central(number_of_elements,
-                     timestep,
-                     epsilon,
-                     H0,
-                     g,
-                     (double*) boundary_flux_sum -> data,
-                     (long*) neighbours -> data,
-                     (long*) neighbour_edges -> data,
-                     (double*) normals -> data,
-                     (double*) edgelengths -> data, 
-                     (double*) radii -> data, 
-                     (double*) areas -> data,
-                     (long*) tri_full_flag -> data,
-                     (double*) stage_edge_values -> data,
-                     (double*) xmom_edge_values -> data,
-                     (double*) ymom_edge_values -> data,
-                     (double*) bed_edge_values -> data,
-                     (double*) height_edge_values -> data,
-                     (double*) stage_boundary_values -> data,
-                     (double*) xmom_boundary_values -> data,
-                     (double*) ymom_boundary_values -> data,
-                     (long*)   boundary_flux_type -> data,
-                     (double*) stage_explicit_update -> data,
-                     (double*) xmom_explicit_update -> data,
-                     (double*) ymom_explicit_update -> data,
-                     (long*) already_computed_flux -> data,
-                     (double*) max_speed_array -> data,
-                     optimise_dry_cells, 
-                     timestep_fluxcalls,
-                     (double*) stage_centroid_values -> data, 
-                     (double*) xmom_centroid_values -> data, 
-                     (double*) ymom_centroid_values -> data, 
-                     (double*) bed_centroid_values -> data,
-                     (double*) height_centroid_values -> data,
-                     //(double*) bed_vertex_values -> data,
-                     (long*)   edge_flux_type-> data,
-                     (double*) riverwall_elevation-> data,
-                     (long*)   riverwall_rowIndex-> data,
-                     ncol_riverwall_hydraulic_properties,
-                     (double*) riverwall_hydraulic_properties ->data, 
-                     (double*) edge_flux_work-> data,
-                     (double*) pressuregrad_work->data,
-                     (double*) edge_timestep->data,
-                     (long*) update_next_flux->data,
-                     (long*) allow_timestep_increase->data
-                     );
+  get_python_domain(&D,domain);
+
+  timestep=_compute_fluxes_central(&D,timestep);
+
   // Return updated flux timestep
   return Py_BuildValue("d", timestep);
@@ -2324 +2081 @@
 
 PyObject *compute_flux_update_frequency(PyObject *self, PyObject *args) {
-  /*Compute how often we should update fluxes and extrapolations (not every timestep)
-
-    python_call
-        compute_flux_update_frequency_ext(
-                self.timestep,
-                self.neighbours,
-                self.neighbour_edges,
-                self.already_computed_flux,
-                self.edge_timestep,
-                self.flux_update_frequency,
-                self.update_extrapolation,
-                self.max_flux_update_frequency,
-                self.update_next_flux)
-
-
+  /*
+    
+    Compute how often we should update fluxes and extrapolations (perhaps not every timestep)
 
   */
 
-
-  PyArrayObject *neighbours, *neighbour_edges,
-    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
-    *edge_timestep,
-    *flux_update_frequency,
-    *update_extrapolation,
-    *update_next_flux,
-    *allow_timestep_increase;
+  struct domain D;
+  PyObject *domain;
+  
     
   double timestep;
   int max_flux_update_frequency;
   
+  if (!PyArg_ParseTuple(args, "Od", &domain, &timestep)) {
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
+  }
     
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dOOOOOOiOO",
-            &timestep,
-            &neighbours,
-            &neighbour_edges,
-            &already_computed_flux,
-            &edge_timestep,
-            &flux_update_frequency,
-            &update_extrapolation,
-            &max_flux_update_frequency,
-            &update_next_flux,
-            &allow_timestep_increase
-            )) {
-    report_python_error(AT, "could not parse input arguments");
-    return NULL;
-  }
-
-  // check that numpy array objects arrays are C contiguous memory
-  CHECK_C_CONTIG(neighbours);
-  CHECK_C_CONTIG(neighbour_edges);
-  CHECK_C_CONTIG(already_computed_flux);
-  CHECK_C_CONTIG(edge_timestep);
-  CHECK_C_CONTIG(flux_update_frequency);
-  CHECK_C_CONTIG(update_extrapolation);
-  CHECK_C_CONTIG(update_next_flux);
-  CHECK_C_CONTIG(allow_timestep_increase);
-
-  int number_of_elements = update_extrapolation -> dimensions[0];
-  // flux_update_frequency determined by rounding edge_timestep/timestep*notSoFast
-  // So notSoFast<1 might increase the number of flux calls a cell has to do, but
-  // can be useful for increasing numerical stability
-  double notSoFast=1.00;
-
-  // Call underlying flux computation routine and update 
-  // the explicit update arrays 
-  _compute_flux_update_frequency(number_of_elements,
-                                 (long*) neighbours->data,
-                                 (long*) neighbour_edges->data,
-                                 (long*) already_computed_flux->data,
-                                 max_flux_update_frequency,
-                                 (double*) edge_timestep->data,
-                                 timestep,
-                                 notSoFast,
-                                 (long*) flux_update_frequency->data,
-                                 (long*) update_extrapolation->data,
-                                 (long*) update_next_flux->data,
-                                 (long*) allow_timestep_increase->data
-                                 );
+  get_python_domain(&D,domain);
+
+  _compute_flux_update_frequency(&D, timestep);
+
   // Return 
   return Py_BuildValue("");
@@ -2412 +2112 @@
     on each triangle
     
-    Python call:
-    extrapolate_second_order_sw(domain.surrogate_neighbours,
-                                domain.number_of_boundaries
-                                domain.centroid_coordinates,
-                                Stage.centroid_values
-                                Xmom.centroid_values
-                                Ymom.centroid_values
-                                domain.edge_coordinates,
-                                Stage.edge_values,
-                                Xmom.edge_values,
-                                Ymom.edge_values)
-
     Post conditions:
         The edges of each triangle have values from a 
@@ -2430 +2118 @@
 
   */
-  PyArrayObject *surrogate_neighbours,
-    *neighbour_edges,
-    *number_of_boundaries,
-    *centroid_coordinates,
-    *stage_centroid_values,
-    *xmom_centroid_values,
-    *ymom_centroid_values,
-    *elevation_centroid_values,
-    *height_centroid_values,
-    *edge_coordinates,
-    *stage_edge_values,
-    *xmom_edge_values,
-    *ymom_edge_values,
-    *elevation_edge_values,
-    *height_edge_values,
-    *stage_vertex_values,
-    *xmom_vertex_values,
-    *ymom_vertex_values,
-    *elevation_vertex_values,
-    *height_vertex_values,
-    *x_centroid_work,
-    *y_centroid_work,
-    *update_extrapolation;
-  
+ 
+  struct domain D; 
   PyObject *domain;
 
-  
-  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;    
-  double minimum_allowed_height, epsilon;
-  int optimise_dry_cells, number_of_elements, extrapolate_velocity_second_order, e, e2;
-  
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOOOOOiiOOO",
-            &domain,
-            &surrogate_neighbours,
-            &neighbour_edges,
-            &number_of_boundaries,
-            &centroid_coordinates,
-            &stage_centroid_values,
-            &xmom_centroid_values,
-            &ymom_centroid_values,
-            &elevation_centroid_values,
-            &height_centroid_values,
-            &edge_coordinates,
-            &stage_edge_values,
-            &xmom_edge_values,
-            &ymom_edge_values,
-            &elevation_edge_values,
-            &height_edge_values,
-            &stage_vertex_values,
-            &xmom_vertex_values,
-            &ymom_vertex_values,
-            &elevation_vertex_values,
-            &height_vertex_values,
-            &optimise_dry_cells,
-            &extrapolate_velocity_second_order,
-            &x_centroid_work,
-            &y_centroid_work,
-            &update_extrapolation)) {         
-
-    report_python_error(AT, "could not parse input arguments");
-    return NULL;
-  }
-
-  // check that numpy array objects arrays are C contiguous memory
-  CHECK_C_CONTIG(surrogate_neighbours);
-  CHECK_C_CONTIG(neighbour_edges);
-  CHECK_C_CONTIG(number_of_boundaries);
-  CHECK_C_CONTIG(centroid_coordinates);
-  CHECK_C_CONTIG(stage_centroid_values);
-  CHECK_C_CONTIG(xmom_centroid_values);
-  CHECK_C_CONTIG(ymom_centroid_values);
-  CHECK_C_CONTIG(elevation_centroid_values);
-  CHECK_C_CONTIG(height_centroid_values);
-  CHECK_C_CONTIG(edge_coordinates);
-  CHECK_C_CONTIG(stage_edge_values);
-  CHECK_C_CONTIG(xmom_edge_values);
-  CHECK_C_CONTIG(ymom_edge_values);
-  CHECK_C_CONTIG(elevation_edge_values);
-  CHECK_C_CONTIG(height_edge_values);
-  CHECK_C_CONTIG(stage_vertex_values);
-  CHECK_C_CONTIG(xmom_vertex_values);
-  CHECK_C_CONTIG(ymom_vertex_values);
-  CHECK_C_CONTIG(elevation_vertex_values);
-  CHECK_C_CONTIG(height_vertex_values);
-  CHECK_C_CONTIG(x_centroid_work);
-  CHECK_C_CONTIG(y_centroid_work);
-  CHECK_C_CONTIG(update_extrapolation);
-  
-  // Get the safety factor beta_w, set in the config.py file. 
-  // This is used in the limiting process
-  
-
-  beta_w                 = get_python_double(domain,"beta_w");
-  beta_w_dry             = get_python_double(domain,"beta_w_dry");
-  beta_uh                = get_python_double(domain,"beta_uh");
-  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
-  beta_vh                = get_python_double(domain,"beta_vh");
-  beta_vh_dry            = get_python_double(domain,"beta_vh_dry");  
-
-  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
-  epsilon                = get_python_double(domain,"epsilon");
-
-  number_of_elements = stage_centroid_values -> dimensions[0];  
-
-  //printf("In C before Extrapolate");
-  //e=1;
-  // Call underlying computational routine
-  e = _extrapolate_second_order_edge_sw(number_of_elements,
-                   epsilon,
-                   minimum_allowed_height,
-                   beta_w,
-                   beta_w_dry,
-                   beta_uh,
-                   beta_uh_dry,
-                   beta_vh,
-                   beta_vh_dry,
-                   (long*) surrogate_neighbours -> data,
-                   (long*) neighbour_edges -> data,
-                   (long*) number_of_boundaries -> data,
-                   (double*) centroid_coordinates -> data,
-                   (double*) stage_centroid_values -> data,
-                   (double*) xmom_centroid_values -> data,
-                   (double*) ymom_centroid_values -> data,
-                   (double*) elevation_centroid_values -> data,
-                   (double*) height_centroid_values -> data,
-                   (double*) edge_coordinates -> data,
-                   (double*) stage_edge_values -> data,
-                   (double*) xmom_edge_values -> data,
-                   (double*) ymom_edge_values -> data,
-                   (double*) elevation_edge_values -> data,
-                   (double*) height_edge_values -> data,
-                   (double*) stage_vertex_values -> data,
-                   (double*) xmom_vertex_values -> data,
-                   (double*) ymom_vertex_values -> data,
-                   (double*) elevation_vertex_values -> data,
-                   (double*) height_vertex_values -> data,
-                   optimise_dry_cells, 
-                   extrapolate_velocity_second_order,
-                   (double*) x_centroid_work -> data,
-                   (double*) y_centroid_work -> data,
-                   (long*) update_extrapolation -> data);
+  int e;
+  
+  if (!PyArg_ParseTuple(args, "O", &domain)) {
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
+  }
+  
+  get_python_domain(&D, domain);
+
+  // Call underlying flux computation routine and update
+  // the explicit update arrays
+  e = _extrapolate_second_order_edge_sw(&D);
 
   if (e == -1) {
@@ -2582 +2144 @@
   
 }// extrapolate_second-order_edge_sw
+
 //========================================================================
 // Protect -- to prevent the water level from falling below the minimum 

trunk/anuga_core/source/anuga/shallow_water/sw_domain.h

@@ -21 +21 @@
     double  minimum_allowed_height;
 
+    long timestep_fluxcalls;
+
     double beta_w;
     double beta_w_dry;
@@ -28 +30 @@
     double beta_vh_dry;
 
+    long max_flux_update_frequency;
+    long ncol_riverwall_hydraulic_properties;
 
     // Changing values in these arrays will change the values in the python object
@@ -38 +42 @@
     double* areas;
 
-
+    long* edge_flux_type;
 
     long*   tri_full_flag;
@@ -53 +57 @@
     double* ymom_edge_values;
     double* bed_edge_values;
+    double* height_edge_values;
 
     double* stage_centroid_values;
@@ -58 +63 @@
     double* ymom_centroid_values;
     double* bed_centroid_values;
+    double* height_centroid_values;
 
     double* stage_vertex_values;
@@ -63 +69 @@
     double* ymom_vertex_values;
     double* bed_vertex_values;
+    double* height_vertex_values;
 
 
@@ -73 +80 @@
     double* xmom_explicit_update;
     double* ymom_explicit_update;
+
+    long* flux_update_frequency;    
+    long* update_next_flux;
+    long* update_extrapolation;
+    double* edge_timestep;
+    double* edge_flux_work;
+    double* pressuregrad_work;
+    double* x_centroid_work;
+    double* y_centroid_work;
+    double* boundary_flux_sum;
+
+    long* allow_timestep_increase;
+
+    double* riverwall_elevation;
+    long* riverwall_rowIndex;
+    double* riverwall_hydraulic_properties;
 };
 
@@ -152 +175 @@
             *radii,
             *areas,
+            *edge_flux_type,
             *tri_full_flag,
             *already_computed_flux,
@@ -159 +183 @@
             *number_of_boundaries,
             *surrogate_neighbours,
-            *max_speed;
+            *max_speed,
+            *flux_update_frequency,
+            *update_next_flux,
+            *update_extrapolation,
+            *allow_timestep_increase,
+            *edge_timestep,
+            *edge_flux_work,
+            *pressuregrad_work,
+            *x_centroid_work,
+            *y_centroid_work,
+            *boundary_flux_sum,
+            *riverwall_elevation,
+            *riverwall_rowIndex,
+            *riverwall_hydraulic_properties;
 
     PyObject *quantities;
+    PyObject *riverwallData;
 
     D->number_of_elements   = get_python_integer(domain, "number_of_elements");
@@ -170 +208 @@
     D->evolve_max_timestep  = get_python_double(domain, "evolve_max_timestep");
     D->minimum_allowed_height = get_python_double(domain, "minimum_allowed_height");
+    D->timestep_fluxcalls = get_python_integer(domain,"timestep_fluxcalls");
+    
 
     D->extrapolate_velocity_second_order  = get_python_integer(domain, "extrapolate_velocity_second_order");
@@ -180 +220 @@
     D->beta_vh_dry = get_python_double(domain, "beta_vh_dry");
 
-
+    D->max_flux_update_frequency = get_python_integer(domain,"max_flux_update_frequency");
     
     neighbours = get_consecutive_array(domain, "neighbours");
@@ -203 +243 @@
     D->areas = (double *) areas->data;
 
+    edge_flux_type = get_consecutive_array(domain, "edge_flux_type");
+    D->edge_flux_type = (long *) edge_flux_type->data;
+
+
     tri_full_flag = get_consecutive_array(domain, "tri_full_flag");
     D->tri_full_flag = (long *) tri_full_flag->data;
@@ -225 +269 @@
     D->number_of_boundaries = (long *) number_of_boundaries->data;
 
-
+    flux_update_frequency = get_consecutive_array(domain, "flux_update_frequency");
+    D->flux_update_frequency = (long*) flux_update_frequency->data;
+    
+    update_next_flux = get_consecutive_array(domain, "update_next_flux");
+    D->update_next_flux = (long*) update_next_flux->data;
+    
+    update_extrapolation = get_consecutive_array(domain, "update_extrapolation");
+    D->update_extrapolation = (long*) update_extrapolation->data;
+    
+    allow_timestep_increase = get_consecutive_array(domain, "allow_timestep_increase");
+    D->allow_timestep_increase = (long*) allow_timestep_increase->data;
+
+    edge_timestep = get_consecutive_array(domain, "edge_timestep");
+    D->edge_timestep = (double*) edge_timestep->data;
+    
+    edge_flux_work = get_consecutive_array(domain, "edge_flux_work");
+    D->edge_flux_work = (double*) edge_flux_work->data;
+    
+    pressuregrad_work = get_consecutive_array(domain, "pressuregrad_work");
+    D->pressuregrad_work = (double*) pressuregrad_work->data;
+    
+    x_centroid_work = get_consecutive_array(domain, "x_centroid_work");
+    D->x_centroid_work = (double*) x_centroid_work->data;
+
+    y_centroid_work = get_consecutive_array(domain, "y_centroid_work");
+    D->y_centroid_work = (double*) y_centroid_work->data;
+    
+    boundary_flux_sum = get_consecutive_array(domain, "boundary_flux_sum");
+    D->boundary_flux_sum = (double*) boundary_flux_sum->data;
 
     quantities = get_python_object(domain, "quantities");
@@ -233 +305 @@
     D->ymom_edge_values      = get_python_array_data_from_dict(quantities, "ymomentum", "edge_values");
     D->bed_edge_values       = get_python_array_data_from_dict(quantities, "elevation", "edge_values");
+    D->height_edge_values    = get_python_array_data_from_dict(quantities, "height", "edge_values");
 
     D->stage_centroid_values     = get_python_array_data_from_dict(quantities, "stage",     "centroid_values");
@@ -238 +311 @@
     D->ymom_centroid_values      = get_python_array_data_from_dict(quantities, "ymomentum", "centroid_values");
     D->bed_centroid_values       = get_python_array_data_from_dict(quantities, "elevation", "centroid_values");
+    D->height_centroid_values    = get_python_array_data_from_dict(quantities, "height", "centroid_values");
 
     D->stage_vertex_values     = get_python_array_data_from_dict(quantities, "stage",     "vertex_values");
@@ -243 +317 @@
     D->ymom_vertex_values      = get_python_array_data_from_dict(quantities, "ymomentum", "vertex_values");
     D->bed_vertex_values       = get_python_array_data_from_dict(quantities, "elevation", "vertex_values");
+    D->height_vertex_values       = get_python_array_data_from_dict(quantities, "height", "vertex_values");
 
     D->stage_boundary_values = get_python_array_data_from_dict(quantities, "stage",     "boundary_values");
@@ -254 +329 @@
 
 
+    riverwallData = get_python_object(domain,"riverwallData");
+
+    riverwall_elevation = get_consecutive_array(riverwallData, "riverwall_elevation");
+    D->riverwall_elevation = (double*) riverwall_elevation->data;
+    
+    riverwall_rowIndex = get_consecutive_array(riverwallData, "hydraulic_properties_rowIndex");
+    D->riverwall_rowIndex = (long*) riverwall_rowIndex->data;
+
+    D->ncol_riverwall_hydraulic_properties = get_python_integer(riverwallData, "ncol_hydraulic_properties");
+
+    riverwall_hydraulic_properties = get_consecutive_array(riverwallData, "hydraulic_properties");
+    D->riverwall_hydraulic_properties = (double*) riverwall_hydraulic_properties->data;
+
     Py_DECREF(quantities);
+    Py_DECREF(riverwallData);
 
     Py_DECREF(neighbours);
@@ -263 +352 @@
     Py_DECREF(radii);
     Py_DECREF(areas);
+    Py_DECREF(edge_flux_type);
     Py_DECREF(tri_full_flag);
     Py_DECREF(already_computed_flux);
@@ -270 +360 @@
     Py_DECREF(max_speed);
     Py_DECREF(number_of_boundaries);
+    Py_DECREF(flux_update_frequency);
+    Py_DECREF(update_next_flux);
+    Py_DECREF(update_extrapolation);
+    Py_DECREF(edge_timestep);
+    Py_DECREF(edge_flux_work);
+    Py_DECREF(pressuregrad_work);
+    Py_DECREF(x_centroid_work);
+    Py_DECREF(y_centroid_work);
+    Py_DECREF(boundary_flux_sum);
+    Py_DECREF(allow_timestep_increase);
 
     return D;
@@ -322 +422 @@
     printf("D->ymom_vertex_values     %p \n", D->ymom_vertex_values);
     printf("D->bed_vertex_values      %p \n", D->bed_vertex_values);
+    printf("D->height_vertex_values      %p \n", D->height_vertex_values);
     printf("D->stage_boundary_values  %p \n", D->stage_boundary_values);
     printf("D->xmom_boundary_values   %p \n", D->xmom_boundary_values);