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

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

Working version of comp_flux_ext

File size: 7.5 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  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  double tmp;
284 
285
286  //Get double from attribute
287  TObject = PyObject_GetAttrString(O, name);
288  if (!TObject) {
289    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_double could not obtain double from object");
290    return 0.0;
291  } 
292 
293  tmp = PyFloat_AsDouble(TObject);
294 
295  Py_DECREF(TObject);
296 
297  return tmp;
298}
299
300int get_python_integer(PyObject *O, char *name) {
301  PyObject *TObject;
302  int tmp;
303 
304
305  //Get double from attribute
306  TObject = PyObject_GetAttrString(O, name);
307  if (!TObject) {
308    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_integer could not obtain double from object");
309    return 0;
310  } 
311 
312  tmp = PyFloat_AsDouble(TObject);
313 
314  Py_DECREF(TObject);
315 
316  return tmp;
317}
318
319
320PyObject *get_python_object(PyObject *O, char *name) {
321  PyObject *Oout;
322
323  Oout = PyObject_GetAttrString(O, name);
324  if (!Oout) {
325    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
326    return NULL;
327  }
328
329  return Oout;
330
331}
332
333
334
Note: See TracBrowser for help on using the repository browser.