source: inundation/utilities/util_ext.h @ 2508

Last change on this file since 2508 was 2508, checked in by ole, 18 years ago

First step towards moving util_ext.h out from pyvolution as per ticket:31

File size: 5.3 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(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_numeric_array(PyArrayObject *x) { 
169  int i, j;
170  for (i=0; i<x->dimensions[0]; i++) { 
171    for (j=0; j<x->dimensions[1]; j++) {
172      printf("%f ", *(double*) (x->data + i*x->strides[0] + j*x->strides[1]));
173    }
174    printf("\n"); 
175  }
176  printf("\n");   
177}
178
179void print_numeric_vector(PyArrayObject *x) { 
180  int i;
181  for (i=0; i<x->dimensions[0]; i++) {
182    printf("%f ", *(double*) (x->data + i*x->strides[0])); 
183  }
184  printf("\n"); 
185}
186
187PyArrayObject *get_consecutive_array(PyObject *O, char *name) {
188  PyArrayObject *A, *B;
189 
190
191  //Get array object from attribute
192 
193  /*
194  //FIXME: THE TEST DOESN't WORK
195  printf("Err = %d\n", PyObject_HasAttrString(O, name));
196  if (PyObject_HasAttrString(O, name) == 1) {
197    B = (PyArrayObject*) PyObject_GetAttrString(O, name);
198    if (!B) return NULL;
199  } else {
200    return NULL;
201    }
202  */
203   
204  B = (PyArrayObject*) PyObject_GetAttrString(O, name);
205  if (!B) return NULL;     
206 
207  //Convert to consecutive array
208  A = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*) B, 
209                                                    B -> descr -> type, 0, 0);
210 
211  Py_DECREF(B); //FIXME: Is this really needed??
212 
213  if (!A) return NULL;
214  return A;
215}
216
Note: See TracBrowser for help on using the repository browser.