source: trunk/anuga_work/development/gareth/balanced_basic/swb2_domain_ext.c @ 8884

// Python - C extension module for shallow_water.py
//
// To compile (Python2.6):
//  gcc -c swb2_domain_ext.c -I/usr/include/python2.6 -o domain_ext.o -Wall -O
//  gcc -shared swb2_domain_ext.o  -o swb2_domain_ext.so
//
// or use python compile.py
//
// See the module swb_domain.py for more documentation on 
// how to use this module
//
//
// Stephen Roberts, ANU 2009
// Ole Nielsen, GA 2004
// Gareth Davies, GA 2011

#include "Python.h"
#include "numpy/arrayobject.h"
#include "math.h"
#include <stdio.h>
//#include "numpy_shim.h"

// Shared code snippets
#include "util_ext.h"


const double pi = 3.14159265358979;

// Computational function for rotation
int _rotate(double *q, double n1, double n2) {
  /*Rotate the last  2 coordinates of q (q[1], q[2])
    from x,y coordinates to coordinates based on normal vector (n1, n2).

    Result is returned in array 2x1 r
    To rotate in opposite direction, call rotate with (q, n1, -n2)

    Contents of q are changed by this function */


  double q1, q2;

  // Shorthands
  q1 = q[1];  // x coordinate
  q2 = q[2];  // y coordinate

  // Rotate
  q[1] =  n1*q1 + n2*q2;
  q[2] = -n2*q1 + n1*q2;

  return 0;
}

// Function to obtain speed from momentum and depth.
// This is used by flux functions
// Input parameters uh and h may be modified by this function.
// Tried to inline, but no speedup was achieved 27th May 2009 (Ole)
//static inline double _compute_speed(double *uh, 
double _compute_speed(double *uh, 
                      double *h, 
                      double epsilon, 
                      double h0,
                      double limiting_threshold) {
  
  double u;

  if (*h < limiting_threshold) {   
    // Apply limiting of speeds according to the ANUGA manual
    if (*h < epsilon) {
      *h = 0.0;  // Could have been negative
      u = 0.0;
    } else {
      u = *uh/(*h + h0/ *h);    
    }
  

    // Adjust momentum to be consistent with speed
    *uh = u * *h;
  } else {
    // We are in deep water - no need for limiting
    u = *uh/ *h;
  }
  
  return u;
}

