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

Last change on this file since 8051 was 8051, checked in by steve, 9 years ago

Added in Lindon's code

File size: 6.9 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, int tot_len, double *centroids, int *neighbours, double *edgelengths, double *edge_midpoints, int *geo_indices, double *geo_values) {
43    int i, edge, edge_counted, j;
44        double dist, this_x, this_y, other_x, other_y, edge_length;
45        edge_counted = 0;
46        for (i=0; i<n; i++) {
47                //The centroid coordinates of triangle i
48                this_x = centroids[2*i];
49                this_y = centroids[2*i+1];
50                for (edge=0; edge<3; edge++) {
51               
52                        j = neighbours[3*i+edge];
53                       
54                        //Get the index and the coordinates of the interacting point
55                        if (j == -1) {
56                                geo_indices[3*i+edge] = n + edge_counted;
57                                edge_counted++;
58                               
59                                other_x = edge_midpoints[2*(3*i+edge)];
60                                other_y = edge_midpoints[2*(3*i+edge)+1];
61                        } else {
62                                geo_indices[3*i+edge] = j;
63                               
64                                other_x = centroids[2*j];
65                                other_y = centroids[2*j+1];
66                        }
67                       
68                        //Compute the interaction
69                        edge_length = edgelengths[3*i+edge];
70                        dist = sqrt((this_x-other_x)*(this_x-other_x) + (this_y-other_y)*(this_y-other_y));
71                        geo_values[3*i+edge] = - edge_length / dist;
72                }
73        }
74        return 0;
75}
76
77int build_operator_matrix(int n, int tot_len, int *geo_indices, double *geo_values, double *h, double *boundary_heights, double *data, int *colind) {
78        int i, k, edge, j[4], *sorted_j, this_index;
79        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
80        for (i=0; i<n; i++) {
81                v_i = 0.0;
82                j[3] = i;
83                //Get the values of each interaction, and the column index at which they occur
84                for (edge=0; edge<3; edge++) {
85                        j[edge] = geo_indices[3*i+edge];
86                        if (j[edge]<n) { //interior
87                                h_j = h[j[edge]];
88                        } else { //boundary
89                                h_j = boundary_heights[j[edge]-n];
90                        }
91                        v[edge] = -0.5*(h[i] + h_j)*geo_values[3*i+edge]; //the negative of the individual interaction
92                        v_i += 0.5*(h[i] + h_j)*geo_values[3*i+edge]; //sum the three interactions
93                }
94                //Organise the set of 4 values/indices into the data and colind arrays
95                sorted_j = quicksort((int *)j,4);
96                for (k=0; k<4; k++) { //loop through the nonzero indices
97                        this_index = sorted_j[k];
98                        if (this_index == i) {
99                                data[4*i+k] = v_i;
100                                colind[4*i+k] = i;
101                        } else if (this_index == j[0]) {
102                                data[4*i+k] = v[0];
103                                colind[4*i+k] = j[0];
104                        } else if (this_index == j[1]) {
105                                data[4*i+k] = v[1];
106                                colind[4*i+k] = j[1];
107                        } else { //this_index == j[2]
108                                data[4*i+k] = v[2];
109                                colind[4*i+k] = j[2];
110                        }
111                }
112        }
113        return 0;
114}
115
116/**
117* Python wrapper methods
118*/
119
120static PyObject *py_build_geo_structure(PyObject *self, PyObject *args) {
121    PyObject *kv_operator, *mesh;
122    int n, tot_len, err;
123    PyArrayObject
124            *centroid_coordinates,
125            *neighbours,
126            *edgelengths,
127            *edge_midpoint_coordinates,
128                        *geo_indices,
129                        *geo_values;
130   
131    //Convert Python arguments to C
132    if (!PyArg_ParseTuple(args, "Oii", &kv_operator, &n, &tot_len)) {
133        PyErr_SetString(PyExc_RuntimeError, "build_geo_structure could not parse input");
134        return NULL;
135    }
136    mesh = PyObject_GetAttrString(kv_operator, "mesh"); //kv_operator.mesh
137    if (!mesh) {
138        PyErr_SetString(PyExc_RuntimeError, "build_geo_structure could not obtain mesh object from kv_operator");
139        return NULL;
140    }
141
142    //Extract parameters
143    centroid_coordinates = get_consecutive_array(mesh,"centroid_coordinates");
144    neighbours = get_consecutive_array(mesh,"neighbours");
145    edgelengths = get_consecutive_array(mesh,"edgelengths");
146    edge_midpoint_coordinates = get_consecutive_array(mesh,"edge_midpoint_coordinates");
147        geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
148        geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
149   
150    //Release
151    Py_DECREF(mesh);
152   
153    err = build_geo_structure(n,tot_len, (double *)centroid_coordinates -> data, (int *)neighbours -> data, (double *)edgelengths->data, (double *)edge_midpoint_coordinates -> data, (int *)geo_indices -> data, (double *)geo_values -> data);
154    if (err != 0) {
155        PyErr_SetString(PyExc_RuntimeError, "Could not build geo structure");
156        return NULL;
157    }
158   
159    //Release the arrays
160    Py_DECREF(centroid_coordinates);
161    Py_DECREF(neighbours);
162    Py_DECREF(edgelengths);
163    Py_DECREF(edge_midpoint_coordinates);
164        Py_DECREF(geo_indices);
165        Py_DECREF(geo_values);
166   
167    return Py_BuildValue("");
168}
169
170static PyObject *py_build_operator_matrix(PyObject *self, PyObject *args) {
171        PyObject *kv_operator;
172        int n, tot_len, err;
173        PyArrayObject
174                *h,
175                *boundary_heights,
176                *geo_indices,
177                *geo_values,
178                *_data,
179                *colind;
180       
181        //Convert Python arguments to C
182    if (!PyArg_ParseTuple(args, "OiiOO", &kv_operator, &n, &tot_len, &h, &boundary_heights)) {
183        PyErr_SetString(PyExc_RuntimeError, "get_stage_height_interactions could not parse input");
184        return NULL;
185    }
186       
187        geo_indices = get_consecutive_array(kv_operator,"geo_structure_indices");
188        geo_values = get_consecutive_array(kv_operator,"geo_structure_values");
189        _data = get_consecutive_array(kv_operator,"operator_data");
190        colind = get_consecutive_array(kv_operator,"operator_colind");
191       
192        err = build_operator_matrix(n,tot_len, (int *)geo_indices -> data, (double *)geo_values -> data, (double *)h -> data, (double *)boundary_heights -> data, (double *)_data -> data, (int *)colind -> data);
193    if (err != 0) {
194        PyErr_SetString(PyExc_RuntimeError, "Could not get stage height interactions");
195        return NULL;
196    }
197       
198        Py_DECREF(geo_indices);
199        Py_DECREF(geo_values);
200        Py_DECREF(_data);
201        Py_DECREF(colind);
202   
203    return Py_BuildValue("");
204}
205
206static struct PyMethodDef MethodTable[] = {
207        {"build_geo_structure",py_build_geo_structure,METH_VARARGS,"Print out"},
208                {"build_operator_matrix",py_build_operator_matrix,METH_VARARGS,"Print out"},
209        {NULL,NULL,0,NULL} // sentinel
210};
211void initkinematic_viscosity_ext(){
212    (void) Py_InitModule("kinematic_viscosity_ext", MethodTable);
213    import_array();
214}
Note: See TracBrowser for help on using the repository browser.