Changeset 1156


Ignore:
Timestamp:
Mar 29, 2005, 3:53:45 PM (19 years ago)
Author:
steve
Message:

Moved zeus files to separate directory.

Hopefully fixed up logging to not casue error if no log.ini file exists in cureent directory

Location:
inundation/ga/storm_surge
Files:
5 added
2 deleted
5 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/cg_solve.py

    r1150 r1156  
    44import logging, logging.config
    55logger = logging.getLogger('cg_solve')
     6logger.setLevel(logging.DEBUG)
    67
    7 logging.config.fileConfig('log.ini')
     8try:
     9    logging.config.fileConfig('log.ini')
     10except:
     11    pass
    812
    9 #logger.setLevel(logging.DEBUG)
     13
     14
     15
     16#
    1017
    1118

  • inundation/ga/storm_surge/pyvolution/log.ini

    r1150 r1156  
    3838# Formats
    3939[formatter_form01]
    40 format=%(name)s %(levelname)s %(message)s
     40format=%(name)s_%(levelname)s: %(message)s
    4141datefmt=
    42 

  • inundation/ga/storm_surge/pyvolution/pyvolution.zpi

    r1065 r1156  
    8181    <file>interpolate_sww.py</file>
    8282    <file>least_squares.py</file>
     83    <file>log.ini</file>
    8384    <file>mesh.py</file>
    8485    <file>mesh_factory.py</file>

  • inundation/ga/storm_surge/pyvolution/quantity_ext.c

    r907 r1156  
    22//
    33// To compile (Python2.3):
    4 //  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O 
    5 //  gcc -shared util_ext.o  -o util_ext.so     
     4//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O
     5//  gcc -shared util_ext.o  -o util_ext.so
    66//
    77// See the module quantity.py
    … …  
    99//
    1010// Ole Nielsen, GA 2004
    11        
     11
    1212#include "Python.h"
    1313#include "Numeric/arrayobject.h"
    … …  
    1515
    1616//Shared code snippets
    17 #include "util_ext.h" 
     17#include "util_ext.h"
    1818//#include "quantity_ext.h" //Obsolete
    1919
    … …  
    2121
    2222int _compute_gradients(int N,
    23                         double* centroids, 
     23                        double* centroids,
    2424                        double* centroid_values,
    2525                        long* number_of_boundaries,
    … …  
    2727                        double* a,
    2828                        double* b) {
    29                        
     29
    3030  int i, k, k0, k1, k2, index3;
    3131  double x0, x1, x2, y0, y1, y2, q0, q1, q2, det;
    32  
     32
    3333
    3434  for (k=0; k<N; k++) {
    3535    index3 = 3*k;
    36    
     36
    3737    if (number_of_boundaries[k] < 2) {
    3838      //Two or three true neighbours
    3939
    40       //Get indices of neighbours (or self when used as surrogate)         
     40      //Get indices of neighbours (or self when used as surrogate)
    4141      //k0, k1, k2 = surrogate_neighbours[k,:]
    42      
     42
    4343      k0 = surrogate_neighbours[index3 + 0];
    4444      k1 = surrogate_neighbours[index3 + 1];
    45       k2 = surrogate_neighbours[index3 + 2];           
    46 
    47 
    48       if (k0 == k1 || k1 == k2) return -1;     
     45      k2 = surrogate_neighbours[index3 + 2];
     46
     47
     48      if (k0 == k1 || k1 == k2) return -1;
    4949
    5050      //Get data
    5151      q0 = centroid_values[k0];
    5252      q1 = centroid_values[k1];
    53       q2 = centroid_values[k2];                   
    54 
    55       x0 = centroids[k0*2]; y0 = centroids[k0*2+1]; 
    56       x1 = centroids[k1*2]; y1 = centroids[k1*2+1]; 
    57       x2 = centroids[k2*2]; y2 = centroids[k2*2+1];             
     53      q2 = centroid_values[k2];
     54
     55      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     56      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
     57      x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
    5858
    5959      //Gradient
    6060      _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
    61      
     61
    6262    } else if (number_of_boundaries[k] == 2) {
    6363      //One true neighbour
    … …  
    7070      }
    7171      if (k0 == k) return -1;
    72      
     72
    7373      k1 = k; //self
    7474
    … …  
    7676      q0 = centroid_values[k0];
    7777      q1 = centroid_values[k1];
    78            
    79       x0 = centroids[k0*2]; y0 = centroids[k0*2+1]; 
    80       x1 = centroids[k1*2]; y1 = centroids[k1*2+1]; 
     78
     79      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     80      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
    8181
    8282      //Gradient
    … …  
    8686        b[k] = (x0*q1 - x1*q0)/det;
    8787      }
    88     } 
     88    }
    8989    //    else:
    90     //        #No true neighbours -     
     90    //        #No true neighbours -
    9191    //        #Fall back to first order scheme
    9292    //        pass
    9393  }
    9494  return 0;
    95 }           
     95}
    9696
    9797
    9898int _extrapolate(int N,
    99                  double* centroids, 
    100                  double* centroid_values,               
     99                 double* centroids,
     100                 double* centroid_values,
    101101                 double* vertex_coordinates,
    102102                 double* vertex_values,
    103103                 double* a,
    104104                 double* b) {
    105  
     105
    106106  int k, k2, k3, k6;
    107107  double x, y, x0, y0, x1, y1, x2, y2;
    108                      
     108
    109109  for (k=0; k<N; k++) {
    110110    k6 = 6*k;
    111     k3 = 3*k;   
     111    k3 = 3*k;
    112112    k2 = 2*k;
    113        
    114     //Centroid coordinates           
    115     x = centroids[k2]; y = centroids[k2+1]; 
     113
     114    //Centroid coordinates
     115    x = centroids[k2]; y = centroids[k2+1];
    116116
    117117    //vertex coordinates
    118118    //x0, y0, x1, y1, x2, y2 = X[k,:]
    119119    x0 = vertex_coordinates[k6 + 0];
    120     y0 = vertex_coordinates[k6 + 1];   
     120    y0 = vertex_coordinates[k6 + 1];
    121121    x1 = vertex_coordinates[k6 + 2];
    122     y1 = vertex_coordinates[k6 + 3];       
     122    y1 = vertex_coordinates[k6 + 3];
    123123    x2 = vertex_coordinates[k6 + 4];
    124     y2 = vertex_coordinates[k6 + 5];           
     124    y2 = vertex_coordinates[k6 + 5];
    125125
    126126    //Extrapolate
    … …  
    128128    vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
    129129    vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
    130    
     130
    131131  }
    132132  return 0;
    … …  
    139139                 double* vertex_values,
    140140                 double* edge_values) {
    141  
     141
    142142  int k, k3;
    143143  double q0, q1, q2;
    144144
    145                      
    146   for (k=0; k<N; k++) {
    147     k3 = 3*k;   
    148    
     145
     146  for (k=0; k<N; k++) {
     147    k3 = 3*k;
     148
    149149    q0 = vertex_values[k3 + 0];
    150150    q1 = vertex_values[k3 + 1];
    151151    q2 = vertex_values[k3 + 2];
    152            
     152
    153153    //printf("%f, %f, %f\n", q0, q1, q2);
    154154    edge_values[k3 + 0] = 0.5*(q1+q2);
    155     edge_values[k3 + 1] = 0.5*(q0+q2); 
     155    edge_values[k3 + 1] = 0.5*(q0+q2);
    156156    edge_values[k3 + 2] = 0.5*(q0+q1);
    157   } 
     157  }
    158158  return 0;
    159159}
    160160
    161 int _update(int N, 
     161int _update(int N,
    162162            double timestep,
    163163            double* centroid_values,
    … …  
    166166  //Update centroid values based on values stored in
    167167  //explicit_update and semi_implicit_update as well as given timestep
    168  
    169  
     168
     169
    170170  int k;
    171171  double denominator, x;
    172172
    173  
    174   //Divide semi_implicit update by conserved quantity 
    175   for (k=0; k<N; k++) { 
     173
     174  //Divide semi_implicit update by conserved quantity
     175  for (k=0; k<N; k++) {
    176176    x = centroid_values[k];
    177177    if (x == 0.0) {
    … …  
    179179    } else {
    180180      semi_implicit_update[k] /= x;
    181     }   
    182   }             
    183  
    184  
     181    }
     182  }
     183
     184
    185185  //Explicit updates
    186186  for (k=0; k<N; k++) {
    187187    centroid_values[k] += timestep*explicit_update[k];
    188188  }
    189            
     189
    190190  //Semi implicit updates
    191   for (k=0; k<N; k++) { 
    192     denominator = 1.0 - timestep*semi_implicit_update[k]; 
    193    
     191  for (k=0; k<N; k++) {
     192    denominator = 1.0 - timestep*semi_implicit_update[k];
     193
    194194    if (denominator == 0.0) {
    195195      return -1;
    … …  
    201201  return 0;
    202202}
    203                      
    204  
     203
     204
    205205/////////////////////////////////////////////////
    206206// Gateways to Python
    207207PyObject *update(PyObject *self, PyObject *args) {
    208  
     208
    209209  PyObject *quantity;
    210   PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update; 
    211  
    212   double timestep; 
     210  PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
     211
     212  double timestep;
    213213  int N, err;
    214  
    215 
    216   // Convert Python arguments to C 
    217   if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep)) 
     214
     215
     216  // Convert Python arguments to C
     217  if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep))
    218218    return NULL;
    219219
    220220  centroid_values = get_consecutive_array(quantity, "centroid_values");
    221   explicit_update = get_consecutive_array(quantity, "explicit_update");   
    222   semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");     
    223 
    224   N = centroid_values -> dimensions[0]; 
    225  
    226  
     221  explicit_update = get_consecutive_array(quantity, "explicit_update");
     222  semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");
     223
     224  N = centroid_values -> dimensions[0];
     225
     226
    227227  err = _update(N, timestep,
    228228                (double*) centroid_values -> data,
    229229                (double*) explicit_update -> data,
    230                 (double*) semi_implicit_update -> data);               
    231                
    232                      
     230                (double*) semi_implicit_update -> data);
     231
     232
    233233  if (err != 0) {
    234     PyErr_SetString(PyExc_RuntimeError, 
     234    PyErr_SetString(PyExc_RuntimeError,
    235235                    "Zero division in semi implicit update - call Stephen :)");
    236236    return NULL;
    237   }                 
    238  
    239   //Release and return               
    240   Py_DECREF(centroid_values);   
    241   Py_DECREF(explicit_update); 
    242   Py_DECREF(semi_implicit_update);   
     237  }
     238
     239  //Release and return
     240  Py_DECREF(centroid_values);
     241  Py_DECREF(explicit_update);
     242  Py_DECREF(semi_implicit_update);
    243243
    244244  return Py_BuildValue("");
    … …  
    247247
    248248PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
    249  
     249
    250250  PyObject *quantity;
    251251  PyArrayObject *vertex_values, *edge_values;
    252    
     252
    253253  int N, err;
    254  
    255   // Convert Python arguments to C 
    256   if (!PyArg_ParseTuple(args, "O", &quantity)) 
     254
     255  // Convert Python arguments to C
     256  if (!PyArg_ParseTuple(args, "O", &quantity))
    257257    return NULL;
    258258
    259259  vertex_values = get_consecutive_array(quantity, "vertex_values");
    260   edge_values = get_consecutive_array(quantity, "edge_values"); 
    261  
    262   N = vertex_values -> dimensions[0]; 
    263  
     260  edge_values = get_consecutive_array(quantity, "edge_values");
     261
     262  N = vertex_values -> dimensions[0];
     263
    264264  err = _interpolate(N,
    265265                     (double*) vertex_values -> data,
    266266                     (double*) edge_values -> data);
    267                      
     267
    268268  if (err != 0) {
    269269    PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed");
    270270    return NULL;
    271   }                 
    272  
    273   //Release and return               
    274   Py_DECREF(vertex_values);   
    275   Py_DECREF(edge_values); 
     271  }
     272
     273  //Release and return
     274  Py_DECREF(vertex_values);
     275  Py_DECREF(edge_values);
    276276
    277277  return Py_BuildValue("");
    … …  
    280280
    281281PyObject *compute_gradients(PyObject *self, PyObject *args) {
    282  
     282
    283283  PyObject *quantity, *domain, *R;
    284   PyArrayObject 
     284  PyArrayObject
    285285    *centroids,            //Coordinates at centroids
    286     *centroid_values,      //Values at centroids   
     286    *centroid_values,      //Values at centroids
    287287    *number_of_boundaries, //Number of boundaries for each triangle
    288288    *surrogate_neighbours, //True neighbours or - if one missing - self
    289289    *a, *b;                //Return values
    290    
     290
    291291  int dimensions[1], N, err;
    292  
    293   // Convert Python arguments to C 
    294   if (!PyArg_ParseTuple(args, "O", &quantity)) 
    295     return NULL;
    296 
    297   domain = PyObject_GetAttrString(quantity, "domain");     
    298   if (!domain) 
    299     return NULL; 
     292
     293  // Convert Python arguments to C
     294  if (!PyArg_ParseTuple(args, "O", &quantity))
     295    return NULL;
     296
     297  domain = PyObject_GetAttrString(quantity, "domain");
     298  if (!domain)
     299    return NULL;
    300300
    301301  //Get pertinent variables
    … …  
    307307
    308308  N = centroid_values -> dimensions[0];
    309    
     309
    310310  //Release
    311   Py_DECREF(domain);     
    312          
    313   //Allocate space for return vectors a and b (don't DECREF) 
     311  Py_DECREF(domain);
     312
     313  //Allocate space for return vectors a and b (don't DECREF)
    314314  dimensions[0] = N;
    315315  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    316   b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 
    317        
    318 
    319  
     316  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     317
     318
     319
    320320  err = _compute_gradients(N,
    321                            (double*) centroids -> data,
    322                            (double*) centroid_values -> data,
    323                            (long*) number_of_boundaries -> data,
    324                            (long*) surrogate_neighbours -> data,
    325                            (double*) a -> data,
    326                            (double*) b -> data);     
    327                            
    328   if (err != 0) {
    329     PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
    330     return NULL;
    331   }
    332                      
    333   //Release                 
    334   Py_DECREF(centroids);   
    335   Py_DECREF(centroid_values);   
    336   Py_DECREF(number_of_boundaries);
    337   Py_DECREF(surrogate_neighbours);           
    338          
    339   //Build result, release and return
    340   R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
    341   Py_DECREF(a);     
    342   Py_DECREF(b);       
    343   return R;
    344 }
    345 
    346 
    347 
    348 PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
    349  
    350   PyObject *quantity, *domain;
    351   PyArrayObject
    352     *centroids,            //Coordinates at centroids
    353     *centroid_values,      //Values at centroids   
    354     *vertex_coordinates,   //Coordinates at vertices           
    355     *vertex_values,        //Values at vertices       
    356     *number_of_boundaries, //Number of boundaries for each triangle
    357     *surrogate_neighbours, //True neighbours or - if one missing - self
    358     *a, *b;                //Gradients
    359    
    360   //int N, err;
    361   int dimensions[1], N, err; 
    362   //double *a, *b;  //Gradients
    363  
    364   // Convert Python arguments to C 
    365   if (!PyArg_ParseTuple(args, "O", &quantity))
    366     return NULL;
    367 
    368   domain = PyObject_GetAttrString(quantity, "domain");     
    369   if (!domain)
    370     return NULL; 
    371 
    372   //Get pertinent variables
    373   centroids = get_consecutive_array(domain, "centroid_coordinates"); 
    374   centroid_values = get_consecutive_array(quantity, "centroid_values");
    375   surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
    376   number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
    377   vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");   
    378   vertex_values = get_consecutive_array(quantity, "vertex_values");     
    379  
    380   N = centroid_values -> dimensions[0];
    381    
    382   //Release
    383   Py_DECREF(domain);     
    384          
    385   //Allocate space for return vectors a and b (don't DECREF)
    386   dimensions[0] = N;
    387   a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    388   b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 
    389 
    390   //FIXME: Odd that I couldn't use normal arrays         
    391   //Allocate space for return vectors a and b (don't DECREF)
    392   //a = (double*) malloc(N * sizeof(double));
    393   //if (!a) return NULL; 
    394   //b = (double*) malloc(N * sizeof(double)); 
    395   //if (!b) return NULL;
    396 
    397  
    398   err = _compute_gradients(N,
    399                            (double*) centroids -> data,
     321                           (double*) centroids -> data,
    400322                           (double*) centroid_values -> data,
    401323                           (long*) number_of_boundaries -> data,
    … …  
    403325                           (double*) a -> data,
    404326                           (double*) b -> data);
    405                            
     327
    406328  if (err != 0) {
    407329    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
    408330    return NULL;
    409331  }
    410  
    411   err = _extrapolate(N,
    412                      (double*) centroids -> data,
    413                      (double*) centroid_values -> data,             
     332
     333  //Release
     334  Py_DECREF(centroids);
     335  Py_DECREF(centroid_values);
     336  Py_DECREF(number_of_boundaries);
     337  Py_DECREF(surrogate_neighbours);
     338
     339  //Build result, release and return
     340  R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
     341  Py_DECREF(a);
     342  Py_DECREF(b);
     343  return R;
     344}
     345
     346
     347
     348PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
     349
     350  PyObject *quantity, *domain;
     351  PyArrayObject
     352    *centroids,            //Coordinates at centroids
     353    *centroid_values,      //Values at centroids
     354    *vertex_coordinates,   //Coordinates at vertices
     355    *vertex_values,        //Values at vertices
     356    *number_of_boundaries, //Number of boundaries for each triangle
     357    *surrogate_neighbours, //True neighbours or - if one missing - self
     358    *a, *b;                //Gradients
     359
     360  //int N, err;
     361  int dimensions[1], N, err;
     362  //double *a, *b;  //Gradients
     363
     364  // Convert Python arguments to C
     365  if (!PyArg_ParseTuple(args, "O", &quantity))
     366    return NULL;
     367
     368  domain = PyObject_GetAttrString(quantity, "domain");
     369  if (!domain)
     370    return NULL;
     371
     372  //Get pertinent variables
     373  centroids = get_consecutive_array(domain, "centroid_coordinates");
     374  centroid_values = get_consecutive_array(quantity, "centroid_values");
     375  surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
     376  number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
     377  vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
     378  vertex_values = get_consecutive_array(quantity, "vertex_values");
     379
     380  N = centroid_values -> dimensions[0];
     381
     382  //Release
     383  Py_DECREF(domain);
     384
     385  //Allocate space for return vectors a and b (don't DECREF)
     386  dimensions[0] = N;
     387  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     388  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     389
     390  //FIXME: Odd that I couldn't use normal arrays
     391  //Allocate space for return vectors a and b (don't DECREF)
     392  //a = (double*) malloc(N * sizeof(double));
     393  //if (!a) return NULL;
     394  //b = (double*) malloc(N * sizeof(double));
     395  //if (!b) return NULL;
     396
     397
     398  err = _compute_gradients(N,
     399                           (double*) centroids -> data,
     400                           (double*) centroid_values -> data,
     401                           (long*) number_of_boundaries -> data,
     402                           (long*) surrogate_neighbours -> data,
     403                           (double*) a -> data,
     404                           (double*) b -> data);
     405
     406  if (err != 0) {
     407    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
     408    return NULL;
     409  }
     410
     411  err = _extrapolate(N,
     412                     (double*) centroids -> data,
     413                     (double*) centroid_values -> data,
    414414                     (double*) vertex_coordinates -> data,
    415415                     (double*) vertex_values -> data,
    416416                     (double*) a -> data,
    417                      (double*) b -> data);                   
    418                      
    419                      
     417                     (double*) b -> data);
     418
     419
    420420  if (err != 0) {
    421     PyErr_SetString(PyExc_RuntimeError, 
     421    PyErr_SetString(PyExc_RuntimeError,
    422422                    "Internal function _extrapolate failed");
    423423    return NULL;
    424   }                 
    425                
    426  
    427 
    428   //Release 
    429   Py_DECREF(centroids);   
    430   Py_DECREF(centroid_values);   
    431   Py_DECREF(number_of_boundaries); 
     424  }
     425
     426
     427
     428  //Release
     429  Py_DECREF(centroids);
     430  Py_DECREF(centroid_values);
     431  Py_DECREF(number_of_boundaries);
    432432  Py_DECREF(surrogate_neighbours);
    433   Py_DECREF(vertex_coordinates);                   
    434   Py_DECREF(vertex_values);                 
    435   Py_DECREF(a);                   
    436   Py_DECREF(b);                     
    437  
    438   return Py_BuildValue("");   
    439 }   
     433  Py_DECREF(vertex_coordinates);
     434  Py_DECREF(vertex_values);
     435  Py_DECREF(a);
     436  Py_DECREF(b);
     437
     438  return Py_BuildValue("");
     439}
    440440
    441441
    442442
    443443PyObject *limit(PyObject *self, PyObject *args) {
    444  
     444
    445445  PyObject *quantity, *domain, *Tmp;
    446   PyArrayObject 
    447     *qv, //Conserved quantities at vertices 
     446  PyArrayObject
     447    *qv, //Conserved quantities at vertices
    448448    *qc, //Conserved quantities at centroids
    449449    *neighbours;
    450    
     450
    451451  int k, i, n, N, k3;
    452452  double beta_w; //Safety factor
    453453  double *qmin, *qmax, qn;
    454  
    455   // Convert Python arguments to C 
    456   if (!PyArg_ParseTuple(args, "O", &quantity)) 
    457     return NULL;
    458 
    459   domain = PyObject_GetAttrString(quantity, "domain");     
    460   if (!domain) 
    461     return NULL; 
    462    
     454
     455  // Convert Python arguments to C
     456  if (!PyArg_ParseTuple(args, "O", &quantity))
     457    return NULL;
     458
     459  domain = PyObject_GetAttrString(quantity, "domain");
     460  if (!domain)
     461    return NULL;
     462
    463463  //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
    464464  neighbours = get_consecutive_array(domain, "neighbours");
    465465
    466466  //Get safety factor beta_w
    467   Tmp = PyObject_GetAttrString(domain, "beta_w"); 
    468   if (!Tmp) 
    469     return NULL;
    470      
    471   beta_w = PyFloat_AsDouble(Tmp); 
    472  
    473   Py_DECREF(Tmp);   
    474   Py_DECREF(domain);     
    475 
    476 
    477   qc = get_consecutive_array(quantity, "centroid_values"); 
    478   qv = get_consecutive_array(quantity, "vertex_values");   
    479    
    480  
     467  Tmp = PyObject_GetAttrString(domain, "beta_w");
     468  if (!Tmp)
     469    return NULL;
     470
     471  beta_w = PyFloat_AsDouble(Tmp);
     472
     473  Py_DECREF(Tmp);
     474  Py_DECREF(domain);
     475
     476
     477  qc = get_consecutive_array(quantity, "centroid_values");
     478  qv = get_consecutive_array(quantity, "vertex_values");
     479
     480
    481481  N = qc -> dimensions[0];
    482    
     482
    483483  //Find min and max of this and neighbour's centroid values
    484484  qmin = malloc(N * sizeof(double));
    485   qmax = malloc(N * sizeof(double)); 
     485  qmax = malloc(N * sizeof(double));
    486486  for (k=0; k<N; k++) {
    487487    k3=k*3;
    488    
     488
    489489    qmin[k] = ((double*) qc -> data)[k];
    490490    qmax[k] = qmin[k];
    491    
     491
    492492    for (i=0; i<3; i++) {
    493493      n = ((long*) neighbours -> data)[k3+i];
    494494      if (n >= 0) {
    495495        qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
    496          
     496
    497497        qmin[k] = min(qmin[k], qn);
    498498        qmax[k] = max(qmax[k], qn);
    … …  
    500500    }
    501501  }
    502  
     502
    503503  // Call underlying routine
    504504  _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
    505          
     505
    506506  free(qmin);
    507   free(qmax); 
    508   return Py_BuildValue(""); 
     507  free(qmax);
     508  return Py_BuildValue("");
    509509}
    510510
    … …  
    513513// Method table for python module
    514514static struct PyMethodDef MethodTable[] = {
    515   {"limit", limit, METH_VARARGS, "Print out"},   
    516   {"update", update, METH_VARARGS, "Print out"},     
    517   {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},     
    518   {"extrapolate_second_order", extrapolate_second_order, 
    519    METH_VARARGS, "Print out"},       
    520   {"interpolate_from_vertices_to_edges", 
     515  {"limit", limit, METH_VARARGS, "Print out"},
     516  {"update", update, METH_VARARGS, "Print out"},
     517  {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
     518  {"extrapolate_second_order", extrapolate_second_order,
     519   METH_VARARGS, "Print out"},
     520  {"interpolate_from_vertices_to_edges",
    521521   interpolate_from_vertices_to_edges,
    522    METH_VARARGS, "Print out"},           
     522   METH_VARARGS, "Print out"},
    523523  {NULL, NULL, 0, NULL}   // sentinel
    524524};
    525525
    526 // Module initialisation   
     526// Module initialisation
    527527void initquantity_ext(void){
    528528  Py_InitModule("quantity_ext", MethodTable);
    529  
    530   import_array();     //Necessary for handling of NumPY structures 
    531 }
    532 
     529
     530  import_array();     //Necessary for handling of NumPY structures
     531}
Note: See TracChangeset for help on using the changeset viewer.