// Innermost flux function (using stage w=z+h)
int _flux_function_central(double *q_left, double *q_right,
                           double z_left, double z_right,
                           double n1, double n2,
                           double epsilon, 
                           double h0,
                           double limiting_threshold, 
                           double g,
                           double *edgeflux, double *max_speed, 
                           double *pressure_flux) 
{

  /*Compute fluxes between volumes for the shallow water wave equation
    cast in terms of the 'stage', w = h+z using
    the 'central scheme' as described in

    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.

    The implemented formula is given in equation (3.15) on page 714
  */

  int i;

  double w_left, h_left, uh_left, vh_left, u_left;
  double w_right, h_right, uh_right, vh_right, u_right;
  double s_min, s_max, soundspeed_left, soundspeed_right;
  double denom, inverse_denominator, z;
  double uint, t1, t2, t3;
  // Workspace (allocate once, use many)
  static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];

  // Copy conserved quantities to protect from modification
  q_left_rotated[0] = q_left[0];
  q_right_rotated[0] = q_right[0];
  q_left_rotated[1] = q_left[1];
  q_right_rotated[1] = q_right[1];
  q_left_rotated[2] = q_left[2];
  q_right_rotated[2] = q_right[2];

  // Align x- and y-momentum with x-axis
  _rotate(q_left_rotated, n1, n2);
  _rotate(q_right_rotated, n1, n2);

  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.

  // Compute speeds in x-direction
  w_left = q_left_rotated[0];          
  h_left = w_left - z;
  uh_left = q_left_rotated[1];
  u_left = _compute_speed(&uh_left, &h_left, 
                          epsilon, h0, limiting_threshold);

  w_right = q_right_rotated[0];
  h_right = w_right - z;
  uh_right = q_right_rotated[1];
  u_right = _compute_speed(&uh_right, &h_right, 
                           epsilon, h0, limiting_threshold);

  // Momentum in y-direction
  vh_left  = q_left_rotated[2];
  vh_right = q_right_rotated[2];

  // Limit y-momentum if necessary 
  // Leaving this out, improves speed significantly (Ole 27/5/2009)
  // All validation tests pass, so do we really need it anymore? 
  _compute_speed(&vh_left, &h_left, 
                 epsilon, h0, limiting_threshold);
  _compute_speed(&vh_right, &h_right, 
                 epsilon, h0, limiting_threshold);

  // Maximal and minimal wave speeds
  soundspeed_left  = sqrt(g*h_left);
  soundspeed_right = sqrt(g*h_right);  
  
  // Code to use fast square root optimisation if desired.
  // Timings on AMD 64 for the Okushiri profile gave the following timings
  //
  // SQRT           Total    Flux
  //=============================
  //
  // Ref            405s     152s
  // Fast (dbl)     453s     173s
  // Fast (sng)     437s     171s
  //
  // Consequently, there is currently (14/5/2009) no reason to use this 
  // approximation.
  
  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
  //soundspeed_right = fast_squareroot_approximation(g*h_right);

  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
  if (s_max < 0.0) 
  {
    s_max = 0.0;
  }

  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
  if (s_min > 0.0)
  {
    s_min = 0.0;
  }
  
  // Flux formulas
  flux_left[0] = u_left*h_left;
  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
  flux_left[2] = u_left*vh_left;

  flux_right[0] = u_right*h_right;
  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
  flux_right[2] = u_right*vh_right;

  // Flux computation
  denom = s_max - s_min;
  if (denom < epsilon) 
  { // FIXME (Ole): Try using h0 here
    memset(edgeflux, 0, 3*sizeof(double));
    //for(i=0;i<3;i++){
    //    edgeflux[i] = _minmod(flux_left[i], flux_right[i]);
    //}
    *max_speed = 0.0;
  } 
  else 
  {
    inverse_denominator = 1.0/denom;
    for (i = 0; i < 3; i++) 
    {
      // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
      // Physics 2:141-163
      //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
      //t1 = (q_right_rotated[i] - uint);
      //t2 = (-q_left_rotated[i] + uint);
      //t3 = _minmod(t1,t2);
      t3 = 0.0;
      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
      //t1 = s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
      //if(fabs(edgeflux[i])>fabs(t1)){ 
      //      edgeflux[i]+=t1;
      //}else{
      //      edgeflux[i]+=fabs(edgeflux[i])*copysign(1.0,t1) ;
      //}

      edgeflux[i] *= inverse_denominator;
    }

    *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
    //if(edgeflux[0]<0.0){
    //    //edgeflux[2] = edgeflux[0]*0.5*(vh_right/h_right + vh_left/h_left);
    //    edgeflux[2] = edgeflux[0]*vh_right/h_right;
    //}else{
    //    edgeflux[2] = edgeflux[0]*vh_left/h_left;
    //}

    // Maximal wavespeed
    *max_speed = max(fabs(s_max), fabs(s_min));

    // Rotate back
    _rotate(edgeflux, n1, -n2);
  }
  
  return 0;
}


