source: anuga_work/development/flow_1d/utilities/util_ext.h @ 7830

Last change on this file since 7830 was 7830, checked in by steve, 14 years ago

Moving channel code to numpy

File size: 8.4 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
264  //printf("B = %p\n",(void*)B);
265  if (!B) {
266    printf("util_ext.h: get_consecutive_array could not obtain python object");
267    printf(" %s\n",name);
268    fflush(stdout);
269    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain python object");
270    return NULL;
271  }     
272 
273  //Convert to consecutive array
274  A = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*) B, 
275                                                    B -> descr -> type, 0, 0);
276 
277  Py_DECREF(B); //FIXME: Is this really needed??
278 
279  if (!A) {
280    printf("util_ext.h: get_consecutive_array could not obtain array object");
281    printf(" %s \n",name);
282    fflush(stdout);
283    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
284    return NULL;
285  }
286
287
288  return A;
289}
290
291double get_python_double(PyObject *O, char *name) {
292  PyObject *TObject;
293  #define BUFFER_SIZE 80
294  char buf[BUFFER_SIZE];
295  double tmp;
296  int n;
297 
298
299  //Get double from attribute
300  TObject = PyObject_GetAttrString(O, name);
301  if (!TObject) {
302        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_double could not obtain double %s.\n", name);
303        //printf("name = %s",name);
304    PyErr_SetString(PyExc_RuntimeError, buf);
305
306    return 0.0;
307  } 
308 
309  tmp = PyFloat_AsDouble(TObject);
310 
311  Py_DECREF(TObject);
312 
313  return tmp;
314}
315
316
317
318
319int get_python_integer(PyObject *O, char *name) {
320  PyObject *TObject;
321  #define BUFFER_SIZE 80
322  char buf[BUFFER_SIZE];
323  long tmp;
324  int n;
325 
326
327  //Get double from attribute
328  TObject = PyObject_GetAttrString(O, name);
329  if (!TObject) {
330        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_integer could not obtain double %s.\n", name);
331        //printf("name = %s",name);
332    PyErr_SetString(PyExc_RuntimeError, buf);
333    return 0;
334  } 
335 
336  tmp = PyInt_AsLong(TObject);
337 
338  Py_DECREF(TObject);
339 
340  return tmp;
341}
342
343
344PyObject *get_python_object(PyObject *O, char *name) {
345  PyObject *Oout;
346
347  Oout = PyObject_GetAttrString(O, name);
348  if (!Oout) {
349    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
350    return NULL;
351  }
352
353  return Oout;
354
355}
356
357
358// check that numpy array objects are C contiguous memory
359#define CHECK_C_CONTIG(varname) if (!PyArray_ISCONTIGUOUS(varname)) { \
360                                    char msg[1024]; \
361                                    sprintf(msg, \
362                                            "%s(): file %s, line %d: " \
363                                            "'%s' object is not C contiguous memory", \
364                                             __func__, __FILE__, __LINE__, #varname); \
365                                    PyErr_SetString(PyExc_RuntimeError, msg); \
366                                    return NULL; \
367                                }
Note: See TracBrowser for help on using the repository browser.