source: inundation/pyvolution/quantity_ext.c @ 2689

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