source: inundation/ga/storm_surge/pyvolution/quantity_ext.c @ 248

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


//FIXME: Should live in util_ext.c
double max(double x, double y) {  
  //Return maximum of two doubles
  
  if (x > y) return x;
  else return y;
}

double min(double x, double y) {  
  //Return minimum of two doubles
  
  if (x < y) return x;
  else return y;
}


void _limit(int N, double beta, double* qc, double* qv, double* qmin, double* qmax) { 

  int k, i, ki;
  double dq, phi, r;
  
  for (k=0; k<N; k++) {
  
    //Find the gradient limiter (phi) across vertices  
    phi = 1.0;
    for (i=0; i<3; i++) {    
      ki = k*3+i;
      r = 1.0;
      
      dq = qv[ki] - qc[k];    //Delta between vertex and centroid values
    
      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
  
      phi = min( min(r*beta, 1.0), phi);    
    }
    //Then update using phi limiter
    for (i=0; i<3; i++) {    
      qv[ki] = qc[k] + phi*(qv[ki] - qc[k]); //FIXME: Could precompute dq (3 elements)
    }
  }

     
  /*            
    #Diffences between centroids and maxima/minima
    dqmax = qmax - qc 
    dqmin = qmin - qc
        
    #Deltas between vertex and centroid values
    dq = zeros(qv.shape, Float)
    for i in range(3):
        dq[:,i] = qv[:,i] - qc

    #Phi limiter    
    for k in range(N):
        
        #Find the gradient limiter (phi) across vertices
        phi = 1.0
        for i in range(3):
            r = 1.0
            if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
            if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
            
            phi = min( min(r*beta, 1), phi )    

        #Then update using phi limiter
        for i in range(3):            
            qv[k,i] = qc[k] + phi*dq[k,i]
            
  */
  
}

// Gateway to Python
//FIXME: DOES NOT WORK YET
PyObject *limit(PyObject *self, PyObject *args) {
  
  PyObject *quantity;
  
  PyObject *Tmp;
  
  PyArrayObject 
    *qv, //Conserved quantities at vertices 
    *qc, //Conserved quantities at centroids
    *neighbours;
    
  int k, i, n, N;
  double beta; //Safety factor
  double *qmin, *qmax, qn;
  
  // Convert Python arguments to C  
  if (!PyArg_ParseTuple(args, "O", &quantity)) 
    return NULL;
  
  //Get safety factor beta
  Tmp = PyObject_GetAttrString(quantity, "beta"); 
  beta = PyFloat_AsDouble(Tmp);  
  Py_DECREF(Tmp);         
  
  
  qc = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); 
  qv = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values");   
  
  Tmp = PyObject_GetAttrString(quantity, "domain");     
  neighbours = (PyArrayObject*) PyObject_GetAttrString(Tmp, "neighbours");     
  Py_DECREF(Tmp);           
  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++) {
    qmax[k] = qmin[k] = ((double*) qc -> data)[k];
    for (i=0; i<3; i++) {
      n = ((int*) neighbours -> data)[k*3+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);
      }
    }
  }
  
  // Call underlying routine
  //_limit(N, beta, (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"},    
  {NULL, NULL, 0, NULL}   /* sentinel */
};

// Module initialisation   
void initquantity_ext(void){
  Py_InitModule("quantity_ext", MethodTable);
  
  import_array();     //Necessary for handling of NumPY structures  
}