// Computational function for flux computation
double _compute_fluxes_central(int number_of_elements,
        double timestep,
        double epsilon,
        double H0,
        double g,
        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* stage_boundary_values,
        double* xmom_boundary_values,
        double* ymom_boundary_values,
        double* stage_explicit_update,
        double* xmom_explicit_update,
        double* ymom_explicit_update,
        long* already_computed_flux,
        double* max_speed_array,
        int optimise_dry_cells, 
        double* stage_centroid_values,
        double* bed_centroid_values,
        double* bed_vertex_values) {
    // Local variables
    double max_speed, length, inv_area, zl, zr;
    double h0 = H0*H0; // This ensures a good balance when h approaches H0.

    double limiting_threshold = 10 * H0; // Avoid applying limiter below this
    // threshold for performance reasons.
    // See ANUGA manual under flux limiting
    int k, i, m, n,j;
    int ki, nm = 0, ki2; // Index shorthands

    // Workspace (making them static actually made function slightly slower (Ole))
    double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
    double stage_edges[3];//Work array
    double bedslope_work;
    int neighbours_wet[3];//Work array
    int useint;
    double usedouble, scale_factor_shallow, bedtop, bedbot, pressure_flux;
    //double depthtemp, max_bed_edge_value; 
    static long call = 1; // Static local variable flagging already computed flux

    // Start computation
    call++; // Flag 'id' of flux calculation for this timestep

    // 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));

    // For all triangles
    for (k = 0; k < number_of_elements; k++) {
        // Loop through neighbours and compute edge flux for each
        for (i = 0; i < 3; i++) {
            ki = k * 3 + i; // Linear index to edge i of triangle k

            if (already_computed_flux[ki] == call) {
                // We've already computed the flux across this edge
                continue;
            }

            // 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];

            // Get right hand side values either from neighbouring triangle
            // or from boundary array (Quantities at neighbour on nearest face).
            n = neighbours[ki];
            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];
                zr = zl; // Extend bed elevation to boundary
            }
            else {
                // Neighbour is a real triangle
                m = neighbour_edges[ki];
                nm = n * 3 + m; // Linear index (triangle n, edge m)

                qr[0] = stage_edge_values[nm];
                qr[1] = xmom_edge_values[nm];
                qr[2] = ymom_edge_values[nm];
                zr = bed_edge_values[nm];
            }

            // Now we have values for this edge - both from left and right side.

            if (optimise_dry_cells) {
                // Check if flux calculation is necessary across this edge
                // This check will exclude dry cells.
                // This will also optimise cases where zl != zr as
                // long as both are dry

                if ((ql[0] - zl) < epsilon &&
                        (qr[0] - zr) < epsilon) {
                    // Cell boundary is dry

                    already_computed_flux[ki] = call; // #k Done
                    if (n >= 0) {
                        already_computed_flux[nm] = call; // #n Done
                    }

                    max_speed = 0.0;
                    continue;
                }
            }


            if (fabs(zl-zr)>1.0e-10) {
                report_python_error(AT,"Discontinuous Elevation");
                return 0.0;
            }
            
            // If both centroids are dry, then the cells should not exchange flux 
            // NOTE: This also means no bedslope term -- which is arguably
            // appropriate, since the depth is zero at the cell centroid
            if( (stage_centroid_values[k]-bed_centroid_values[k]<=0.0)){
                if(n>=0){
                    if(stage_centroid_values[n]-bed_centroid_values[n]<=0.0){
                        continue;
                    }
                }else{
                    // Treat boundaries -- if the boundary stage < zl, there is
                    // no flux
                    if(qr[0]< zl){
                        continue;
                    }    
                }
            }  
           
            //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
            // unless the local centroid value is smaller 
            if(n>=0){
                if((stage_centroid_values[k]-bed_centroid_values[k]<= 0.0)){
                          ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
                }
                if(stage_centroid_values[n]-bed_centroid_values[n]<= 0.0){
                        qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
                }
            }else{
                // Treat the boundary case
                if((stage_centroid_values[k]-bed_centroid_values[k]<=0.0)){
                          ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
                }

            }
            
            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
            ki2 = 2 * ki; //k*6 + i*2

            // Edge flux computation (triangle k, edge i)
            _flux_function_central(ql, qr, zl, zr,
                    normals[ki2], normals[ki2 + 1],
                    epsilon, h0, limiting_threshold, g,
                    edgeflux, &max_speed, &pressure_flux);
            //_flux_function_fraccaro(ql, qr, zl, zr,
            //        normals[ki2], normals[ki2 + 1],
            //        epsilon, h0, limiting_threshold, g,
            //        edgeflux, &max_speed);
   
            //Find the bed vertex values associated with the edge 
            //if(i==0){
            //    bedtop = max(bed_vertex_values[3*k+1], bed_vertex_values[3*k+2]);
            //    bedbot = min(bed_vertex_values[3*k+1], bed_vertex_values[3*k+2]);
            //}else{
            //    if(i==1){
            //    bedtop = max(bed_vertex_values[3*k+0], bed_vertex_values[3*k+2]);
            //    bedbot = min(bed_vertex_values[3*k+0], bed_vertex_values[3*k+2]);
            //    }else{
            //    // i==2
            //    bedtop = max(bed_vertex_values[3*k+1], bed_vertex_values[3*k+0]);
            //    bedbot = min(bed_vertex_values[3*k+1], bed_vertex_values[3*k+0]);
            //    }
            //
            //}       

            //if(fabs(bedtop-bedbot)<1.0e-10){
            //    scale_factor_shallow = 1.0;
            //}else{
            //    scale_factor_shallow = min( (bedtop - bedbot)/(min(stage_centroid_values[k], bedtop) - bedbot),10.0);
            //}
            
            // Multiply edgeflux by edgelength
            length = edgelengths[ki];
            //length = edgelengths[ki]*scale_factor_shallow;
            edgeflux[0] *= length;
            edgeflux[1] *= length;
            edgeflux[2] *= length;
           
            // To correctly implement reflective boundary conditions would be something like this? 
            //if(n<0){
            //   _rotate(edgeflux,normals[ki2],normals[ki2+1]);
            //   edgeflux[1] = 0.0;
            //   _rotate(edgeflux,-normals[ki2+1],normals[ki2]);
            //}

            // Update triangle k with flux from edge i
            stage_explicit_update[k] -= edgeflux[0];
            xmom_explicit_update[k] -= edgeflux[1];
            ymom_explicit_update[k] -= edgeflux[2];

            // Calculate the bed slope term, using the 'divergence form' from 
            // Alessandro Valiani and Lorenzo Begnudelli, Divergence Form for Bed Slope Source Term in Shallow Water Equations
            // Journal of Hydraulic Engineering, 2006, 132:652-665
            if((stage_centroid_values[k]>bed_centroid_values[k]+0.0)){
                bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
            }else{
                // Treat nearly dry cells 
                //bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
                bedslope_work = -pressure_flux*length;
                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
            }
                // The triangle is (partially?) dry
                // Find the average stage of the wet neighbours edges -- set this as the within triangle stage for the bedslope term
                // FIXME: Note that we keep recomputing this -- could improve efficiency 
                //for(j=0;j<3;j++){ 
                //    // For each edge, check if each neighbouring triangle is wet
                //    useint = neighbours[3*k+j];
                //    if(useint>0){
                //    // neighbour is not a boundary condition
                //        if(stage_centroid_values[useint] > bed_centroid_values[useint]){
                //            //wet
                //            neighbours_wet[j]=1;
                //            useint = useint*3 + neighbour_edges[3*k+j]; // Index of neighbouring edge
                //            stage_edges[j] = stage_edge_values[useint]; // stage of neighbouring edge
                //        }else{
                //            //dry
                //            neighbours_wet[j]=0;
                //            stage_edges[j]=0.0; // Arbitrary
                //        }
                //    }else{ 
                //        // neighbour is a boundary condition, ASSUME IT IS WET
                //        useint = -useint -1;
                //        neighbours_wet[j]=1;
                //        stage_edges[j] = stage_boundary_values[useint];
                //        }         
                //}
                ////calculate average neighbouring edge stage, only considering wet triangles
                //if(max(neighbours_wet[0], max(neighbours_wet[1],neighbours_wet[2]))>0){
                //    //usedouble == average edge stage
                //    usedouble = stage_edges[0]*neighbours_wet[0] + stage_edges[1]*neighbours_wet[1] + stage_edges[2]*neighbours_wet[2];
                //    usedouble /= 1.0*(neighbours_wet[0]+neighbours_wet[1]+neighbours_wet[2]);
                //}else{
                //    usedouble = zl;
                //}
                //usedouble=min(usedouble,stage_centroid_values[k]);
                //xmom_explicit_update[k] -= -0.5*g*max(usedouble-zl,0.0)*(usedouble-zl)*normals[ki2]*length;
                //ymom_explicit_update[k] -= -0.5*g*max(usedouble-zl,0.0)*(usedouble-zl)*normals[ki2+1]*length;
                ////xmom_explicit_update[k] -= -0.5*g*abs(usedouble-zl)*(usedouble-zl)*normals[ki2]*length;
                ////ymom_explicit_update[k] -= -0.5*g*abs(usedouble-zl)*(usedouble-zl)*normals[ki2+1]*length;

            //}

            already_computed_flux[ki] = call; // #k Done


            // Update neighbour n with same flux but reversed sign
            if (n >= 0) {
                stage_explicit_update[n] += edgeflux[0];
                xmom_explicit_update[n] += edgeflux[1];
                ymom_explicit_update[n] += edgeflux[2];
                //Add bed slope term here
                if((stage_centroid_values[n]>bed_centroid_values[n]+0.0)){
                    bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
                }else{
                    // Treat nearly dry cells
                    //bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
                    bedslope_work = -pressure_flux*length;
                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
                }
                    // The triangle is (partially?) dry
                    // Find the average stage of the wet neighbours edges -- set this as the within triangle stage for the bedslope term
                    // FIXME: Note that we keep recomputing this -- could improve efficiency
                    //for(j=0;j<3;j++){
                    //    // For each edge, check if each neighbouring triangle is wet
                    //    useint = neighbours[3*n+j];
                    //    if(useint>0){
                    //    // neighbour is not a boundary condition
                    //        if(stage_centroid_values[useint] > bed_centroid_values[useint]){
                    //            //wet
                    //            neighbours_wet[j]=1;
                    //            useint = useint*3 + neighbour_edges[3*n+j]; // Index of neighbouring edge
                    //            stage_edges[j] = stage_edge_values[useint]; // stage of neighbouring edge
                    //        }else{
                    //            //dry
                    //            neighbours_wet[j]=0;
                    //            stage_edges[j]=0.0; //Arbitrary
                    //        }
                    //    }else{ 
                    //        // neighbour is a boundary condition, ASSUME IT IS WET
                    //        useint = -useint -1;
                    //        neighbours_wet[j]=1;
                    //        stage_edges[j] = stage_boundary_values[useint];
                    //        }         
                    //}
                    ////calculate average neighbouring edge stage, only considering wet triangles
                    //if(max(neighbours_wet[0],max(neighbours_wet[1],neighbours_wet[2]))>0){
                    //    //average edge stage
                    //    usedouble = stage_edges[0]*neighbours_wet[0] + stage_edges[1]*neighbours_wet[1] + stage_edges[2]*neighbours_wet[2];
                    //    usedouble /= 1.0*(neighbours_wet[0]+neighbours_wet[1]+neighbours_wet[2]);
                    //}else{
                    //    usedouble = zr;
                    //}
                    //usedouble=min(usedouble,stage_centroid_values[n]);
                    //xmom_explicit_update[n] += -0.5*g*max(usedouble-zr,0.0)*(usedouble-zr)*normals[ki2]*length;
                    //ymom_explicit_update[n] += -0.5*g*max(usedouble-zr,0.0)*(usedouble-zr)*normals[ki2+1]*length;
                    ////xmom_explicit_update[n] += -0.5*g*abs(usedouble-zr)*(usedouble-zr)*normals[ki2]*length;
                    ////ymom_explicit_update[n] += -0.5*g*abs(usedouble-zr)*(usedouble-zr)*normals[ki2+1]*length;

               //}

                already_computed_flux[nm] = call; // #n Done
            }

            // Update timestep based on edge i and possibly neighbour n
            if (tri_full_flag[k] == 1) {
                if (max_speed > epsilon) {
                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)

                    // CFL for triangle k
                    timestep = min(timestep, radii[k] / max_speed);

                    if (n >= 0) {
                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
                        timestep = min(timestep, radii[n] / max_speed);
                    }

                    // Ted Rigby's suggested less conservative version
                    //if(n>=0){
                    //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
                    //}else{
                    //  timestep = min(timestep, radii[k]/max_speed);
                    //}
                }
            }

        } // End edge i (and neighbour n)


        // 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;


        // Keep track of maximal speeds
        max_speed_array[k] = max_speed;

    } // End triangle k


    // Try to adjust stage_explicit_update to account for the negative stages,
    // using an approach assumes a very small discharge is enough to allow 
    // the depth to increase from negative to zero
    //for(k=0;k<number_of_elements;k++){
        //depthtemp = stage_centroid_values[k]-bed_centroid_values[k];
    //    depthtemp = stage_centroid_values[k]-max(stage_edge_values[3*k+1], max(stage_edge_values[3*k+2], stage_edge_values[3*k]));
    //   if(depthtemp< 0.0){
    //       depthtemp = -depthtemp; // stage defecit
    //       usedouble = 2.0; //1.0e+2; //Large constant
    //       if(stage_explicit_update[k]*timestep*usedouble > depthtemp){
    //            stage_explicit_update[k] += (depthtemp - depthtemp/usedouble)/timestep;
    //        }else{
    //            stage_explicit_update[k] *=usedouble ;
    //        } 
    //    }
    //}

    return timestep;
}

