source: trunk/anuga_work/development/2010-projects/kv/kinematic_viscosity_ext.c @ 8052

Last change on this file since 8052 was 8052, checked in by steve, 13 years ago

Modified C code to take long arrays instead of int (caused error on 64 bit machine)

File size: 7.5 KB
Line 
1#include "Python.h"
2#include "numpy/arrayobject.h"
3#include <math.h>
4#include <stdio.h>
5#include <unistd.h>
6#include "util_ext.h"
7
8//Rough quicksort implementation (for build_operator_matrix)
9int *quicksort(int *array, int n) {
10        int *less, *more, pivot, i, num_less, num_more, *sorted;
11        less = malloc(sizeof(int)*n);
12        more = malloc(sizeof(int)*n);
13        sorted = malloc(sizeof(int)*n);
14        if (n<2) return array;
15        pivot = array[n-1]; //choose the last element as a pivot
16        num_less = 0; num_more = 0;
17        for (i=0; i<n-1; i++) {
18                if (array[i] <= pivot) {
19                        less[num_less] = array[i];
20                        num_less++;
21                } else {
22                        more[num_more] = array[i];
23                        num_more++;
24                }
25        }
26       
27        less = quicksort(less,num_less);
28        more = quicksort(more,num_more);
29       
30        for (i=0; i<n; i++) {
31                if (i<num_less) {
32                        sorted[i] = less[i];
33                } else if (i==num_less) {
34                        sorted[i] = pivot;
35                } else {
36                        sorted[i] = more[i-num_less-1];
37                }
38        }
39        return sorted;
40}
41
42int build_geo_structure(int n, 
43                        int tot_len, 
44                        double *centroids, 
45                        long *neighbours, 
46                        double *edgelengths, 
47                        double *edge_midpoints, 
48                        long *geo_indices, 
49                        double *geo_values) {
50    int i, edge, edge_counted, j;
51        double dist, this_x, this_y, other_x, other_y, edge_length;
52        edge_counted = 0;
53        for (i=0; i<n; i++) {
54                //The centroid coordinates of triangle i
55                this_x = centroids[2*i];
56                this_y = centroids[2*i+1];
57                for (edge=0; edge<3; edge++) {
58               
59                        j = neighbours[3*i+edge];
60                       
61                        //Get the index and the coordinates of the interacting point
62                        if (j == -1) {
63                                geo_indices[3*i+edge] = n + edge_counted;
64                                edge_counted++;
65                               
66                                other_x = edge_midpoints[2*(3*i+edge)];
67                                other_y = edge_midpoints[2*(3*i+edge)+1];
68                        } else {
69                                geo_indices[3*i+edge] = j;
70                               
71                                other_x = centroids[2*j];
72                                other_y = centroids[2*j+1];
73                        }
74                       
75                        //Compute the interaction
76                        edge_length = edgelengths[3*i+edge];
77                        dist = sqrt((this_x-other_x)*(this_x-other_x) + (this_y-other_y)*(this_y-other_y));
78                        geo_values[3*i+edge] = - edge_length / dist;
79                }
80        }
81        return 0;
82}
83
84int build_operator_matrix(int n, 
85                          int tot_len, 
86                          long *geo_indices, 
87                          double *geo_values, 
88                          double *h, 
89                          double *boundary_heights, 
90                          double *data, 
91                          long *colind) {
92        int i, k, edge, j[4], *sorted_j, this_index;
93        double h_j, v[3], v_i; //v[k] = value of the interaction of edge k in a given triangle, v_i = (i,i) entry
94        for (i=0; i<n; i++) {
95                v_i = 0.0;
96                j[3] = i;
97                //Get the values of each interaction, and the column index at which they occur
98                for (edge=0; edge<3; edge++) {
99                        j[edge] = geo_indices[3*i+edge];
100                        if (j[edge]<n) { //interior
101                                h_j = h[j[edge]];
102                        } else { //boundary
103                                h_j = boundary_heights[j[edge]-n];
104                        }
105                        v[edge] = -0.5*(h[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
106                        v_i += 0.5*(h[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
107                }
108                //Organise the set of 4 values/indices into the data and colind arrays
109                sorted_j = quicksort((int *)j,4);
110                for (k=0; k<4; k++) { //loop through the nonzero indices
111                        this_index = sorted_j[k];
112                        if (this_index == i) {
113                                data[4*i+k] = v_i;
114                                colind[4*i+k] = i;
115                        } else if (this_index == j[0]) {
116                                data[4*i+k] = v[0];
117                                colind[4*i+k] = j[0];
118                        } else if (this_index == j[1]) {
119                                data[4*i+k] = v[1];
120                                colind[4*i+k] = j[1];
121                        } else { //this_index == j[2]
122                                data[4*i+k] = v[2];
123                                colind[4*i+k] = j[2];
124                        }
125                }
126        }
127        return 0;
128}
129
130/**
131* Python wrapper methods
132*/
133
134static PyObject *py_build_geo_structure(PyObject *self, PyObject *args) {
135    PyObject *kv_operator, *mesh;
136    int n, tot_len, err;
137    PyArrayObject
138            *centroid_coordinates,
139            *neighbours,
140            *edgelengths,
141            *edge_midpoint_coordinates,
142                        *geo_indices,
143                        *geo_values;
144   
145    //Convert Python arguments to C
146    if (!PyArg_ParseTuple(args, "Oii", &kv_operator, &n, &tot_len)) {
147        PyErr_SetString(PyExc_RuntimeError, "build_geo_structure could not parse input");
148        return NULL;
149    }
150    mesh = PyObject_GetAttrString(kv_operator, "mesh"); //kv_operator.mesh
151    if (!mesh) {
152        PyErr_SetString(PyExc_RuntimeError, "build_geo_structure could not obtain mesh object from kv_operator");
153        return NULL;
154    }
155
156    //Extract parameters
157    centroid_coordinates = get_consecutive_array(mesh,"centroid_coordinates");
158    neighbours = get_consecutive_array(mesh,"neighbours");
159    edgelengths = get_consecutive_array(mesh,"edgelengths");
160    edge_midpoint_coordinates = get_consecutive_array(mesh,"edge_midpoint_coordinates");
161        geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
162        geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
163   
164    //Release
165    Py_DECREF(mesh);
166   
167    err = build_geo_structure(n,tot_len, 
168                             (double *)centroid_coordinates -> data, 
169                             (long *)neighbours -> data, 
170                             (double *)edgelengths->data, 
171                             (double *)edge_midpoint_coordinates -> data, 
172                             (long *)geo_indices -> data, 
173                             (double *)geo_values -> data);
174    if (err != 0) {
175        PyErr_SetString(PyExc_RuntimeError, "Could not build geo structure");
176        return NULL;
177    }
178   
179    //Release the arrays
180    Py_DECREF(centroid_coordinates);
181    Py_DECREF(neighbours);
182    Py_DECREF(edgelengths);
183    Py_DECREF(edge_midpoint_coordinates);
184        Py_DECREF(geo_indices);
185        Py_DECREF(geo_values);
186   
187    return Py_BuildValue("");
188}
189
190static PyObject *py_build_operator_matrix(PyObject *self, PyObject *args) {
191        PyObject *kv_operator;
192        int n, tot_len, err;
193        PyArrayObject
194                *h,
195                *boundary_heights,
196                *geo_indices,
197                *geo_values,
198                *_data,
199                *colind;
200       
201        //Convert Python arguments to C
202    if (!PyArg_ParseTuple(args, "OiiOO", &kv_operator, &n, &tot_len, &h, &boundary_heights)) {
203        PyErr_SetString(PyExc_RuntimeError, "get_stage_height_interactions could not parse input");
204        return NULL;
205    }
206       
207        geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
208        geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
209        _data = get_consecutive_array(kv_operator,"operator_data");
210        colind = get_consecutive_array(kv_operator,"operator_colind");
211       
212        err = build_operator_matrix(n,tot_len, 
213                                                                (long *)geo_indices -> data, 
214                                                                (double *)geo_values -> data, 
215                                                                (double *)h -> data, 
216                                                                (double *)boundary_heights -> data, 
217                                                                (double *)_data -> data, 
218                                                                (long *)colind -> data);
219    if (err != 0) {
220        PyErr_SetString(PyExc_RuntimeError, "Could not get stage height interactions");
221        return NULL;
222    }
223       
224        Py_DECREF(geo_indices);
225        Py_DECREF(geo_values);
226        Py_DECREF(_data);
227        Py_DECREF(colind);
228   
229    return Py_BuildValue("");
230}
231
232static struct PyMethodDef MethodTable[] = {
233        {"build_geo_structure",py_build_geo_structure,METH_VARARGS,"Print out"},
234                {"build_operator_matrix",py_build_operator_matrix,METH_VARARGS,"Print out"},
235        {NULL,NULL,0,NULL} // sentinel
236};
237void initkinematic_viscosity_ext(){
238    (void) Py_InitModule("kinematic_viscosity_ext", MethodTable);
239    import_array();
240}
Note: See TracBrowser for help on using the repository browser.