source: anuga_work/development/anuga_1d/util_ext.h @ 6042

Last change on this file since 6042 was 6042, checked in by steve, 15 years ago

Found bug associated with changing max_speed_array to an integer

File size: 8.0 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 "Numeric/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
291
292
293double get_python_double(PyObject *O, char *name) {
294  PyObject *TObject;
295  #define BUFFER_SIZE 80
296  char buf[BUFFER_SIZE];
297  double tmp;
298  int n;
299 
300
301  //Get double from attribute
302  TObject = PyObject_GetAttrString(O, name);
303  if (!TObject) {
304        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_double could not obtain double %s.\n", name);
305        //printf("name = %s",name);
306    PyErr_SetString(PyExc_RuntimeError, buf);
307
308    return 0.0;
309  } 
310 
311  tmp = PyFloat_AsDouble(TObject);
312 
313  Py_DECREF(TObject);
314 
315  return tmp;
316}
317
318
319
320
321int get_python_integer(PyObject *O, char *name) {
322  PyObject *TObject;
323  #define BUFFER_SIZE 80
324  char buf[BUFFER_SIZE];
325  long tmp;
326  int n;
327 
328
329  //Get double from attribute
330  TObject = PyObject_GetAttrString(O, name);
331  if (!TObject) {
332        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_integer could not obtain double %s.\n", name);
333        //printf("name = %s",name);
334    PyErr_SetString(PyExc_RuntimeError, buf);
335    return 0;
336  } 
337 
338  tmp = PyInt_AsLong(TObject);
339 
340  Py_DECREF(TObject);
341 
342  return tmp;
343}
344
345
346PyObject *get_python_object(PyObject *O, char *name) {
347  PyObject *Oout;
348
349  Oout = PyObject_GetAttrString(O, name);
350  if (!Oout) {
351    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
352    return NULL;
353  }
354
355  return Oout;
356
357}
358
359
360
Note: See TracBrowser for help on using the repository browser.