// Protect against the water elevation falling below the triangle bed
int _protect(int N,
         double minimum_allowed_height,
         double maximum_allowed_speed,
         double epsilon,
         double* wc,
         double* zc,
         double* zv,
         double* xmomc,
         double* ymomc) {

  int k;
  double hc, wmin;
  double u, v, reduced_speed;


  // Protect against initesimal and negative heights  
  if (maximum_allowed_speed < epsilon) {
    for (k=0; k<N; k++) {
      hc = wc[k] - zc[k];

      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){
             // Definine the minimum possible stage value as the minimum bed edge value on triangle k.
             wmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
             wc[k] = max(wc[k], wmin);
        }
      }
    }
  } else {
   
    // Protect against initesimal and negative heights
    for (k=0; k<N; k++) {
      hc = wc[k] - zc[k];
      
      if (hc < minimum_allowed_height) {

        //New code: Adjust momentum to guarantee speeds are physical
        //          ensure h is non negative
              
        if (hc <= 0.0) {
            wc[k] = zc[k];
        xmomc[k] = 0.0;
        ymomc[k] = 0.0;
        } else {
          //Reduce excessive speeds derived from division by small hc
          //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;
          }
          
          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;
          }
        }
      }
    }
  }
  return 0;
}

//=========================================================================
// Python Glue
//=========================================================================


