source: anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c @ 4712

// Python - C extension for quantity module.
//
// To compile (Python2.3):
//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O
//  gcc -shared util_ext.o  -o util_ext.so
//
// See the module quantity.py
//
//
// Ole Nielsen, GA 2004

#include "Python.h"
#include "Numeric/arrayobject.h"
#include "math.h"

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



int _compute_gradients(int N,
                        double* centroids,
                        double* centroid_values,
                        long* number_of_boundaries,
                        long* surrogate_neighbours,
                        double* a,
                        double* b){

        int i, k, k0, k1, k2, index3;
        double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;


        for (k=0; k<N; k++) {
                index3 = 3*k;

                if (number_of_boundaries[k] < 2) {
                        //Two or three true neighbours

                        //Get indices of neighbours (or self when used as surrogate)
                        //k0, k1, k2 = surrogate_neighbours[k,:]

                        k0 = surrogate_neighbours[index3 + 0];
                        k1 = surrogate_neighbours[index3 + 1];
                        k2 = surrogate_neighbours[index3 + 2];


                        if (k0 == k1 || k1 == k2) return -1;

                        //Get data
                        q0 = centroid_values[k0];
                        q1 = centroid_values[k1];
                        q2 = centroid_values[k2];

                        x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
                        x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
                        x2 = centroids[k2*2]; y2 = centroids[k2*2+1];

                        //Gradient
                        _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);

                } else if (number_of_boundaries[k] == 2) {
                        //One true neighbour

                        //Get index of the one neighbour
                        i=0; k0 = k;
                        while (i<3 && k0==k) {
                                k0 = surrogate_neighbours[index3 + i];
                                i++;
                        }
                        if (k0 == k) return -1;

                        k1 = k; //self

                        //Get data
                        q0 = centroid_values[k0];
                        q1 = centroid_values[k1];

                        x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
                        x1 = centroids[k1*2]; y1 = centroids[k1*2+1];

                        //Two point gradient
                        _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);


                        //Old (wrong code)
                        //det = x0*y1 - x1*y0;
                        //if (det != 0.0) {
                        //      a[k] = (y1*q0 - y0*q1)/det;
                        //      b[k] = (x0*q1 - x1*q0)/det;
                        //}
                }
                //    else:
                //        #No true neighbours -
                //        #Fall back to first order scheme
                //        pass
        }
        return 0;
}



int _extrapolate(int N,
                 double* centroids,
                 double* centroid_values,
                 double* vertex_coordinates,
                 double* vertex_values,
                 double* a,
                 double* b) {

        int k, k2, k3, k6;
        double x, y, x0, y0, x1, y1, x2, y2;

        for (k=0; k<N; k++){
                k6 = 6*k;
                k3 = 3*k;
                k2 = 2*k;

                //Centroid coordinates
                x = centroids[k2]; y = centroids[k2+1];

                //vertex coordinates
                //x0, y0, x1, y1, x2, y2 = X[k,:]
                x0 = vertex_coordinates[k6 + 0];
                y0 = vertex_coordinates[k6 + 1];
                x1 = vertex_coordinates[k6 + 2];
                y1 = vertex_coordinates[k6 + 3];
                x2 = vertex_coordinates[k6 + 4];
                y2 = vertex_coordinates[k6 + 5];

                //Extrapolate
                vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
                vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
                vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);

        }
        return 0;
}




int _interpolate(int N,
                 double* vertex_values,
                 double* edge_values) {

        int k, k3;
        double q0, q1, q2;


        for (k=0; k<N; k++) {
                k3 = 3*k;

                q0 = vertex_values[k3 + 0];
                q1 = vertex_values[k3 + 1];
                q2 = vertex_values[k3 + 2];

                //printf("%f, %f, %f\n", q0, q1, q2);
                edge_values[k3 + 0] = 0.5*(q1+q2);
                edge_values[k3 + 1] = 0.5*(q0+q2);
                edge_values[k3 + 2] = 0.5*(q0+q1);
        }
        return 0;
}

int _backup_centroid_values(int N,
                            double* centroid_values,
                            double* centroid_backup_values) {
    //Backup centroid values


    int k;

    for (k=0; k<N; k++) {
        centroid_backup_values[k] = centroid_values[k];
    }


    return 0;
}


int _saxpy_centroid_values(int N,
                           double a,
                           double b,
                           double* centroid_values,
                           double* centroid_backup_values) {
    //saxby centroid values


    int k;


    for (k=0; k<N; k++) {
        centroid_values[k] = a*centroid_values[k] + b*centroid_backup_values[k];
    }


    return 0;
}


