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

Last change on this file since 4883 was 4883, checked in by steve, 16 years ago

Starting to develop edge limiting and vertex limiting in
quantity (and later in shallow_water_domain

File size: 20.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* edge_values,
101                 double* a,
102                 double* b) {
103
104        int k, k2, k3, k6;
105        double x, y, x0, y0, x1, y1, x2, y2;
106
107        for (k=0; k<N; k++){
108                k6 = 6*k;
109                k3 = 3*k;
110                k2 = 2*k;
111
112                // Centroid coordinates
113                x = centroids[k2]; y = centroids[k2+1];
114
115                // vertex coordinates
116                // x0, y0, x1, y1, x2, y2 = X[k,:]
117                x0 = vertex_coordinates[k6 + 0];
118                y0 = vertex_coordinates[k6 + 1];
119                x1 = vertex_coordinates[k6 + 2];
120                y1 = vertex_coordinates[k6 + 3];
121                x2 = vertex_coordinates[k6 + 4];
122                y2 = vertex_coordinates[k6 + 5];
123
124                // Extrapolate
125                vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
126                vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
127                vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
128
129
130                edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
131                edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
132                edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
133
134        }
135        return 0;
136}
137
138
139
140
141int _interpolate(int N,
142                 double* vertex_values,
143                 double* edge_values) {
144
145        int k, k3;
146        double q0, q1, q2;
147
148
149        for (k=0; k<N; k++) {
150                k3 = 3*k;
151
152                q0 = vertex_values[k3 + 0];
153                q1 = vertex_values[k3 + 1];
154                q2 = vertex_values[k3 + 2];
155
156                edge_values[k3 + 0] = 0.5*(q1+q2);
157                edge_values[k3 + 1] = 0.5*(q0+q2);
158                edge_values[k3 + 2] = 0.5*(q0+q1);
159        }
160        return 0;
161}
162
163int _backup_centroid_values(int N,
164                            double* centroid_values,
165                            double* centroid_backup_values) {
166    // Backup centroid values
167
168
169    int k;
170
171    for (k=0; k<N; k++) {
172        centroid_backup_values[k] = centroid_values[k];
173    }
174
175
176    return 0;
177}
178
179
180int _saxpy_centroid_values(int N,
181                           double a,
182                           double b,
183                           double* centroid_values,
184                           double* centroid_backup_values) {
185    // Saxby centroid values
186
187
188    int k;
189
190
191    for (k=0; k<N; k++) {
192        centroid_values[k] = a*centroid_values[k] + b*centroid_backup_values[k];
193    }
194
195
196    return 0;
197}
198
199
200int _update(int N,
201            double timestep,
202            double* centroid_values,
203            double* explicit_update,
204            double* semi_implicit_update) {
205        // Update centroid values based on values stored in
206        // explicit_update and semi_implicit_update as well as given timestep
207
208
209        int k;
210        double denominator, x;
211
212
213        // Divide semi_implicit update by conserved quantity
214        for (k=0; k<N; k++) {
215                x = centroid_values[k];
216                if (x == 0.0) {
217                        semi_implicit_update[k] = 0.0;
218                } else {
219                        semi_implicit_update[k] /= x;
220                }
221        }
222
223
224        // Semi implicit updates
225        for (k=0; k<N; k++) {
226                denominator = 1.0 - timestep*semi_implicit_update[k];
227                if (denominator == 0.0) {
228                        return -1;
229                } else {
230                        //Update conserved_quantities from semi implicit updates
231                        centroid_values[k] /= denominator;
232                }
233        }
234
235
236        // Explicit updates
237        for (k=0; k<N; k++) {
238                centroid_values[k] += timestep*explicit_update[k];
239        }
240
241
242        // MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
243        for (k=0;k<N;k++){
244                semi_implicit_update[k]=0.0;
245        }
246
247        return 0;
248}
249
250
251int _average_vertex_values(int N,
252                           long* vertex_value_indices,
253                           long* number_of_triangles_per_node,
254                           double* vertex_values, 
255                           double* A) {
256  // Average vertex values to obtain one value per node
257
258  int i, index; 
259  int k = 0; //Track triangles touching each node
260  int current_node = 0;
261  double total = 0.0;
262
263  for (i=0; i<N; i++) {
264 
265    // if (current_node == N) {
266    //   printf("Current node exceeding number of nodes (%d)", N);   
267    //   return 1;
268    // }
269   
270    index = vertex_value_indices[i];
271    k += 1;
272           
273    // volume_id = index / 3
274    // vertex_id = index % 3
275    // total += self.vertex_values[volume_id, vertex_id]
276    total += vertex_values[index];
277     
278    // printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);
279    if (number_of_triangles_per_node[current_node] == k) {
280      A[current_node] = total/k;
281               
282      // Move on to next node
283      total = 0.0;
284      k = 0;
285      current_node += 1;
286    }
287  }
288                           
289  return 0;
290}
291
292
293/////////////////////////////////////////////////
294// Gateways to Python
295PyObject *update(PyObject *self, PyObject *args) {
296  // FIXME (Ole): It would be great to turn this text into a Python DOC string
297
298  /*"""Update centroid values based on values stored in
299    explicit_update and semi_implicit_update as well as given timestep
300
301    Function implementing forcing terms must take on argument
302    which is the domain and they must update either explicit
303    or implicit updates, e,g,:
304
305    def gravity(domain):
306        ....
307        domain.quantities['xmomentum'].explicit_update = ...
308        domain.quantities['ymomentum'].explicit_update = ...
309
310
311
312    Explicit terms must have the form
313
314        G(q, t)
315
316    and explicit scheme is
317
318       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
319
320
321    Semi implicit forcing terms are assumed to have the form
322
323       G(q, t) = H(q, t) q
324
325    and the semi implicit scheme will then be
326
327      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
328
329  */
330
331        PyObject *quantity;
332        PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
333
334        double timestep;
335        int N, err;
336
337
338        // Convert Python arguments to C
339        if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep)) {
340          PyErr_SetString(PyExc_RuntimeError, 
341                          "quantity_ext.c: update could not parse input");
342          return NULL;
343        }
344
345        centroid_values = get_consecutive_array(quantity, "centroid_values");
346        explicit_update = get_consecutive_array(quantity, "explicit_update");
347        semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");
348
349        N = centroid_values -> dimensions[0];
350
351        err = _update(N, timestep,
352                      (double*) centroid_values -> data,
353                      (double*) explicit_update -> data,
354                      (double*) semi_implicit_update -> data);
355
356
357        if (err != 0) {
358          PyErr_SetString(PyExc_RuntimeError,
359                          "Zero division in semi implicit update - call Stephen :)");
360          return NULL;
361        }
362
363        // Release and return
364        Py_DECREF(centroid_values);
365        Py_DECREF(explicit_update);
366        Py_DECREF(semi_implicit_update);
367
368        return Py_BuildValue("");
369}
370
371
372PyObject *backup_centroid_values(PyObject *self, PyObject *args) {
373
374        PyObject *quantity;
375        PyArrayObject *centroid_values, *centroid_backup_values;
376
377        int N, err;
378
379
380        // Convert Python arguments to C
381        if (!PyArg_ParseTuple(args, "O", &quantity)) {
382          PyErr_SetString(PyExc_RuntimeError, 
383                          "quantity_ext.c: backup_centroid_values could not parse input");
384          return NULL;
385        }
386
387        centroid_values        = get_consecutive_array(quantity, "centroid_values");
388        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
389
390        N = centroid_values -> dimensions[0];
391
392        err = _backup_centroid_values(N,
393                      (double*) centroid_values -> data,
394                      (double*) centroid_backup_values -> data);
395
396
397        // Release and return
398        Py_DECREF(centroid_values);
399        Py_DECREF(centroid_backup_values);
400
401        return Py_BuildValue("");
402}
403
404PyObject *saxpy_centroid_values(PyObject *self, PyObject *args) {
405
406        PyObject *quantity;
407        PyArrayObject *centroid_values, *centroid_backup_values;
408
409        double a,b;
410        int N, err;
411
412
413        // Convert Python arguments to C
414        if (!PyArg_ParseTuple(args, "Odd", &quantity, &a, &b)) {
415          PyErr_SetString(PyExc_RuntimeError, 
416                          "quantity_ext.c: saxpy_centroid_values could not parse input");
417          return NULL;
418        }
419
420        centroid_values        = get_consecutive_array(quantity, "centroid_values");
421        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
422
423        N = centroid_values -> dimensions[0];
424
425        err = _saxpy_centroid_values(N,a,b,
426                      (double*) centroid_values -> data,
427                      (double*) centroid_backup_values -> data);
428
429
430        // Release and return
431        Py_DECREF(centroid_values);
432        Py_DECREF(centroid_backup_values);
433
434        return Py_BuildValue("");
435}
436
437
438PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
439        //
440        //Compute edge values from vertex values using linear interpolation
441        //
442       
443        PyObject *quantity;
444        PyArrayObject *vertex_values, *edge_values;
445
446        int N, err;
447
448        // Convert Python arguments to C
449        if (!PyArg_ParseTuple(args, "O", &quantity)) {
450          PyErr_SetString(PyExc_RuntimeError, 
451                          "quantity_ext.c: interpolate_from_vertices_to_edges could not parse input");
452          return NULL;
453        }
454       
455        vertex_values = get_consecutive_array(quantity, "vertex_values");
456        edge_values = get_consecutive_array(quantity, "edge_values");
457
458        N = vertex_values -> dimensions[0];
459
460        err = _interpolate(N,
461                           (double*) vertex_values -> data,
462                           (double*) edge_values -> data);
463
464        if (err != 0) {
465          PyErr_SetString(PyExc_RuntimeError, 
466                          "Interpolate could not be computed");
467          return NULL;
468        }
469
470        // Release and return
471        Py_DECREF(vertex_values);
472        Py_DECREF(edge_values);
473
474        return Py_BuildValue("");
475}
476
477
478PyObject *average_vertex_values(PyObject *self, PyObject *args) {
479
480        PyArrayObject
481          *vertex_value_indices, 
482          *number_of_triangles_per_node,
483          *vertex_values,
484          *A;
485       
486
487        int N, err;
488
489        // Convert Python arguments to C
490        if (!PyArg_ParseTuple(args, "OOOO",
491                              &vertex_value_indices, 
492                              &number_of_triangles_per_node,
493                              &vertex_values, 
494                              &A)) {
495          PyErr_SetString(PyExc_RuntimeError, 
496                          "quantity_ext.c: average_vertex_values could not parse input");
497          return NULL;
498        }
499       
500        N = vertex_value_indices -> dimensions[0];
501        // printf("Got parameters, N=%d\n", N);
502        err = _average_vertex_values(N,
503                                     (long*) vertex_value_indices -> data,
504                                     (long*) number_of_triangles_per_node -> data,
505                                     (double*) vertex_values -> data, 
506                                     (double*) A -> data);
507
508        //printf("Error %d", err);
509        if (err != 0) {
510          PyErr_SetString(PyExc_RuntimeError, 
511                          "average_vertex_values could not be computed");
512          return NULL;
513        }
514
515        return Py_BuildValue("");
516}
517
518
519
520PyObject *compute_gradients(PyObject *self, PyObject *args) {
521  //"""Compute gradients of triangle surfaces defined by centroids of
522  //neighbouring volumes.
523  //If one edge is on the boundary, use own centroid as neighbour centroid.
524  //If two or more are on the boundary, fall back to first order scheme.
525  //"""
526
527
528        PyObject *quantity, *domain, *R;
529        PyArrayObject
530                *centroids,            //Coordinates at centroids
531                *centroid_values,      //Values at centroids
532                *number_of_boundaries, //Number of boundaries for each triangle
533                *surrogate_neighbours, //True neighbours or - if one missing - self
534                *a, *b;                //Return values
535
536        int dimensions[1], N, err;
537
538        // Convert Python arguments to C
539        if (!PyArg_ParseTuple(args, "O", &quantity)) {
540          PyErr_SetString(PyExc_RuntimeError, 
541                          "quantity_ext.c: compute_gradients could not parse input");   
542          return NULL;
543        }
544
545        domain = PyObject_GetAttrString(quantity, "domain");
546        if (!domain) {
547          PyErr_SetString(PyExc_RuntimeError, 
548                          "compute_gradients could not obtain domain object from quantity");
549          return NULL;
550        }
551
552        // Get pertinent variables
553
554        centroids = get_consecutive_array(domain, "centroid_coordinates");
555        centroid_values = get_consecutive_array(quantity, "centroid_values");
556        surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
557        number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
558
559        N = centroid_values -> dimensions[0];
560
561        // Release
562        Py_DECREF(domain);
563
564        // Allocate space for return vectors a and b (don't DECREF)
565        dimensions[0] = N;
566        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
567        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
568
569
570
571        err = _compute_gradients(N,
572                        (double*) centroids -> data,
573                        (double*) centroid_values -> data,
574                        (long*) number_of_boundaries -> data,
575                        (long*) surrogate_neighbours -> data,
576                        (double*) a -> data,
577                        (double*) b -> data);
578
579        if (err != 0) {
580          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
581          return NULL;
582        }
583
584        // Release
585        Py_DECREF(centroids);
586        Py_DECREF(centroid_values);
587        Py_DECREF(number_of_boundaries);
588        Py_DECREF(surrogate_neighbours);
589
590        // Build result, release and return
591        R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
592        Py_DECREF(a);
593        Py_DECREF(b);
594        return R;
595}
596
597
598
599PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
600
601        PyObject *quantity, *domain;
602        PyArrayObject
603            *centroids,            //Coordinates at centroids
604            *centroid_values,      //Values at centroids
605            *vertex_coordinates,   //Coordinates at vertices
606            *vertex_values,        //Values at vertices
607            *edge_values,          //Values at edges
608            *number_of_boundaries, //Number of boundaries for each triangle
609            *surrogate_neighbours, //True neighbours or - if one missing - self
610            *a, *b;                //Gradients
611
612        //int N, err;
613        int dimensions[1], N, err;
614        //double *a, *b;  //Gradients
615
616        // Convert Python arguments to C
617        if (!PyArg_ParseTuple(args, "O", &quantity)) {
618          PyErr_SetString(PyExc_RuntimeError, 
619                          "extrapolate_second_order could not parse input");   
620          return NULL;
621        }
622
623        domain = PyObject_GetAttrString(quantity, "domain");
624        if (!domain) {
625          PyErr_SetString(PyExc_RuntimeError, 
626                          "extrapolate_second_order could not obtain domain object from quantity");     
627          return NULL;
628        }
629
630        // Get pertinent variables
631        centroids            = get_consecutive_array(domain,   "centroid_coordinates");
632        centroid_values      = get_consecutive_array(quantity, "centroid_values");
633        surrogate_neighbours = get_consecutive_array(domain,   "surrogate_neighbours");
634        number_of_boundaries = get_consecutive_array(domain,   "number_of_boundaries");
635        vertex_coordinates   = get_consecutive_array(domain,   "vertex_coordinates");
636        vertex_values        = get_consecutive_array(quantity, "vertex_values");
637        edge_values          = get_consecutive_array(quantity, "edge_values");
638
639        N = centroid_values -> dimensions[0];
640
641        // Release
642        Py_DECREF(domain);
643
644        //Allocate space for return vectors a and b (don't DECREF)
645        dimensions[0] = N;
646        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
647        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
648
649        //FIXME: Odd that I couldn't use normal arrays
650        //Allocate space for return vectors a and b (don't DECREF)
651        //a = (double*) malloc(N * sizeof(double));
652        //if (!a) return NULL;
653        //b = (double*) malloc(N * sizeof(double));
654        //if (!b) return NULL;
655
656
657        err = _compute_gradients(N,
658                        (double*) centroids -> data,
659                        (double*) centroid_values -> data,
660                        (long*) number_of_boundaries -> data,
661                        (long*) surrogate_neighbours -> data,
662                        (double*) a -> data,
663                        (double*) b -> data);
664
665        if (err != 0) {
666          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
667          return NULL;
668        }
669
670        err = _extrapolate(N,
671                        (double*) centroids -> data,
672                        (double*) centroid_values -> data,
673                        (double*) vertex_coordinates -> data,
674                        (double*) vertex_values -> data,
675                        (double*) edge_values -> data,
676                        (double*) a -> data,
677                        (double*) b -> data);
678
679
680        if (err != 0) {
681          PyErr_SetString(PyExc_RuntimeError,
682                          "Internal function _extrapolate failed");
683          return NULL;
684        }
685
686
687
688        // Release
689        Py_DECREF(centroids);
690        Py_DECREF(centroid_values);
691        Py_DECREF(number_of_boundaries);
692        Py_DECREF(surrogate_neighbours);
693        Py_DECREF(vertex_coordinates);
694        Py_DECREF(vertex_values);
695        Py_DECREF(edge_values);
696        Py_DECREF(a);
697        Py_DECREF(b);
698
699        return Py_BuildValue("");
700}
701
702
703
704PyObject *limit_old(PyObject *self, PyObject *args) {
705  //Limit slopes for each volume to eliminate artificial variance
706  //introduced by e.g. second order extrapolator
707
708  //This is an unsophisticated limiter as it does not take into
709  //account dependencies among quantities.
710
711  //precondition:
712  //  vertex values are estimated from gradient
713  //postcondition:
714  //  vertex values are updated
715  //
716
717        PyObject *quantity, *domain, *Tmp;
718        PyArrayObject
719            *qv, //Conserved quantities at vertices
720            *qc, //Conserved quantities at centroids
721            *neighbours;
722
723        int k, i, n, N, k3;
724        double beta_w; //Safety factor
725        double *qmin, *qmax, qn;
726
727        // Convert Python arguments to C
728        if (!PyArg_ParseTuple(args, "O", &quantity)) {
729          PyErr_SetString(PyExc_RuntimeError, 
730                          "quantity_ext.c: limit_old could not parse input");
731          return NULL;
732        }
733
734        domain = PyObject_GetAttrString(quantity, "domain");
735        if (!domain) {
736          PyErr_SetString(PyExc_RuntimeError, 
737                          "quantity_ext.c: limit_old could not obtain domain object from quantity");                   
738         
739          return NULL;
740        }
741
742        //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
743        neighbours = get_consecutive_array(domain, "neighbours");
744
745        // Get safety factor beta_w
746        Tmp = PyObject_GetAttrString(domain, "beta_w");
747        if (!Tmp) {
748          PyErr_SetString(PyExc_RuntimeError, 
749                          "quantity_ext.c: limit_old could not obtain beta_w object from domain");                     
750         
751          return NULL;
752        }       
753
754        beta_w = PyFloat_AsDouble(Tmp);
755
756        Py_DECREF(Tmp);
757        Py_DECREF(domain);
758
759
760        qc = get_consecutive_array(quantity, "centroid_values");
761        qv = get_consecutive_array(quantity, "vertex_values");
762
763
764        N = qc -> dimensions[0];
765
766        // Find min and max of this and neighbour's centroid values
767        qmin = malloc(N * sizeof(double));
768        qmax = malloc(N * sizeof(double));
769        for (k=0; k<N; k++) {
770                k3=k*3;
771
772                qmin[k] = ((double*) qc -> data)[k];
773                qmax[k] = qmin[k];
774
775                for (i=0; i<3; i++) {
776                        n = ((long*) neighbours -> data)[k3+i];
777                        if (n >= 0) {
778                                qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
779
780                                qmin[k] = min(qmin[k], qn);
781                                qmax[k] = max(qmax[k], qn);
782                        }
783                        //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]);
784                        //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]);
785                }
786        }
787
788        // Call underlying routine
789        _limit_old(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
790
791        free(qmin);
792        free(qmax);
793        return Py_BuildValue("");
794}
795
796
797
798// Method table for python module
799static struct PyMethodDef MethodTable[] = {
800        {"limit_old", limit_old, METH_VARARGS, "Print out"},
801        {"update", update, METH_VARARGS, "Print out"},
802        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
803        {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"},
804        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
805        {"extrapolate_second_order", extrapolate_second_order,
806                METH_VARARGS, "Print out"},
807        {"interpolate_from_vertices_to_edges",
808                interpolate_from_vertices_to_edges,
809                METH_VARARGS, "Print out"},
810        {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"},           
811        {NULL, NULL, 0, NULL}   // sentinel
812};
813
814// Module initialisation
815void initquantity_ext(void){
816  Py_InitModule("quantity_ext", MethodTable);
817
818  import_array(); // Necessary for handling of NumPY structures
819}
Note: See TracBrowser for help on using the repository browser.