source: trunk/anuga_work/development/2010-projects/anuga_1d/utilities/util_ext.h @ 7884

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

Moving 2010 project

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