//========================================================================
// Compute fluxes
//========================================================================

// Modified central flux function

PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
  /*Compute all fluxes and the timestep suitable for all volumes
    in domain.

    Compute total flux for each conserved quantity using "flux_function_central"

    Fluxes across each edge are scaled by edgelengths and summed up
    Resulting flux is then scaled by area and stored in
    explicit_update for each of the three conserved quantities
    stage, xmomentum and ymomentum

    The maximal allowable speed computed by the flux_function for each volume
    is converted to a timestep that must not be exceeded. The minimum of
    those is computed as the next overall timestep.

    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)


    Post conditions:
      domain.explicit_update is reset to computed flux values
      domain.timestep is set to the largest step satisfying all volumes.


  */


  PyArrayObject *neighbours, *neighbour_edges,
    *normals, *edgelengths, *radii, *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, //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,
    *bed_centroid_values,
    *bed_vertex_values;
    
  double timestep, epsilon, H0, g;
  int optimise_dry_cells;
    
  // Convert Python arguments to C
  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOiOOO",
            &timestep,
            &epsilon,
            &H0,
            &g,
            &neighbours,
            &neighbour_edges,
            &normals,
            &edgelengths, &radii, &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,
            &max_speed_array,
            &optimise_dry_cells,
            &stage_centroid_values,
            &bed_centroid_values,
            &bed_vertex_values)) {
    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(stage_boundary_values);
  CHECK_C_CONTIG(xmom_boundary_values);
  CHECK_C_CONTIG(ymom_boundary_values);
  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(bed_centroid_values);
  CHECK_C_CONTIG(bed_vertex_values);

  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,
                     (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*) stage_boundary_values -> data,
                     (double*) xmom_boundary_values -> data,
                     (double*) ymom_boundary_values -> 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, 
                     (double*) stage_centroid_values -> data, 
                     (double*) bed_centroid_values -> data,
                     (double*) bed_vertex_values -> data);
  
  // Return updated flux timestep
  return Py_BuildValue("d", timestep);
}


