Ignore:
Timestamp:
Sep 1, 2004, 4:16:18 AM (20 years ago)
Author:
ole
Message:

Compute_gradients c implementation

Location:
inundation/ga/storm_surge/pyvolution
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r260 r261  
    338338            k0, k1, k2 = surrogate_neighbours[k,:]
    339339
     340            #Get data       
    340341            q0 = centroid_values[k0]
    341342            q1 = centroid_values[k1]
    … …  
    450451    #Replace python version with c implementations
    451452       
    452     from quantity_ext import limit #, compute_gradients #, extrapolate_second_order       
     453    from quantity_ext import limit, compute_gradients #, extrapolate_second_order       

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

    r260 r261  
    1818#include "quantity_ext.h"
    1919
     20
     21
     22int _compute_gradients(int N,
     23                        double* centroids,
     24                        double* centroid_values,
     25                        int* number_of_boundaries,
     26                        int* surrogate_neighbours,
     27                        double* a,
     28                        double* b) {
     29                       
     30  int i, k, k0, k1, k2, index3;
     31  double x0, x1, x2, y0, y1, y2, q0, q1, q2, det;
     32 
     33 
     34  for (k=0; k<N; k++) {
     35    index3 = 3*k;
     36   
     37    if (number_of_boundaries[k] < 2) {
     38      //Two or three true neighbours
     39
     40      //Get indices of neighbours (or self when used as surrogate)         
     41      //k0, k1, k2 = surrogate_neighbours[k,:]
     42     
     43      k0 = surrogate_neighbours[index3 + 0];
     44      k1 = surrogate_neighbours[index3 + 1];
     45      k2 = surrogate_neighbours[index3 + 2];           
     46      if (k0 == k1 || k1 == k2) return -1;     
     47
     48      //Get data
     49      q0 = centroid_values[k0];
     50      q1 = centroid_values[k1];
     51      q2 = centroid_values[k2];                   
     52
     53      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     54      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
     55      x2 = centroids[k2*2]; y2 = centroids[k2*2+1];             
     56
     57      //Gradient
     58      _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
     59     
     60    } else if (number_of_boundaries[k] == 2) {
     61      //One true neighbour
     62
     63      //#Get index of the one neighbour
     64      i=0; k0 = k;
     65      while (i<3 && k0==k) {
     66        k0 = surrogate_neighbours[index3 + i];
     67        i++;
     68      }
     69      if (k0 == k) return -1;
     70     
     71      k1 = k; //self
     72
     73      //Get data
     74      q0 = centroid_values[k0];
     75      q1 = centroid_values[k1];
     76           
     77      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     78      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
     79
     80      //Gradient
     81      det = x0*y1 - x1*y0;
     82      if (det != 0.0) {
     83        a[k] = (y1*q0 - y0*q1)/det;
     84        b[k] = (x0*q1 - x1*q0)/det;
     85      }
     86    }
     87    //    else:
     88    //        #No true neighbours -     
     89    //        #Fall back to first order scheme
     90    //        pass
     91  }
     92  return 0;
     93}           
     94
     95
     96
     97
    2098/////////////////////////////////////////////////
    2199// Gateways to Python
    … …  
    23101
    24102
    25 
    26 
    27 /*
    28 
    29 def compute_gradients(quantity):                   
    30     """Compute gradients of triangle surfaces defined by centroids of
    31     neighbouring volumes.
    32     If one edge is on the boundary, use own centroid as neighbour centroid.
    33     If two or more are on the boundary, fall back to first order scheme.
    34     """
    35 
    36     from Numeric import zeros, Float
    37     from util import gradient
    38    
    39     centroids = quantity.domain.centroids
    40     surrogate_neighbours = quantity.domain.surrogate_neighbours   
    41     centroid_values = quantity.centroid_values   
    42     number_of_boundaries = quantity.domain.number_of_boundaries     
    43    
    44     N = centroid_values.shape[0]
    45 
    46     a = zeros(N, Float)
    47     b = zeros(N, Float)
    48    
    49     for k in range(N):
    50         if number_of_boundaries[k] < 2:
    51             #Two or three true neighbours
    52 
    53             #Get indices of neighbours (or self when used as surrogate)     
    54             k0, k1, k2 = surrogate_neighbours[k,:]
    55 
    56             q0 = centroid_values[k0]
    57             q1 = centroid_values[k1]
    58             q2 = centroid_values[k2]                   
    59 
    60             x0, y0 = centroids[k0] #V0 centroid
    61             x1, y1 = centroids[k1] #V1 centroid
    62             x2, y2 = centroids[k2] #V2 centroid
    63 
    64             #Gradient
    65             a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    66        
    67         elif number_of_boundaries[k] == 2:
    68             #One true neighbour
    69 
    70             #Get index of the one neighbour
    71             for k0 in surrogate_neighbours[k,:]:
    72                 if k0 != k: break
    73             assert k0 != k
    74 
    75             k1 = k  #self
    76 
    77             #Get data
    78             q0 = centroid_values[k0]
    79             q1 = centroid_values[k1]
    80 
    81             x0, y0 = centroids[k0] #V0 centroid
    82             x1, y1 = centroids[k1] #V1 centroid       
    83 
    84             #Gradient
    85             det = x0*y1 - x1*y0
    86             if det != 0.0:
    87                 a[k] = (y1*q0 - y0*q1)/det
    88                 b[k] = (x0*q1 - x1*q0)/det
    89 
    90         else:
    91             #No true neighbours -       
    92             #Fall back to first order scheme
    93             pass
    94        
    95    
    96     return a, b
    97        
    98 */
    99 
    100 
    101 
    102 
    103 
    104 
    105 
    106 
    107 
     103PyObject *compute_gradients(PyObject *self, PyObject *args) {
     104 
     105  PyObject *quantity, *domain;
     106  PyArrayObject
     107    *centroids,            //Coordinates at centroids
     108    *centroid_values,      //Values at centroids   
     109    *number_of_boundaries, //Number of boundaries for each triangle
     110    *surrogate_neighbours, //True neighbours or - if one missing - self
     111    *a, *b;                //Return values
     112   
     113  int dimensions[1], N, err;
     114 
     115  // Convert Python arguments to C 
     116  if (!PyArg_ParseTuple(args, "O", &quantity))
     117    return NULL;
     118
     119  domain = PyObject_GetAttrString(quantity, "domain");     
     120  if (!domain)
     121    return NULL; 
     122
     123  //Get pertinent variables
     124  centroids = (PyArrayObject*) PyObject_GetAttrString(domain, "centroids");
     125  if (!centroids) return NULL;   
     126 
     127  centroid_values = (PyArrayObject*)
     128    PyObject_GetAttrString(quantity, "centroid_values");   
     129  if (!centroid_values) return NULL;       
     130 
     131  surrogate_neighbours = (PyArrayObject*)
     132    PyObject_GetAttrString(domain, "surrogate_neighbours");
     133  if (!surrogate_neighbours) return NULL;       
     134 
     135  number_of_boundaries = (PyArrayObject*)
     136    PyObject_GetAttrString(domain, "number_of_boundaries");           
     137  if (!number_of_boundaries) return NULL;           
     138 
     139  N = centroid_values -> dimensions[0];
     140   
     141  //Release
     142  Py_DECREF(domain);     
     143         
     144  //Allocate space for return vectors a and b (don't DECREF)
     145  dimensions[0] = N;
     146  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     147  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 
     148       
     149
     150 
     151  err = _compute_gradients(N,
     152                           (double*) centroids -> data,
     153                           (double*) centroid_values -> data,
     154                           (int*) number_of_boundaries -> data,
     155                           (int*) surrogate_neighbours -> data,
     156                           (double*) a -> data,
     157                           (double*) b -> data);     
     158                           
     159  if (err != 0) {
     160    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
     161    return NULL;
     162  }
     163                     
     164 
     165//Return values
     166  return Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
     167}
    108168
    109169
    110170/*
    111171PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
     172
     173
     174
     175PyObject *compute_gradients(PyObject *self, PyObject *args) {
     176 
     177  PyObject *quantity, *domain;
     178  PyArrayObject
     179    *centroids,            //Coordinates at centroids
     180    *centroid_values,      //Values at centroids   
     181    *number_of_boundaries, //Number of boundaries for each triangle
     182    *surrogate_neighbours, //True neighbours or - if one missing - self
     183    *a, *b;                //Return values
     184   
     185  int dimensions[1], N, err;
     186 
     187  // Convert Python arguments to C 
     188  if (!PyArg_ParseTuple(args, "O", &quantity))
     189    return NULL;
     190
     191  domain = PyObject_GetAttrString(quantity, "domain");     
     192  if (!domain)
     193    return NULL; 
     194
     195  //Get pertinent variables
     196  centroids = (PyArrayObject*) PyObject_GetAttrString(domain, "centroids");
     197  if (!centroids) return NULL;   
     198 
     199  centroid_values = (PyArrayObject*)
     200    PyObject_GetAttrString(quantity, "centroid_values");   
     201  if (!centroid_values) return NULL;       
     202 
     203  surrogate_neighbours = (PyArrayObject*)
     204    PyObject_GetAttrString(domain, "surrogate_neighbours");
     205  if (!surrogate_neighbours) return NULL;       
     206 
     207  number_of_boundaries = (PyArrayObject*)
     208    PyObject_GetAttrString(domain, "number_of_boundaries");           
     209  if (!number_of_boundaries) return NULL;           
     210 
     211  N = centroid_values -> dimensions[0];
     212   
     213  //Release
     214  Py_DECREF(domain);     
     215         
     216  //Allocate space for return vectors a and b (don't DECREF)
     217  dimensions[0] = N;
     218  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     219  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 
     220       
     221
     222 
     223  err = _compute_gradients(N,
     224                           (double*) centroids -> data,
     225                           (double*) centroid_values -> data,
     226                           (int*) number_of_boundaries -> data,
     227                           (int*) surrogate_neighbours -> data,
     228                           (double*) a -> data,
     229                           (double*) b -> data);     
     230                           
     231  if (err != 0) {
     232    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
     233    return NULL;
     234  }
     235
    112236
    113237  PyObject *quantity, *domain, *Tmp;
    … …  
    234358static struct PyMethodDef MethodTable[] = {
    235359  {"limit", limit, METH_VARARGS, "Print out"},   
     360  {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},     
    236361  {NULL, NULL, 0, NULL}   /* sentinel */
    237362};
Note: See TracChangeset for help on using the changeset viewer.