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

Last change on this file since 628 was 587, checked in by steve, 20 years ago

testing sparse
File size: 15.2 KB
Line 
1// Python - C extension for quantity module.
2//
3// 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     
6//
7// See the module quantity.py
8//
9//
10// Ole Nielsen, GA 2004
11       
12#include "Python.h"
13#include "Numeric/arrayobject.h"
14#include "math.h"
15
16//Shared code snippets
17#include "util_ext.h"
18#include "quantity_ext.h"
19
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
96int _extrapolate(int N,
97                 double* centroids, 
98                 double* centroid_values,               
99                 double* vertex_coordinates,
100                 double* vertex_values,
101                 double* a,
102                 double* b) {
103 
104  int k, k2, k3, k6;
105  double x, y, x0, y0, x1, y1, x2, y2;
106                     
107  for (k=0; k<N; k++) {
108    k6 = 6*k;
109    k3 = 3*k;   
110    k2 = 2*k;
111       
112    //Centroid coordinates           
113    x = centroids[k2]; y = centroids[k2+1]; 
114
115    //vertex coordinates
116    //x0, y0, x1, y1, x2, y2 = X[k,:]
117    x0 = vertex_coordinates[k6 + 0];
118    y0 = vertex_coordinates[k6 + 1];   
119    x1 = vertex_coordinates[k6 + 2];
120    y1 = vertex_coordinates[k6 + 3];       
121    x2 = vertex_coordinates[k6 + 4];
122    y2 = vertex_coordinates[k6 + 5];           
123
124    //Extrapolate
125    vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
126    vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
127    vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
128   
129  }
130  return 0;
131}
132
133
134
135
136int _interpolate(int N,
137                 double* vertex_values,
138                 double* edge_values) {
139 
140  int k, k3;
141  double q0, q1, q2;
142
143                     
144  for (k=0; k<N; k++) {
145    k3 = 3*k;   
146   
147    q0 = vertex_values[k3 + 0];
148    q1 = vertex_values[k3 + 1];
149    q2 = vertex_values[k3 + 2];
150           
151    //printf("%f, %f, %f\n", q0, q1, q2);
152    edge_values[k3 + 0] = 0.5*(q1+q2);
153    edge_values[k3 + 1] = 0.5*(q0+q2); 
154    edge_values[k3 + 2] = 0.5*(q0+q1);
155  } 
156  return 0;
157}
158
159int _update(int N, 
160            double timestep,
161            double* centroid_values,
162            double* explicit_update,
163            double* semi_implicit_update) {
164  //Update centroid values based on values stored in
165  //explicit_update and semi_implicit_update as well as given timestep
166 
167 
168  int k;
169  double denominator, x;
170
171 
172  //Divide semi_implicit update by conserved quantity 
173  for (k=0; k<N; k++) { 
174    x = centroid_values[k];
175    if (x == 0.0) {
176      //FIXME: Is this right
177      semi_implicit_update[k] = 0.0;
178    } else {
179      semi_implicit_update[k] /= x;
180    }   
181  }             
182 
183 
184  //Explicit updates
185  for (k=0; k<N; k++) {
186    centroid_values[k] += timestep*explicit_update[k];
187  }
188           
189  //Semi implicit updates
190  for (k=0; k<N; k++) { 
191    denominator = 1.0 - timestep*semi_implicit_update[k]; 
192   
193    if (denominator == 0.0) {
194      return -1;
195    } else {
196      //Update conserved_quantities from semi implicit updates
197      centroid_values[k] /= denominator;
198    }
199  }
200  return 0;
201}
202                     
203 
204/////////////////////////////////////////////////
205// Gateways to Python
206PyObject *update(PyObject *self, PyObject *args) {
207 
208  PyObject *quantity;
209  PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update; 
210 
211  double timestep; 
212  int N, err;
213 
214
215  // Convert Python arguments to C 
216  if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep)) 
217    return NULL;
218
219  centroid_values = (PyArrayObject*) 
220    PyObject_GetAttrString(quantity, "centroid_values");           
221  if (!centroid_values) return NULL;               
222
223  explicit_update = (PyArrayObject*) 
224    PyObject_GetAttrString(quantity, "explicit_update");   
225  if (!explicit_update) return NULL;         
226 
227  semi_implicit_update = (PyArrayObject*) 
228    PyObject_GetAttrString(quantity, "semi_implicit_update");   
229  if (!semi_implicit_update) return NULL;           
230
231  N = centroid_values -> dimensions[0]; 
232 
233 
234  err = _update(N, timestep,
235                (double*) centroid_values -> data,
236                (double*) explicit_update -> data,
237                (double*) semi_implicit_update -> data);               
238               
239                     
240  if (err != 0) {
241    PyErr_SetString(PyExc_RuntimeError, 
242                    "Zero division in semi implicit update - call Stephen :)");
243    return NULL;
244  }                 
245 
246  //Release and return               
247  Py_DECREF(centroid_values);   
248  Py_DECREF(explicit_update); 
249  Py_DECREF(semi_implicit_update);   
250
251  return Py_BuildValue("");
252}
253
254
255PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
256 
257  PyObject *quantity;
258  PyArrayObject *vertex_values, *edge_values;
259   
260  int N, err;
261 
262  // Convert Python arguments to C 
263  if (!PyArg_ParseTuple(args, "O", &quantity)) 
264    return NULL;
265
266  vertex_values = (PyArrayObject*) 
267    PyObject_GetAttrString(quantity, "vertex_values");           
268  if (!vertex_values) return NULL;               
269
270  edge_values = (PyArrayObject*) 
271    PyObject_GetAttrString(quantity, "edge_values");   
272  if (!edge_values) return NULL;         
273
274  N = vertex_values -> dimensions[0]; 
275 
276  err = _interpolate(N,
277                     (double*) vertex_values -> data,
278                     (double*) edge_values -> data);
279                     
280  if (err != 0) {
281    PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed");
282    return NULL;
283  }                 
284 
285  //Release and return               
286  Py_DECREF(vertex_values);   
287  Py_DECREF(edge_values); 
288
289  return Py_BuildValue("");
290}
291
292
293PyObject *compute_gradients(PyObject *self, PyObject *args) {
294 
295  PyObject *quantity, *domain, *R;
296  PyArrayObject
297    *centroids,            //Coordinates at centroids
298    *centroid_values,      //Values at centroids   
299    *number_of_boundaries, //Number of boundaries for each triangle
300    *surrogate_neighbours, //True neighbours or - if one missing - self
301    *a, *b;                //Return values
302   
303  int dimensions[1], N, err;
304 
305  // Convert Python arguments to C 
306  if (!PyArg_ParseTuple(args, "O", &quantity)) 
307    return NULL;
308
309  domain = PyObject_GetAttrString(quantity, "domain");     
310  if (!domain) 
311    return NULL; 
312
313  //Get pertinent variables
314  centroids = (PyArrayObject*) 
315    PyObject_GetAttrString(domain, "centroid_coordinates"); 
316  if (!centroids) return NULL;   
317 
318  centroid_values = (PyArrayObject*) 
319    PyObject_GetAttrString(quantity, "centroid_values");   
320  if (!centroid_values) return NULL;       
321 
322  surrogate_neighbours = (PyArrayObject*) 
323    PyObject_GetAttrString(domain, "surrogate_neighbours");
324  if (!surrogate_neighbours) return NULL;       
325 
326  number_of_boundaries = (PyArrayObject*) 
327    PyObject_GetAttrString(domain, "number_of_boundaries");           
328  if (!number_of_boundaries) return NULL;           
329 
330  N = centroid_values -> dimensions[0];
331   
332  //Release
333  Py_DECREF(domain);     
334         
335  //Allocate space for return vectors a and b (don't DECREF)
336  dimensions[0] = N;
337  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
338  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 
339       
340
341 
342  err = _compute_gradients(N,
343                           (double*) centroids -> data, 
344                           (double*) centroid_values -> data,
345                           (int*) number_of_boundaries -> data,
346                           (int*) surrogate_neighbours -> data,
347                           (double*) a -> data,
348                           (double*) b -> data);     
349                           
350  if (err != 0) {
351    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
352    return NULL;
353  }
354                     
355  //Release                 
356  Py_DECREF(centroids);   
357  Py_DECREF(centroid_values);   
358  Py_DECREF(number_of_boundaries); 
359  Py_DECREF(surrogate_neighbours);           
360         
361  //Build result, release and return
362  R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); 
363  Py_DECREF(a);     
364  Py_DECREF(b);       
365  return R;
366}
367
368
369
370PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
371 
372  PyObject *quantity, *domain;
373  PyArrayObject
374    *centroids,            //Coordinates at centroids
375    *centroid_values,      //Values at centroids   
376    *vertex_coordinates,   //Coordinates at vertices           
377    *vertex_values,        //Values at vertices       
378    *number_of_boundaries, //Number of boundaries for each triangle
379    *surrogate_neighbours, //True neighbours or - if one missing - self
380    *a, *b;                //Gradients
381   
382  //int N, err;
383  int dimensions[1], N, err; 
384  //double *a, *b;  //Gradients
385 
386  // Convert Python arguments to C 
387  if (!PyArg_ParseTuple(args, "O", &quantity)) 
388    return NULL;
389
390  domain = PyObject_GetAttrString(quantity, "domain");     
391  if (!domain) 
392    return NULL; 
393
394  //Get pertinent variables
395  centroids = (PyArrayObject*) 
396    PyObject_GetAttrString(domain, "centroid_coordinates"); 
397  if (!centroids) return NULL;   
398 
399  centroid_values = (PyArrayObject*) 
400    PyObject_GetAttrString(quantity, "centroid_values");   
401  if (!centroid_values) return NULL;       
402 
403  surrogate_neighbours = (PyArrayObject*) 
404    PyObject_GetAttrString(domain, "surrogate_neighbours");
405  if (!surrogate_neighbours) return NULL;       
406 
407  number_of_boundaries = (PyArrayObject*) 
408    PyObject_GetAttrString(domain, "number_of_boundaries");           
409  if (!number_of_boundaries) return NULL;           
410 
411  vertex_coordinates = (PyArrayObject*) 
412    PyObject_GetAttrString(domain, "vertex_coordinates");           
413  if (!vertex_coordinates) return NULL;             
414 
415  vertex_values = (PyArrayObject*) 
416    PyObject_GetAttrString(quantity, "vertex_values");           
417  if (!vertex_values) return NULL;               
418
419 
420  /*
421  printf("In extrapolate C routine\n"); 
422  printf("d0=%d, d1=%d\n",
423         vertex_values -> dimensions[0],
424         vertex_values -> dimensions[1]);
425  */
426 
427  vertex_values = (PyArrayObject*) 
428    PyObject_GetAttrString(quantity, "vertex_values");           
429  if (!vertex_values) return NULL;               
430 
431  N = centroid_values -> dimensions[0];
432   
433  //Release
434  Py_DECREF(domain);     
435         
436  //Allocate space for return vectors a and b (don't DECREF)
437  dimensions[0] = N;
438  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
439  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 
440
441  //FIXME: Odd that I couldn't use normal arrays         
442  //Allocate space for return vectors a and b (don't DECREF)
443  //a = (double*) malloc(N * sizeof(double));
444  //if (!a) return NULL; 
445  //b = (double*) malloc(N * sizeof(double)); 
446  //if (!b) return NULL;
447
448 
449  err = _compute_gradients(N,
450                           (double*) centroids -> data, 
451                           (double*) centroid_values -> data,
452                           (int*) number_of_boundaries -> data,
453                           (int*) surrogate_neighbours -> data,
454                           (double*) a -> data,
455                           (double*) b -> data);
456                           
457  if (err != 0) {
458    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
459    return NULL;
460  }
461 
462  err = _extrapolate(N, 
463                     (double*) centroids -> data, 
464                     (double*) centroid_values -> data,             
465                     (double*) vertex_coordinates -> data,
466                     (double*) vertex_values -> data,
467                     (double*) a -> data,
468                     (double*) b -> data);                   
469  //a, b);
470                     
471                     
472  if (err != 0) {
473    PyErr_SetString(PyExc_RuntimeError, 
474                    "Internal function _extrapolate failed");
475    return NULL;
476  }                 
477               
478 
479
480  //Release 
481  Py_DECREF(centroids);   
482  Py_DECREF(centroid_values);   
483  Py_DECREF(number_of_boundaries); 
484  Py_DECREF(surrogate_neighbours);
485  Py_DECREF(vertex_coordinates);                   
486  Py_DECREF(vertex_values);                 
487  Py_DECREF(a);                   
488  Py_DECREF(b);                     
489  //free(a);     
490  //free(b);       
491 
492  return Py_BuildValue("");   
493}   
494
495
496
497PyObject *limit(PyObject *self, PyObject *args) {
498 
499  PyObject *quantity, *domain, *Tmp;
500  PyArrayObject
501    *qv, //Conserved quantities at vertices
502    *qc, //Conserved quantities at centroids
503    *neighbours;
504   
505  int k, i, n, N, k3;
506  double beta; //Safety factor
507  double *qmin, *qmax, qn;
508 
509  // Convert Python arguments to C 
510  if (!PyArg_ParseTuple(args, "O", &quantity)) 
511    return NULL;
512
513  domain = PyObject_GetAttrString(quantity, "domain");     
514  if (!domain) 
515    return NULL; 
516   
517  neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
518
519  //Get safety factor beta
520  Tmp = PyObject_GetAttrString(domain, "beta"); 
521  if (!Tmp) 
522    return NULL;
523     
524  beta = PyFloat_AsDouble(Tmp); 
525 
526  Py_DECREF(Tmp);   
527  Py_DECREF(domain);     
528
529  qc = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); 
530  qv = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values");   
531  N = qc -> dimensions[0];
532   
533  //Find min and max of this and neighbour's centroid values
534  qmin = malloc(N * sizeof(double));
535  qmax = malloc(N * sizeof(double)); 
536  for (k=0; k<N; k++) {
537    k3=k*3;
538   
539    qmin[k] = ((double*) qc -> data)[k];
540    qmax[k] = qmin[k];
541   
542    for (i=0; i<3; i++) {
543      n = ((int*) neighbours -> data)[k3+i];
544      if (n >= 0) {
545        qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
546         
547        qmin[k] = min(qmin[k], qn);
548        qmax[k] = max(qmax[k], qn);
549      }
550    }
551  }
552 
553  // Call underlying routine
554  _limit(N, beta, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
555         
556  free(qmin);
557  free(qmax); 
558  return Py_BuildValue(""); 
559}
560
561
562
563// Method table for python module
564static struct PyMethodDef MethodTable[] = {
565  {"limit", limit, METH_VARARGS, "Print out"},   
566  {"update", update, METH_VARARGS, "Print out"},     
567  {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},     
568  {"extrapolate_second_order", extrapolate_second_order, 
569   METH_VARARGS, "Print out"},       
570  {"interpolate_from_vertices_to_edges", 
571   interpolate_from_vertices_to_edges,
572   METH_VARARGS, "Print out"},           
573  {NULL, NULL, 0, NULL}   /* sentinel */
574};
575
576// Module initialisation   
577void initquantity_ext(void){
578  Py_InitModule("quantity_ext", MethodTable);
579 
580  import_array();     //Necessary for handling of NumPY structures 
581}
582
Note: See TracBrowser for help on using the repository browser.