PyObject *flux_function_central(PyObject *self, PyObject *args) {
  //
  // Gateway to innermost flux function.
  // This is only used by the unit tests as the c implementation is
  // normally called by compute_fluxes in this module.


  PyArrayObject *normal, *ql, *qr,  *edgeflux;
  double g, epsilon, max_speed, H0, zl, zr;
  double h0, limiting_threshold, pressure_flux;

  if (!PyArg_ParseTuple(args, "OOOddOddd",
            &normal, &ql, &qr, &zl, &zr, &edgeflux,
            &epsilon, &g, &H0)) {
      report_python_error(AT, "could not parse input arguments");
      return NULL;
  }

  
  h0 = H0*H0; // This ensures a good balance when h approaches H0.
              // But evidence suggests that h0 can be as little as
              // epsilon!
              
  limiting_threshold = 10*H0; // Avoid applying limiter below this
                              // threshold for performance reasons.
                              // See ANUGA manual under flux limiting  
 
  pressure_flux = 0.0; // Store the water pressure related component of the flux 
  _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,
             &pressure_flux);
  
  return Py_BuildValue("d", max_speed);  
}

//========================================================================
// Gravity
//========================================================================

PyObject *gravity(PyObject *self, PyObject *args) {
  //
  //  gravity(g, h, v, x, xmom, ymom)
  //
  
  
  PyArrayObject *h, *v, *x, *xmom, *ymom;
  int k, N, k3, k6;
  double g, avg_h, zx, zy;
  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
  //double epsilon;
  
  if (!PyArg_ParseTuple(args, "dOOOOO",
                        &g, &h, &v, &x,
                        &xmom, &ymom)) {
    //&epsilon)) {
    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
    return NULL;
  }
  
  // check that numpy array objects arrays are C contiguous memory
  CHECK_C_CONTIG(h);
  CHECK_C_CONTIG(v);
  CHECK_C_CONTIG(x);
  CHECK_C_CONTIG(xmom);
  CHECK_C_CONTIG(ymom);
  
  N = h -> dimensions[0];
  for (k=0; k<N; k++) {
    k3 = 3*k;  // base index
    
    // Get bathymetry
    z0 = ((double*) v -> data)[k3 + 0];
    z1 = ((double*) v -> data)[k3 + 1];
    z2 = ((double*) v -> data)[k3 + 2];
    
    // Optimise for flat bed
    // Note (Ole): This didn't produce measurable speed up.
    // Revisit later
    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
    //  continue;
    //} 
    
    // Get average depth from centroid values
    avg_h = ((double *) h -> data)[k];
    
    // Compute bed slope
    k6 = 6*k;  // base index
    
    x0 = ((double*) x -> data)[k6 + 0];
    y0 = ((double*) x -> data)[k6 + 1];
    x1 = ((double*) x -> data)[k6 + 2];
    y1 = ((double*) x -> data)[k6 + 3];
    x2 = ((double*) x -> data)[k6 + 4];
    y2 = ((double*) x -> data)[k6 + 5];
    
    
    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
    
    // Update momentum
    ((double*) xmom -> data)[k] += -g*zx*avg_h;
    ((double*) ymom -> data)[k] += -g*zy*avg_h;
  }
  
  return Py_BuildValue("");
}