int _update(int N,
            double timestep,
            double* centroid_values,
            double* explicit_update,
            double* semi_implicit_update) {
        //Update centroid values based on values stored in
        //explicit_update and semi_implicit_update as well as given timestep


        int k;
        double denominator, x;


        //Divide semi_implicit update by conserved quantity
        for (k=0; k<N; k++) {
                x = centroid_values[k];
                if (x == 0.0) {
                        semi_implicit_update[k] = 0.0;
                } else {
                        semi_implicit_update[k] /= x;
                }
        }


        //Semi implicit updates
        for (k=0; k<N; k++) {
                denominator = 1.0 - timestep*semi_implicit_update[k];
                if (denominator == 0.0) {
                        return -1;
                } else {
                        //Update conserved_quantities from semi implicit updates
                        centroid_values[k] /= denominator;
                }
        }

        /*  for (k=0; k<N; k++) {*/
        /*    centroid_values[k] = exp(timestep*semi_implicit_update[k])*centroid_values[k];*/
        /*  }*/


        //Explicit updates
        for (k=0; k<N; k++) {
                centroid_values[k] += timestep*explicit_update[k];
        }


        //MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
        for (k=0;k<N;k++){
                semi_implicit_update[k]=0.0;
        }

        return 0;
}


int _average_vertex_values(int N,
                           long* vertex_value_indices,
                           long* number_of_triangles_per_node,
                           double* vertex_values, 
                           double* A) {
  //Average vertex values to obtain one value per node

  int i, index; 
  int k = 0; //Track triangles touching each node
  int current_node = 0;
  double total = 0.0;

  for (i=0; i<N; i++) {
  
    //if (current_node == N) {
    //  printf("Current node exceeding number of nodes (%d)", N);    
    //  return 1;
    // }
    
    index = vertex_value_indices[i];
    k += 1;
            
    //volume_id = index / 3
    //vertex_id = index % 3
    //total += self.vertex_values[volume_id, vertex_id]
    total += vertex_values[index];
      
    //printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);
    if (number_of_triangles_per_node[current_node] == k) {
      A[current_node] = total/k;
                
      // Move on to next node
      total = 0.0;
      k = 0;
      current_node += 1;
    }
  }
                           
  return 0;
}


/////////////////////////////////////////////////
// Gateways to Python
PyObject *update(PyObject *self, PyObject *args) {

        PyObject *quantity;
        PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;

        double timestep;
        int N, err;


        // Convert Python arguments to C
        if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep)) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "quantity_ext.c: update could not parse input");
          return NULL;
        }

        centroid_values = get_consecutive_array(quantity, "centroid_values");
        explicit_update = get_consecutive_array(quantity, "explicit_update");
        semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");

        N = centroid_values -> dimensions[0];

        err = _update(N, timestep,
                      (double*) centroid_values -> data,
                      (double*) explicit_update -> data,
                      (double*) semi_implicit_update -> data);


        if (err != 0) {
          PyErr_SetString(PyExc_RuntimeError,
                          "Zero division in semi implicit update - call Stephen :)");
          return NULL;
        }

        //Release and return
        Py_DECREF(centroid_values);
        Py_DECREF(explicit_update);
        Py_DECREF(semi_implicit_update);

        return Py_BuildValue("");
}


PyObject *backup_centroid_values(PyObject *self, PyObject *args) {

        PyObject *quantity;
        PyArrayObject *centroid_values, *centroid_backup_values;

        int N, err;


        // Convert Python arguments to C
        if (!PyArg_ParseTuple(args, "O", &quantity)) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "quantity_ext.c: backup_centroid_values could not parse input");
          return NULL;
        }

        centroid_values        = get_consecutive_array(quantity, "centroid_values");
        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");

        N = centroid_values -> dimensions[0];

        err = _backup_centroid_values(N,
                      (double*) centroid_values -> data,
                      (double*) centroid_backup_values -> data);


        //Release and return
        Py_DECREF(centroid_values);
        Py_DECREF(centroid_backup_values);

        return Py_BuildValue("");
}

PyObject *saxpy_centroid_values(PyObject *self, PyObject *args) {

        PyObject *quantity;
        PyArrayObject *centroid_values, *centroid_backup_values;

        double a,b;
        int N, err;


        // Convert Python arguments to C
        if (!PyArg_ParseTuple(args, "Odd", &quantity, &a, &b)) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "quantity_ext.c: saxpy_centroid_values could not parse input");
          return NULL;
        }

        centroid_values        = get_consecutive_array(quantity, "centroid_values");
        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");

        N = centroid_values -> dimensions[0];

        err = _saxpy_centroid_values(N,a,b,
                      (double*) centroid_values -> data,
                      (double*) centroid_backup_values -> data);


        //Release and return
        Py_DECREF(centroid_values);
        Py_DECREF(centroid_backup_values);

        return Py_BuildValue("");
}


PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {

        PyObject *quantity;
        PyArrayObject *vertex_values, *edge_values;

        int N, err;

        // Convert Python arguments to C
        if (!PyArg_ParseTuple(args, "O", &quantity)) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "quantity_ext.c: interpolate_from_vertices_to_edges could not parse input");
          return NULL;
        }
        
        vertex_values = get_consecutive_array(quantity, "vertex_values");
        edge_values = get_consecutive_array(quantity, "edge_values");

        N = vertex_values -> dimensions[0];

        err = _interpolate(N,
                           (double*) vertex_values -> data,
                           (double*) edge_values -> data);

        if (err != 0) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "Interpolate could not be computed");
          return NULL;
        }

        //Release and return
        Py_DECREF(vertex_values);
        Py_DECREF(edge_values);

        return Py_BuildValue("");
}


PyObject *average_vertex_values(PyObject *self, PyObject *args) {

        PyArrayObject 
          *vertex_value_indices, 
          *number_of_triangles_per_node,
          *vertex_values,
          *A;
        

        int N, err;

        // Convert Python arguments to C
        if (!PyArg_ParseTuple(args, "OOOO",
                              &vertex_value_indices, 
                              &number_of_triangles_per_node,
                              &vertex_values, 
                              &A)) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "quantity_ext.c: average_vertex_values could not parse input");
          return NULL;
        }
        
        N = vertex_value_indices -> dimensions[0];
        // printf("Got parameters, N=%d\n", N);
        err = _average_vertex_values(N,
                                     (long*) vertex_value_indices -> data,
                                     (long*) number_of_triangles_per_node -> data,
                                     (double*) vertex_values -> data, 
                                     (double*) A -> data);

        //printf("Error %d", err);
        if (err != 0) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "average_vertex_values could not be computed");
          return NULL;
        }

        return Py_BuildValue("");
}



PyObject *compute_gradients(PyObject *self, PyObject *args) {

        PyObject *quantity, *domain, *R;
        PyArrayObject
                *centroids,            //Coordinates at centroids
                *centroid_values,      //Values at centroids
                *number_of_boundaries, //Number of boundaries for each triangle
                *surrogate_neighbours, //True neighbours or - if one missing - self
                *a, *b;                //Return values

        int dimensions[1], N, err;

        // Convert Python arguments to C
        if (!PyArg_ParseTuple(args, "O", &quantity)) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "quantity_ext.c: compute_gradients could not parse input");   
          return NULL;
        }

        domain = PyObject_GetAttrString(quantity, "domain");
        if (!domain) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "compute_gradients could not obtain domain object from quantity");
          return NULL;
        }

        //Get pertinent variables

        centroids = get_consecutive_array(domain, "centroid_coordinates");
        centroid_values = get_consecutive_array(quantity, "centroid_values");
        surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
        number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");

        N = centroid_values -> dimensions[0];

        //Release
        Py_DECREF(domain);

        //Allocate space for return vectors a and b (don't DECREF)
        dimensions[0] = N;
        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);



        err = _compute_gradients(N,
                        (double*) centroids -> data,
                        (double*) centroid_values -> data,
                        (long*) number_of_boundaries -> data,
                        (long*) surrogate_neighbours -> data,
                        (double*) a -> data,
                        (double*) b -> data);

        if (err != 0) {
          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
          return NULL;
        }

        //Release
        Py_DECREF(centroids);
        Py_DECREF(centroid_values);
        Py_DECREF(number_of_boundaries);
        Py_DECREF(surrogate_neighbours);

        //Build result, release and return
        R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
        Py_DECREF(a);
        Py_DECREF(b);
        return R;
}



PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {

        PyObject *quantity, *domain;
        PyArrayObject
                *centroids,            //Coordinates at centroids
                *centroid_values,      //Values at centroids
                *vertex_coordinates,   //Coordinates at vertices
                *vertex_values,        //Values at vertices
                *number_of_boundaries, //Number of boundaries for each triangle
                *surrogate_neighbours, //True neighbours or - if one missing - self
                *a, *b;                //Gradients

        //int N, err;
        int dimensions[1], N, err;
        //double *a, *b;  //Gradients

        // Convert Python arguments to C
        if (!PyArg_ParseTuple(args, "O", &quantity)) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "extrapolate_second_order could not parse input");    
          return NULL;
        }

        domain = PyObject_GetAttrString(quantity, "domain");
        if (!domain) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "extrapolate_second_order could not obtain domain object from quantity");     
          return NULL;
        }

        //Get pertinent variables
        centroids = get_consecutive_array(domain, "centroid_coordinates");
        centroid_values = get_consecutive_array(quantity, "centroid_values");
        surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
        number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
        vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
        vertex_values = get_consecutive_array(quantity, "vertex_values");

        N = centroid_values -> dimensions[0];

        //Release
        Py_DECREF(domain);

        //Allocate space for return vectors a and b (don't DECREF)
        dimensions[0] = N;
        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);

        //FIXME: Odd that I couldn't use normal arrays
        //Allocate space for return vectors a and b (don't DECREF)
        //a = (double*) malloc(N * sizeof(double));
        //if (!a) return NULL;
        //b = (double*) malloc(N * sizeof(double));
        //if (!b) return NULL;


        err = _compute_gradients(N,
                        (double*) centroids -> data,
                        (double*) centroid_values -> data,
                        (long*) number_of_boundaries -> data,
                        (long*) surrogate_neighbours -> data,
                        (double*) a -> data,
                        (double*) b -> data);

        if (err != 0) {
          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
          return NULL;
        }

        err = _extrapolate(N,
                        (double*) centroids -> data,
                        (double*) centroid_values -> data,
                        (double*) vertex_coordinates -> data,
                        (double*) vertex_values -> data,
                        (double*) a -> data,
                        (double*) b -> data);


        if (err != 0) {
          PyErr_SetString(PyExc_RuntimeError,
                          "Internal function _extrapolate failed");
          return NULL;
        }



        //Release
        Py_DECREF(centroids);
        Py_DECREF(centroid_values);
        Py_DECREF(number_of_boundaries);
        Py_DECREF(surrogate_neighbours);
        Py_DECREF(vertex_coordinates);
        Py_DECREF(vertex_values);
        Py_DECREF(a);
        Py_DECREF(b);

        return Py_BuildValue("");
}



PyObject *limit(PyObject *self, PyObject *args) {

        PyObject *quantity, *domain, *Tmp;
        PyArrayObject
                *qv, //Conserved quantities at vertices
                *qc, //Conserved quantities at centroids
                *neighbours;

        int k, i, n, N, k3;
        double beta_w; //Safety factor
        double *qmin, *qmax, qn;

        // Convert Python arguments to C
        if (!PyArg_ParseTuple(args, "O", &quantity)) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "quantity_ext.c: limit could not parse input");
          return NULL;
        }

        domain = PyObject_GetAttrString(quantity, "domain");
        if (!domain) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "quantity_ext.c: limit could not obtain domain object from quantity");                        
          
          return NULL;
        }

        //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
        neighbours = get_consecutive_array(domain, "neighbours");

        //Get safety factor beta_w
        Tmp = PyObject_GetAttrString(domain, "beta_w");
        if (!Tmp) {
          PyErr_SetString(PyExc_RuntimeError, 
                          "quantity_ext.c: limit could not obtain beta_w object from domain");                  
          
          return NULL;
        }       

        beta_w = PyFloat_AsDouble(Tmp);

        Py_DECREF(Tmp);
        Py_DECREF(domain);


        qc = get_consecutive_array(quantity, "centroid_values");
        qv = get_consecutive_array(quantity, "vertex_values");


        N = qc -> dimensions[0];

        //Find min and max of this and neighbour's centroid values
        qmin = malloc(N * sizeof(double));
        qmax = malloc(N * sizeof(double));
        for (k=0; k<N; k++) {
                k3=k*3;

                qmin[k] = ((double*) qc -> data)[k];
                qmax[k] = qmin[k];

                for (i=0; i<3; i++) {
                        n = ((long*) neighbours -> data)[k3+i];
                        if (n >= 0) {
                                qn = ((double*) qc -> data)[n]; //Neighbour's centroid value

                                qmin[k] = min(qmin[k], qn);
                                qmax[k] = max(qmax[k], qn);
                        }
                        //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]);
                        //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]);
                }
        }

        // Call underlying routine
        _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);

        free(qmin);
        free(qmax);
        return Py_BuildValue("");
}



// Method table for python module
static struct PyMethodDef MethodTable[] = {
        {"limit", limit, METH_VARARGS, "Print out"},
        {"update", update, METH_VARARGS, "Print out"},
        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
        {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"},
        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
        {"extrapolate_second_order", extrapolate_second_order,
                METH_VARARGS, "Print out"},
        {"interpolate_from_vertices_to_edges",
                interpolate_from_vertices_to_edges,
                METH_VARARGS, "Print out"},
        {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"},            
        {NULL, NULL, 0, NULL}   // sentinel
};

// Module initialisation
void initquantity_ext(void){
  Py_InitModule("quantity_ext", MethodTable);

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