source: anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c @ 4728

Last change on this file since 4728 was 4728, checked in by ole, 17 years ago

Cosmetics

File size: 18.1 KB
Line 
1// Python - C extension for quantity 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 quantity.py
8//
9//
10// Ole Nielsen, GA 2004
11
12#include "Python.h"
13#include "Numeric/arrayobject.h"
14#include "math.h"
15
16//Shared code snippets
17#include "util_ext.h"
18//#include "quantity_ext.h" //Obsolete
19
20
21
22int _compute_gradients(int N,
23                        double* centroids,
24                        double* centroid_values,
25                        long* number_of_boundaries,
26                        long* surrogate_neighbours,
27                        double* a,
28                        double* b){
29
30        int i, k, k0, k1, k2, index3;
31        double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
32
33
34        for (k=0; k<N; k++) {
35                index3 = 3*k;
36
37                if (number_of_boundaries[k] < 2) {
38                        // Two or three true neighbours
39
40                        // Get indices of neighbours (or self when used as surrogate)
41                        // k0, k1, k2 = surrogate_neighbours[k,:]
42
43                        k0 = surrogate_neighbours[index3 + 0];
44                        k1 = surrogate_neighbours[index3 + 1];
45                        k2 = surrogate_neighbours[index3 + 2];
46
47
48                        if (k0 == k1 || k1 == k2) return -1;
49
50                        // Get data
51                        q0 = centroid_values[k0];
52                        q1 = centroid_values[k1];
53                        q2 = centroid_values[k2];
54
55                        x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
56                        x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
57                        x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
58
59                        // Gradient
60                        _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
61
62                } else if (number_of_boundaries[k] == 2) {
63                        // One true neighbour
64
65                        // Get index of the one neighbour
66                        i=0; k0 = k;
67                        while (i<3 && k0==k) {
68                                k0 = surrogate_neighbours[index3 + i];
69                                i++;
70                        }
71                        if (k0 == k) return -1;
72
73                        k1 = k; //self
74
75                        // Get data
76                        q0 = centroid_values[k0];
77                        q1 = centroid_values[k1];
78
79                        x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
80                        x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
81
82                        // Two point gradient
83                        _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
84
85                }
86                //    else:
87                //        #No true neighbours -
88                //        #Fall back to first order scheme
89        }
90        return 0;
91}
92
93
94
95int _extrapolate(int N,
96                 double* centroids,
97                 double* centroid_values,
98                 double* vertex_coordinates,
99                 double* vertex_values,
100                 double* a,
101                 double* b) {
102
103        int k, k2, k3, k6;
104        double x, y, x0, y0, x1, y1, x2, y2;
105
106        for (k=0; k<N; k++){
107                k6 = 6*k;
108                k3 = 3*k;
109                k2 = 2*k;
110
111                // Centroid coordinates
112                x = centroids[k2]; y = centroids[k2+1];
113
114                // vertex coordinates
115                // x0, y0, x1, y1, x2, y2 = X[k,:]
116                x0 = vertex_coordinates[k6 + 0];
117                y0 = vertex_coordinates[k6 + 1];
118                x1 = vertex_coordinates[k6 + 2];
119                y1 = vertex_coordinates[k6 + 3];
120                x2 = vertex_coordinates[k6 + 4];
121                y2 = vertex_coordinates[k6 + 5];
122
123                // Extrapolate
124                vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
125                vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
126                vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
127
128        }
129        return 0;
130}
131
132
133
134
135int _interpolate(int N,
136                 double* vertex_values,
137                 double* edge_values) {
138
139        int k, k3;
140        double q0, q1, q2;
141
142
143        for (k=0; k<N; k++) {
144                k3 = 3*k;
145
146                q0 = vertex_values[k3 + 0];
147                q1 = vertex_values[k3 + 1];
148                q2 = vertex_values[k3 + 2];
149
150                edge_values[k3 + 0] = 0.5*(q1+q2);
151                edge_values[k3 + 1] = 0.5*(q0+q2);
152                edge_values[k3 + 2] = 0.5*(q0+q1);
153        }
154        return 0;
155}
156
157int _backup_centroid_values(int N,
158                            double* centroid_values,
159                            double* centroid_backup_values) {
160    // Backup centroid values
161
162
163    int k;
164
165    for (k=0; k<N; k++) {
166        centroid_backup_values[k] = centroid_values[k];
167    }
168
169
170    return 0;
171}
172
173
174int _saxpy_centroid_values(int N,
175                           double a,
176                           double b,
177                           double* centroid_values,
178                           double* centroid_backup_values) {
179    // Saxby centroid values
180
181
182    int k;
183
184
185    for (k=0; k<N; k++) {
186        centroid_values[k] = a*centroid_values[k] + b*centroid_backup_values[k];
187    }
188
189
190    return 0;
191}
192
193
194int _update(int N,
195            double timestep,
196            double* centroid_values,
197            double* explicit_update,
198            double* semi_implicit_update) {
199        // Update centroid values based on values stored in
200        // explicit_update and semi_implicit_update as well as given timestep
201
202
203        int k;
204        double denominator, x;
205
206
207        // Divide semi_implicit update by conserved quantity
208        for (k=0; k<N; k++) {
209                x = centroid_values[k];
210                if (x == 0.0) {
211                        semi_implicit_update[k] = 0.0;
212                } else {
213                        semi_implicit_update[k] /= x;
214                }
215        }
216
217
218        // Semi implicit updates
219        for (k=0; k<N; k++) {
220                denominator = 1.0 - timestep*semi_implicit_update[k];
221                if (denominator == 0.0) {
222                        return -1;
223                } else {
224                        //Update conserved_quantities from semi implicit updates
225                        centroid_values[k] /= denominator;
226                }
227        }
228
229
230        // Explicit updates
231        for (k=0; k<N; k++) {
232                centroid_values[k] += timestep*explicit_update[k];
233        }
234
235
236        // MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
237        for (k=0;k<N;k++){
238                semi_implicit_update[k]=0.0;
239        }
240
241        return 0;
242}
243
244
245int _average_vertex_values(int N,
246                           long* vertex_value_indices,
247                           long* number_of_triangles_per_node,
248                           double* vertex_values, 
249                           double* A) {
250  // Average vertex values to obtain one value per node
251
252  int i, index; 
253  int k = 0; //Track triangles touching each node
254  int current_node = 0;
255  double total = 0.0;
256
257  for (i=0; i<N; i++) {
258 
259    // if (current_node == N) {
260    //   printf("Current node exceeding number of nodes (%d)", N);   
261    //   return 1;
262    // }
263   
264    index = vertex_value_indices[i];
265    k += 1;
266           
267    // volume_id = index / 3
268    // vertex_id = index % 3
269    // total += self.vertex_values[volume_id, vertex_id]
270    total += vertex_values[index];
271     
272    // printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);
273    if (number_of_triangles_per_node[current_node] == k) {
274      A[current_node] = total/k;
275               
276      // Move on to next node
277      total = 0.0;
278      k = 0;
279      current_node += 1;
280    }
281  }
282                           
283  return 0;
284}
285
286
287/////////////////////////////////////////////////
288// Gateways to Python
289PyObject *update(PyObject *self, PyObject *args) {
290
291        PyObject *quantity;
292        PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
293
294        double timestep;
295        int N, err;
296
297
298        // Convert Python arguments to C
299        if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep)) {
300          PyErr_SetString(PyExc_RuntimeError, 
301                          "quantity_ext.c: update could not parse input");
302          return NULL;
303        }
304
305        centroid_values = get_consecutive_array(quantity, "centroid_values");
306        explicit_update = get_consecutive_array(quantity, "explicit_update");
307        semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");
308
309        N = centroid_values -> dimensions[0];
310
311        err = _update(N, timestep,
312                      (double*) centroid_values -> data,
313                      (double*) explicit_update -> data,
314                      (double*) semi_implicit_update -> data);
315
316
317        if (err != 0) {
318          PyErr_SetString(PyExc_RuntimeError,
319                          "Zero division in semi implicit update - call Stephen :)");
320          return NULL;
321        }
322
323        // Release and return
324        Py_DECREF(centroid_values);
325        Py_DECREF(explicit_update);
326        Py_DECREF(semi_implicit_update);
327
328        return Py_BuildValue("");
329}
330
331
332PyObject *backup_centroid_values(PyObject *self, PyObject *args) {
333
334        PyObject *quantity;
335        PyArrayObject *centroid_values, *centroid_backup_values;
336
337        int N, err;
338
339
340        // Convert Python arguments to C
341        if (!PyArg_ParseTuple(args, "O", &quantity)) {
342          PyErr_SetString(PyExc_RuntimeError, 
343                          "quantity_ext.c: backup_centroid_values could not parse input");
344          return NULL;
345        }
346
347        centroid_values        = get_consecutive_array(quantity, "centroid_values");
348        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
349
350        N = centroid_values -> dimensions[0];
351
352        err = _backup_centroid_values(N,
353                      (double*) centroid_values -> data,
354                      (double*) centroid_backup_values -> data);
355
356
357        // Release and return
358        Py_DECREF(centroid_values);
359        Py_DECREF(centroid_backup_values);
360
361        return Py_BuildValue("");
362}
363
364PyObject *saxpy_centroid_values(PyObject *self, PyObject *args) {
365
366        PyObject *quantity;
367        PyArrayObject *centroid_values, *centroid_backup_values;
368
369        double a,b;
370        int N, err;
371
372
373        // Convert Python arguments to C
374        if (!PyArg_ParseTuple(args, "Odd", &quantity, &a, &b)) {
375          PyErr_SetString(PyExc_RuntimeError, 
376                          "quantity_ext.c: saxpy_centroid_values could not parse input");
377          return NULL;
378        }
379
380        centroid_values        = get_consecutive_array(quantity, "centroid_values");
381        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
382
383        N = centroid_values -> dimensions[0];
384
385        err = _saxpy_centroid_values(N,a,b,
386                      (double*) centroid_values -> data,
387                      (double*) centroid_backup_values -> data);
388
389
390        // Release and return
391        Py_DECREF(centroid_values);
392        Py_DECREF(centroid_backup_values);
393
394        return Py_BuildValue("");
395}
396
397
398PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
399
400        PyObject *quantity;
401        PyArrayObject *vertex_values, *edge_values;
402
403        int N, err;
404
405        // Convert Python arguments to C
406        if (!PyArg_ParseTuple(args, "O", &quantity)) {
407          PyErr_SetString(PyExc_RuntimeError, 
408                          "quantity_ext.c: interpolate_from_vertices_to_edges could not parse input");
409          return NULL;
410        }
411       
412        vertex_values = get_consecutive_array(quantity, "vertex_values");
413        edge_values = get_consecutive_array(quantity, "edge_values");
414
415        N = vertex_values -> dimensions[0];
416
417        err = _interpolate(N,
418                           (double*) vertex_values -> data,
419                           (double*) edge_values -> data);
420
421        if (err != 0) {
422          PyErr_SetString(PyExc_RuntimeError, 
423                          "Interpolate could not be computed");
424          return NULL;
425        }
426
427        // Release and return
428        Py_DECREF(vertex_values);
429        Py_DECREF(edge_values);
430
431        return Py_BuildValue("");
432}
433
434
435PyObject *average_vertex_values(PyObject *self, PyObject *args) {
436
437        PyArrayObject
438          *vertex_value_indices, 
439          *number_of_triangles_per_node,
440          *vertex_values,
441          *A;
442       
443
444        int N, err;
445
446        // Convert Python arguments to C
447        if (!PyArg_ParseTuple(args, "OOOO",
448                              &vertex_value_indices, 
449                              &number_of_triangles_per_node,
450                              &vertex_values, 
451                              &A)) {
452          PyErr_SetString(PyExc_RuntimeError, 
453                          "quantity_ext.c: average_vertex_values could not parse input");
454          return NULL;
455        }
456       
457        N = vertex_value_indices -> dimensions[0];
458        // printf("Got parameters, N=%d\n", N);
459        err = _average_vertex_values(N,
460                                     (long*) vertex_value_indices -> data,
461                                     (long*) number_of_triangles_per_node -> data,
462                                     (double*) vertex_values -> data, 
463                                     (double*) A -> data);
464
465        //printf("Error %d", err);
466        if (err != 0) {
467          PyErr_SetString(PyExc_RuntimeError, 
468                          "average_vertex_values could not be computed");
469          return NULL;
470        }
471
472        return Py_BuildValue("");
473}
474
475
476
477PyObject *compute_gradients(PyObject *self, PyObject *args) {
478
479        PyObject *quantity, *domain, *R;
480        PyArrayObject
481                *centroids,            //Coordinates at centroids
482                *centroid_values,      //Values at centroids
483                *number_of_boundaries, //Number of boundaries for each triangle
484                *surrogate_neighbours, //True neighbours or - if one missing - self
485                *a, *b;                //Return values
486
487        int dimensions[1], N, err;
488
489        // Convert Python arguments to C
490        if (!PyArg_ParseTuple(args, "O", &quantity)) {
491          PyErr_SetString(PyExc_RuntimeError, 
492                          "quantity_ext.c: compute_gradients could not parse input");   
493          return NULL;
494        }
495
496        domain = PyObject_GetAttrString(quantity, "domain");
497        if (!domain) {
498          PyErr_SetString(PyExc_RuntimeError, 
499                          "compute_gradients could not obtain domain object from quantity");
500          return NULL;
501        }
502
503        // Get pertinent variables
504
505        centroids = get_consecutive_array(domain, "centroid_coordinates");
506        centroid_values = get_consecutive_array(quantity, "centroid_values");
507        surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
508        number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
509
510        N = centroid_values -> dimensions[0];
511
512        // Release
513        Py_DECREF(domain);
514
515        // Allocate space for return vectors a and b (don't DECREF)
516        dimensions[0] = N;
517        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
518        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
519
520
521
522        err = _compute_gradients(N,
523                        (double*) centroids -> data,
524                        (double*) centroid_values -> data,
525                        (long*) number_of_boundaries -> data,
526                        (long*) surrogate_neighbours -> data,
527                        (double*) a -> data,
528                        (double*) b -> data);
529
530        if (err != 0) {
531          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
532          return NULL;
533        }
534
535        // Release
536        Py_DECREF(centroids);
537        Py_DECREF(centroid_values);
538        Py_DECREF(number_of_boundaries);
539        Py_DECREF(surrogate_neighbours);
540
541        // Build result, release and return
542        R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
543        Py_DECREF(a);
544        Py_DECREF(b);
545        return R;
546}
547
548
549
550PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
551
552        PyObject *quantity, *domain;
553        PyArrayObject
554                *centroids,            //Coordinates at centroids
555                *centroid_values,      //Values at centroids
556                *vertex_coordinates,   //Coordinates at vertices
557                *vertex_values,        //Values at vertices
558                *number_of_boundaries, //Number of boundaries for each triangle
559                *surrogate_neighbours, //True neighbours or - if one missing - self
560                *a, *b;                //Gradients
561
562        //int N, err;
563        int dimensions[1], N, err;
564        //double *a, *b;  //Gradients
565
566        // Convert Python arguments to C
567        if (!PyArg_ParseTuple(args, "O", &quantity)) {
568          PyErr_SetString(PyExc_RuntimeError, 
569                          "extrapolate_second_order could not parse input");   
570          return NULL;
571        }
572
573        domain = PyObject_GetAttrString(quantity, "domain");
574        if (!domain) {
575          PyErr_SetString(PyExc_RuntimeError, 
576                          "extrapolate_second_order could not obtain domain object from quantity");     
577          return NULL;
578        }
579
580        // Get pertinent variables
581        centroids = get_consecutive_array(domain, "centroid_coordinates");
582        centroid_values = get_consecutive_array(quantity, "centroid_values");
583        surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
584        number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
585        vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
586        vertex_values = get_consecutive_array(quantity, "vertex_values");
587
588        N = centroid_values -> dimensions[0];
589
590        // Release
591        Py_DECREF(domain);
592
593        //Allocate space for return vectors a and b (don't DECREF)
594        dimensions[0] = N;
595        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
596        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
597
598        //FIXME: Odd that I couldn't use normal arrays
599        //Allocate space for return vectors a and b (don't DECREF)
600        //a = (double*) malloc(N * sizeof(double));
601        //if (!a) return NULL;
602        //b = (double*) malloc(N * sizeof(double));
603        //if (!b) return NULL;
604
605
606        err = _compute_gradients(N,
607                        (double*) centroids -> data,
608                        (double*) centroid_values -> data,
609                        (long*) number_of_boundaries -> data,
610                        (long*) surrogate_neighbours -> data,
611                        (double*) a -> data,
612                        (double*) b -> data);
613
614        if (err != 0) {
615          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
616          return NULL;
617        }
618
619        err = _extrapolate(N,
620                        (double*) centroids -> data,
621                        (double*) centroid_values -> data,
622                        (double*) vertex_coordinates -> data,
623                        (double*) vertex_values -> data,
624                        (double*) a -> data,
625                        (double*) b -> data);
626
627
628        if (err != 0) {
629          PyErr_SetString(PyExc_RuntimeError,
630                          "Internal function _extrapolate failed");
631          return NULL;
632        }
633
634
635
636        // Release
637        Py_DECREF(centroids);
638        Py_DECREF(centroid_values);
639        Py_DECREF(number_of_boundaries);
640        Py_DECREF(surrogate_neighbours);
641        Py_DECREF(vertex_coordinates);
642        Py_DECREF(vertex_values);
643        Py_DECREF(a);
644        Py_DECREF(b);
645
646        return Py_BuildValue("");
647}
648
649
650
651PyObject *limit(PyObject *self, PyObject *args) {
652
653        PyObject *quantity, *domain, *Tmp;
654        PyArrayObject
655                *qv, //Conserved quantities at vertices
656                *qc, //Conserved quantities at centroids
657                *neighbours;
658
659        int k, i, n, N, k3;
660        double beta_w; //Safety factor
661        double *qmin, *qmax, qn;
662
663        // Convert Python arguments to C
664        if (!PyArg_ParseTuple(args, "O", &quantity)) {
665          PyErr_SetString(PyExc_RuntimeError, 
666                          "quantity_ext.c: limit could not parse input");
667          return NULL;
668        }
669
670        domain = PyObject_GetAttrString(quantity, "domain");
671        if (!domain) {
672          PyErr_SetString(PyExc_RuntimeError, 
673                          "quantity_ext.c: limit could not obtain domain object from quantity");                       
674         
675          return NULL;
676        }
677
678        //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
679        neighbours = get_consecutive_array(domain, "neighbours");
680
681        // Get safety factor beta_w
682        Tmp = PyObject_GetAttrString(domain, "beta_w");
683        if (!Tmp) {
684          PyErr_SetString(PyExc_RuntimeError, 
685                          "quantity_ext.c: limit could not obtain beta_w object from domain");                 
686         
687          return NULL;
688        }       
689
690        beta_w = PyFloat_AsDouble(Tmp);
691
692        Py_DECREF(Tmp);
693        Py_DECREF(domain);
694
695
696        qc = get_consecutive_array(quantity, "centroid_values");
697        qv = get_consecutive_array(quantity, "vertex_values");
698
699
700        N = qc -> dimensions[0];
701
702        // Find min and max of this and neighbour's centroid values
703        qmin = malloc(N * sizeof(double));
704        qmax = malloc(N * sizeof(double));
705        for (k=0; k<N; k++) {
706                k3=k*3;
707
708                qmin[k] = ((double*) qc -> data)[k];
709                qmax[k] = qmin[k];
710
711                for (i=0; i<3; i++) {
712                        n = ((long*) neighbours -> data)[k3+i];
713                        if (n >= 0) {
714                                qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
715
716                                qmin[k] = min(qmin[k], qn);
717                                qmax[k] = max(qmax[k], qn);
718                        }
719                        //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]);
720                        //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]);
721                }
722        }
723
724        // Call underlying routine
725        _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
726
727        free(qmin);
728        free(qmax);
729        return Py_BuildValue("");
730}
731
732
733
734// Method table for python module
735static struct PyMethodDef MethodTable[] = {
736        {"limit", limit, METH_VARARGS, "Print out"},
737        {"update", update, METH_VARARGS, "Print out"},
738        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
739        {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"},
740        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
741        {"extrapolate_second_order", extrapolate_second_order,
742                METH_VARARGS, "Print out"},
743        {"interpolate_from_vertices_to_edges",
744                interpolate_from_vertices_to_edges,
745                METH_VARARGS, "Print out"},
746        {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"},           
747        {NULL, NULL, 0, NULL}   // sentinel
748};
749
750// Module initialisation
751void initquantity_ext(void){
752  Py_InitModule("quantity_ext", MethodTable);
753
754  import_array(); // Necessary for handling of NumPY structures
755}
Note: See TracBrowser for help on using the repository browser.