source: branches/numpy/anuga/utilities/util_ext.h @ 7091

Last change on this file since 7091 was 6689, checked in by rwilson, 15 years ago

Back-merge from Numeric trunk to numpy branch.

File size: 8.1 KB
Line 
1// Python - C extension for finite_volumes util 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 util.py
8//
9//
10// Ole Nielsen, GA 2004
11       
12#include "Python.h"     
13#include "numpy/arrayobject.h"
14#include "math.h"
15
16
17double max(double x, double y) { 
18  //Return maximum of two doubles
19 
20  if (x > y) return x;
21  else return y;
22}
23
24
25double min(double x, double y) { 
26  //Return minimum of two doubles
27 
28  if (x < y) return x;
29  else return y;
30}
31
32
33int _gradient(double x0, double y0, 
34              double x1, double y1, 
35              double x2, double y2, 
36              double q0, double q1, double q2, 
37              double *a, double *b) {
38             
39  /*Compute gradient (a,b) based on three points (x0,y0), (x1,y1) and (x2,y2)
40  with values q0, q1 and q2.
41 
42  Extrapolation formula (q0 is selected as an arbitrary origin)
43    q(x,y) = q0 + a*(x-x0) + b*(y-y0)                    (1)
44 
45  Substituting the known values for q1 and q2 into (1) yield the
46  equations for a and b
47 
48      q1-q0 = a*(x1-x0) + b*(y1-y0)                      (2)
49      q2-q0 = a*(x2-x0) + b*(y2-y0)                      (3)     
50     
51  or in matrix form
52 
53  /               \  /   \   /       \ 
54  |  x1-x0  y1-y0 |  | a |   | q1-q0 |
55  |               |  |   | = |       |
56  |  x2-x0  y2-y0 |  | b |   | q2-q0 |
57  \               /  \   /   \       /
58   
59  which is solved using the standard determinant technique   
60     
61  */
62             
63
64  double det;
65 
66  det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0);
67
68  *a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0);
69  *a /= det;
70
71  *b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0);
72  *b /= det;
73
74  return 0;
75}
76
77
78int _gradient2(double x0, double y0, 
79               double x1, double y1, 
80               double q0, double q1, 
81               double *a, double *b) {
82  /*Compute gradient (a,b) between two points (x0,y0) and (x1,y1)
83  with values q0 and q1 such that the plane is constant in the direction
84  orthogonal to (x1-x0, y1-y0).
85 
86  Extrapolation formula
87    q(x,y) = q0 + a*(x-x0) + b*(y-y0)                    (1)
88 
89  Substituting the known values for q1 into (1) yields an
90  under determined  equation for a and b
91      q1-q0 = a*(x1-x0) + b*(y1-y0)                      (2)
92     
93     
94  Now add the additional requirement that the gradient in the direction
95  orthogonal to (x1-x0, y1-y0) should be zero. The orthogonal direction
96  is given by the vector (y0-y1, x1-x0).
97 
98  Define the point (x2, y2) = (x0 + y0-y1, y0 + x1-x0) on the orthognal line.
99  Then we know that the corresponding value q2 should be equal to q0 in order
100  to obtain the zero gradient, hence applying (1) again   
101      q0 = q2 = q(x2, y2) = q0 + a*(x2-x0) + b*(y2-y0)
102                          = q0 + a*(x0 + y0-y1-x0) + b*(y0 + x1-x0 - y0)
103                          = q0 + a*(y0-y1) + b*(x1-x0)
104                         
105  leads to the orthogonality constraint
106     a*(y0-y1) + b*(x1-x0) = 0                           (3)
107     
108  which closes the system and yields
109 
110  /               \  /   \   /       \ 
111  |  x1-x0  y1-y0 |  | a |   | q1-q0 |
112  |               |  |   | = |       |
113  |  y0-y1  x1-x0 |  | b |   |   0   |
114  \               /  \   /   \       /
115   
116  which is solved using the standard determinant technique   
117     
118  */
119
120  double det, xx, yy, qq;
121 
122  xx = x1-x0;
123  yy = y1-y0;
124  qq = q1-q0;
125   
126  det = xx*xx + yy*yy;  //FIXME  catch det == 0
127  *a = xx*qq/det;
128  *b = yy*qq/det;
129       
130  return 0;
131}
132
133
134void _limit_old(int N, double beta, double* qc, double* qv, 
135            double* qmin, double* qmax) { 
136
137  //N are the number of elements
138  int k, i, k3;
139  double dq, dqa[3], phi, r;
140 
141  //printf("INSIDE\n");
142  for (k=0; k<N; k++) {
143    k3 = k*3;
144   
145    //Find the gradient limiter (phi) across vertices 
146    phi = 1.0;
147    for (i=0; i<3; i++) {   
148      r = 1.0;
149     
150      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
151      dqa[i] = dq;              //Save dq for use in the next loop
152     
153      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
154      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;     
155 
156 
157      phi = min( min(r*beta, 1.0), phi);   
158    }
159   
160    //Then update using phi limiter
161    for (i=0; i<3; i++) {   
162      qv[k3+i] = qc[k] + phi*dqa[i];
163    }
164  }
165}
166
167
168void  print_double_array(char* name, double* array, int n, int m){
169
170    int k,i,km;
171
172    printf("%s = [",name);
173    for (k=0; k<n; k++){
174        km = m*k;
175        printf("[");
176        for (i=0; i<m ; i++){
177            printf("%g ",array[km+i]);
178        }
179        if (k==(n-1))
180            printf("]");
181        else
182            printf("]\n");
183    }
184    printf("]\n");
185}
186
187void  print_int_array(char* name, int* array, int n, int m){
188
189    int k,i,km;
190
191    printf("%s = [",name);
192    for (k=0; k<n; k++){
193        km = m*k;
194        printf("[");
195        for (i=0; i<m ; i++){
196            printf("%i ",array[km+i]);
197        }
198        if (k==(n-1))
199            printf("]");
200        else
201            printf("]\n");
202    }
203    printf("]\n");
204}
205
206
207void  print_long_array(char* name, long* array, int n, int m){
208
209    int k,i,km;
210
211    printf("%s = [",name);
212    for (k=0; k<n; k++){
213        km = m*k;
214        printf("[");
215        for (i=0; i<m ; i++){
216          printf("%i ",(int) array[km+i]);
217        }
218        if (k==(n-1))
219            printf("]");
220        else
221            printf("]\n");
222    }
223    printf("]\n");
224}
225
226void print_numeric_array(PyArrayObject *x) { 
227  int i, j;
228  for (i=0; i<x->dimensions[0]; i++) { 
229    for (j=0; j<x->dimensions[1]; j++) {
230      printf("%f ", *(double*) (x->data + i*x->strides[0] + j*x->strides[1]));
231    }
232    printf("\n"); 
233  }
234  printf("\n");   
235}
236
237void print_numeric_vector(PyArrayObject *x) { 
238  int i;
239  for (i=0; i<x->dimensions[0]; i++) {
240    printf("%f ", *(double*) (x->data + i*x->strides[0])); 
241  }
242  printf("\n"); 
243}
244
245PyArrayObject *get_consecutive_array(PyObject *O, char *name) {
246  PyArrayObject *A, *B;
247 
248
249  //Get array object from attribute
250 
251  /*
252  //FIXME: THE TEST DOESN't WORK
253  printf("Err = %d\n", PyObject_HasAttrString(O, name));
254  if (PyObject_HasAttrString(O, name) == 1) {
255    B = (PyArrayObject*) PyObject_GetAttrString(O, name);
256    if (!B) return NULL;
257  } else {
258    return NULL;
259    }
260  */
261   
262  B = (PyArrayObject*) PyObject_GetAttrString(O, name);
263  if (!B) {
264    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
265    return NULL;
266  }     
267 
268  //Convert to consecutive array
269  A = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*) B, 
270                                                    B -> descr -> type, 0, 0);
271 
272  Py_DECREF(B); //FIXME: Is this really needed??
273 
274  if (!A) {
275    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
276    return NULL;
277  }
278  return A;
279}
280
281double get_python_double(PyObject *O, char *name) {
282  PyObject *TObject;
283  #define BUFFER_SIZE 80
284  char buf[BUFFER_SIZE];
285  double tmp;
286  int n;
287 
288
289  //Get double from attribute
290  TObject = PyObject_GetAttrString(O, name);
291  if (!TObject) {
292        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_double could not obtain double %s.\n", name);
293        //printf("name = %s",name);
294    PyErr_SetString(PyExc_RuntimeError, buf);
295
296    return 0.0;
297  } 
298 
299  tmp = PyFloat_AsDouble(TObject);
300 
301  Py_DECREF(TObject);
302 
303  return tmp;
304}
305
306
307
308
309int get_python_integer(PyObject *O, char *name) {
310  PyObject *TObject;
311  #define BUFFER_SIZE 80
312  char buf[BUFFER_SIZE];
313  long tmp;
314  int n;
315 
316
317  //Get double from attribute
318  TObject = PyObject_GetAttrString(O, name);
319  if (!TObject) {
320        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_integer could not obtain double %s.\n", name);
321        //printf("name = %s",name);
322    PyErr_SetString(PyExc_RuntimeError, buf);
323    return 0;
324  } 
325 
326  tmp = PyInt_AsLong(TObject);
327 
328  Py_DECREF(TObject);
329 
330  return tmp;
331}
332
333
334PyObject *get_python_object(PyObject *O, char *name) {
335  PyObject *Oout;
336
337  Oout = PyObject_GetAttrString(O, name);
338  if (!Oout) {
339    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
340    return NULL;
341  }
342
343  return Oout;
344
345}
346
347
348// check that numpy array objects are C contiguous memory
349#define CHECK_C_CONTIG(varname) if (!PyArray_ISCONTIGUOUS(varname)) { \
350                                    char msg[1024]; \
351                                    sprintf(msg, \
352                                            "%s(): file %s, line %d: " \
353                                            "'%s' object is not C contiguous memory", \
354                                             __func__, __FILE__, __LINE__, #varname); \
355                                    PyErr_SetString(PyExc_RuntimeError, msg); \
356                                    return NULL; \
357                                }
Note: See TracBrowser for help on using the repository browser.