//========================================================================
// Protect -- to prevent the water level from falling below the minimum 
// bed_edge_value
//========================================================================

PyObject *protect(PyObject *self, PyObject *args) {
  //
  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)


  PyArrayObject
  *wc,            //Stage at centroids
  *zc,            //Elevation at centroids
  *zv,            //Elevation at vertices
  *xmomc,         //Momentums at centroids
  *ymomc;


  int N;
  double minimum_allowed_height, maximum_allowed_speed, epsilon;

  // Convert Python arguments to C
  if (!PyArg_ParseTuple(args, "dddOOOOO",
            &minimum_allowed_height,
            &maximum_allowed_speed,
            &epsilon,
            &wc, &zc, &zv, &xmomc, &ymomc)) {
    report_python_error(AT, "could not parse input arguments");
    return NULL;
  }  

  N = wc -> dimensions[0];

  _protect(N,
       minimum_allowed_height,
       maximum_allowed_speed,
       epsilon,
       (double*) wc -> data,
       (double*) zc -> data,
       (double*) zv -> data,
       (double*) xmomc -> data,
       (double*) ymomc -> data);

  return Py_BuildValue("");
}

//========================================================================
// Method table for python module
//========================================================================

static struct PyMethodDef MethodTable[] = {
  /* The cast of the function is necessary since PyCFunction values
   * only take two PyObject* parameters, and rotate() takes
   * three.
   */
  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
  {"gravity_c",        gravity,            METH_VARARGS, "Print out"},
  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},  
  {"protect",          protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
  {NULL, NULL, 0, NULL}
};

// Module initialisation
void initswb2_domain_ext(void){
  Py_InitModule("swb2_domain_ext", MethodTable);

  import_array(); // Necessary for handling of NumPY structures
}