source: branches/numpy/anuga/abstract_2d_finite_volumes/quantity_ext.c @ 6883

Last change on this file since 6883 was 6410, checked in by rwilson, 16 years ago

numpy changes.
File size: 56.0 KB
RevLine 
[5897]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"
[6304]13#include "numpy/arrayobject.h"
[5897]14#include "math.h"
15
16//Shared code snippets
17#include "util_ext.h"
18
19
20//-------------------------------------------
21// Low level routines (called from wrappers)
22//------------------------------------------
23
24int _compute_gradients(int N,
25                        double* centroids,
26                        double* centroid_values,
27                        long* number_of_boundaries,
28                        long* surrogate_neighbours,
29                        double* a,
30                        double* b){
31
32        int i, k, k0, k1, k2, index3;
33        double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
34
35
36        for (k=0; k<N; k++) {
37                index3 = 3*k;
38
39                if (number_of_boundaries[k] < 2) {
40                        // Two or three true neighbours
41
42                        // Get indices of neighbours (or self when used as surrogate)
43                        // k0, k1, k2 = surrogate_neighbours[k,:]
44
45                        k0 = surrogate_neighbours[index3 + 0];
46                        k1 = surrogate_neighbours[index3 + 1];
47                        k2 = surrogate_neighbours[index3 + 2];
48
49
50                        if (k0 == k1 || k1 == k2) return -1;
51
52                        // Get data
53                        q0 = centroid_values[k0];
54                        q1 = centroid_values[k1];
55                        q2 = centroid_values[k2];
56
57                        x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
58                        x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
59                        x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
60
61                        // Gradient
62                        _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
63
64                } else if (number_of_boundaries[k] == 2) {
65                        // One true neighbour
66
67                        // Get index of the one neighbour
68                        i=0; k0 = k;
69                        while (i<3 && k0==k) {
70                                k0 = surrogate_neighbours[index3 + i];
71                                i++;
72                        }
73                        if (k0 == k) return -1;
74
75                        k1 = k; //self
76
77                        // Get data
78                        q0 = centroid_values[k0];
79                        q1 = centroid_values[k1];
80
81                        x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
82                        x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
83
84                        // Two point gradient
85                        _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
86
87                }
88                //    else:
89                //        #No true neighbours -
90                //        #Fall back to first order scheme
91        }
92        return 0;
93}
94
95
96int _extrapolate_from_gradient(int N,
97                 double* centroids,
98                 double* centroid_values,
99                 double* vertex_coordinates,
100                 double* vertex_values,
101                 double* edge_values,
102                 double* a,
103                 double* b) {
104
105        int k, k2, k3, k6;
106        double x, y, x0, y0, x1, y1, x2, y2;
107
108        for (k=0; k<N; k++){
109                k6 = 6*k;
110                k3 = 3*k;
111                k2 = 2*k;
112
113                // Centroid coordinates
114                x = centroids[k2]; y = centroids[k2+1];
115
116                // vertex coordinates
117                // x0, y0, x1, y1, x2, y2 = X[k,:]
118                x0 = vertex_coordinates[k6 + 0];
119                y0 = vertex_coordinates[k6 + 1];
120                x1 = vertex_coordinates[k6 + 2];
121                y1 = vertex_coordinates[k6 + 3];
122                x2 = vertex_coordinates[k6 + 4];
123                y2 = vertex_coordinates[k6 + 5];
124
125                // Extrapolate to Vertices
126                vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
127                vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
128                vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
129
130                // Extrapolate to Edges (midpoints)
131                edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
132                edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
133                edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
134
135        }
136        return 0;
137}
138
139
140int _extrapolate_and_limit_from_gradient(int N,double beta,
141                                         double* centroids,
142                                         long*   neighbours,
143                                         double* centroid_values,
144                                         double* vertex_coordinates,
145                                         double* vertex_values,
146                                         double* edge_values,
147                                         double* phi,
148                                         double* x_gradient,
149                                         double* y_gradient) {
150
151        int i, k, k2, k3, k6;
152        double x, y, x0, y0, x1, y1, x2, y2;
153        long n;
154        double qmin, qmax, qc;
155        double qn[3];
156        double dq, dqa[3], r;
157
158        for (k=0; k<N; k++){
159                k6 = 6*k;
160                k3 = 3*k;
161                k2 = 2*k;
162
163                // Centroid coordinates
164                x = centroids[k2+0]; 
165                y = centroids[k2+1];
166
167                // vertex coordinates
168                // x0, y0, x1, y1, x2, y2 = X[k,:]
169                x0 = vertex_coordinates[k6 + 0];
170                y0 = vertex_coordinates[k6 + 1];
171                x1 = vertex_coordinates[k6 + 2];
172                y1 = vertex_coordinates[k6 + 3];
173                x2 = vertex_coordinates[k6 + 4];
174                y2 = vertex_coordinates[k6 + 5];
175
176                // Extrapolate to Vertices
177                vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
178                vertex_values[k3+1] = centroid_values[k] + x_gradient[k]*(x1-x) + y_gradient[k]*(y1-y);
179                vertex_values[k3+2] = centroid_values[k] + x_gradient[k]*(x2-x) + y_gradient[k]*(y2-y);
180
181                // Extrapolate to Edges (midpoints)
182                edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
183                edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
184                edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
185        }
186
187
188
189        for (k=0; k<N; k++){
190                k6 = 6*k;
191                k3 = 3*k;
192                k2 = 2*k;
193
194
195                qc = centroid_values[k];
196               
197                qmin = qc;
198                qmax = qc;
199               
200                for (i=0; i<3; i++) {
201                    n = neighbours[k3+i];
202                    if (n < 0) {
203                        qn[i] = qc;
204                    } else {
205                        qn[i] = centroid_values[n];
206                    }
207
208                    qmin = min(qmin, qn[i]);
209                    qmax = max(qmax, qn[i]);
210                }
211
212                //qtmin = min(min(min(qn[0],qn[1]),qn[2]),qc);
213                //qtmax = max(max(max(qn[0],qn[1]),qn[2]),qc);
214
215/*              for (i=0; i<3; i++) { */
216/*                  n = neighbours[k3+i]; */
217/*                  if (n < 0) { */
218/*                      qn[i] = qc; */
219/*                      qmin[i] = qtmin; */
220/*                      qmax[i] = qtmax; */
221/*                  }  */
222/*              } */
223                       
224                phi[k] = 1.0;
225
226                for (i=0; i<3; i++) {
227                    dq = edge_values[k3+i] - qc;      //Delta between edge and centroid values
228                    dqa[i] = dq;                      //Save dq for use in updating vertex values
229
230                    r = 1.0;
231     
232                    if (dq > 0.0) r = (qmax - qc)/dq;
233                    if (dq < 0.0) r = (qmin - qc)/dq;
234                   
235                    phi[k] = min( min(r*beta, 1.0), phi[k]);
236                   
237                }
238
239
240
241                //Update gradient, edge and vertex values using phi limiter
242                x_gradient[k] = x_gradient[k]*phi[k];
243                y_gradient[k] = y_gradient[k]*phi[k];
244
245                edge_values[k3+0] = qc + phi[k]*dqa[0];
246                edge_values[k3+1] = qc + phi[k]*dqa[1];
247                edge_values[k3+2] = qc + phi[k]*dqa[2];
248
249
250                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
251                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
252                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
253               
254
255        }
256
257        return 0;
258
259}
260
261
262
263
264int _limit_vertices_by_all_neighbours(int N, double beta,
265                     double* centroid_values,
266                     double* vertex_values,
267                     double* edge_values,
268                     long*   neighbours,
269                     double* x_gradient,
270                     double* y_gradient) {
271
272
273        int i, k, k2, k3, k6;
274        long n;
275        double qmin, qmax, qn, qc;
276        double dq, dqa[3], phi, r;
277
278        for (k=0; k<N; k++){
279                k6 = 6*k;
280                k3 = 3*k;
281                k2 = 2*k;
282
283                qc = centroid_values[k];
284                qmin = qc;
285                qmax = qc;
286
287                for (i=0; i<3; i++) {
288                    n = neighbours[k3+i];
289                    if (n >= 0) {
290                        qn = centroid_values[n]; //Neighbour's centroid value
291
292                        qmin = min(qmin, qn);
293                        qmax = max(qmax, qn);
294                    }
295                }
296
297                phi = 1.0;
298                for (i=0; i<3; i++) {   
299                    r = 1.0;
300     
301                    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
302                    dqa[i] = dq;                      //Save dq for use in updating vertex values
303     
304                    if (dq > 0.0) r = (qmax - qc)/dq;
305                    if (dq < 0.0) r = (qmin - qc)/dq;     
306 
307                   
308                    phi = min( min(r*beta, 1.0), phi);   
309                }
310
311                //Update gradient, vertex and edge values using phi limiter
312                x_gradient[k] = x_gradient[k]*phi;
313                y_gradient[k] = y_gradient[k]*phi;
314   
315                vertex_values[k3+0] = qc + phi*dqa[0];
316                vertex_values[k3+1] = qc + phi*dqa[1];
317                vertex_values[k3+2] = qc + phi*dqa[2];
318               
319                edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
320                edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
321                edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
322
323        }
324
325        return 0;
326}
327
328
329
330
331int _limit_edges_by_all_neighbours(int N, double beta,
332                                   double* centroid_values,
333                                   double* vertex_values,
334                                   double* edge_values,
335                                   long*   neighbours,
336                                   double* x_gradient,
337                                   double* y_gradient) {
338
339        int i, k, k2, k3, k6;
340        long n;
341        double qmin, qmax, qn, qc;
342        double dq, dqa[3], phi, r;
343
344        for (k=0; k<N; k++){
345                k6 = 6*k;
346                k3 = 3*k;
347                k2 = 2*k;
348
349                qc = centroid_values[k];
350                qmin = qc;
351                qmax = qc;
352
353                for (i=0; i<3; i++) {
354                    n = neighbours[k3+i];
355                    if (n >= 0) {
356                        qn = centroid_values[n]; //Neighbour's centroid value
357
358                        qmin = min(qmin, qn);
359                        qmax = max(qmax, qn);
360                    }
361                }
362
363                phi = 1.0;
364                for (i=0; i<3; i++) {   
365                    r = 1.0;
366     
367                    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
368                    dqa[i] = dq;                      //Save dq for use in updating vertex values
369     
370                    if (dq > 0.0) r = (qmax - qc)/dq;
371                    if (dq < 0.0) r = (qmin - qc)/dq;     
372 
373                   
374                    phi = min( min(r*beta, 1.0), phi);   
375                }
376   
377                //Update gradient, vertex and edge values using phi limiter
378                x_gradient[k] = x_gradient[k]*phi;
379                y_gradient[k] = y_gradient[k]*phi;
380
381                edge_values[k3+0] = qc + phi*dqa[0];
382                edge_values[k3+1] = qc + phi*dqa[1];
383                edge_values[k3+2] = qc + phi*dqa[2];
384               
385                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
386                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
387                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
388
389        }
390
391        return 0;
392}
393
394
395int _limit_edges_by_neighbour(int N, double beta,
396                     double* centroid_values,
397                     double* vertex_values,
398                     double* edge_values,
399                     long*   neighbours) {
400
401        int i, k, k2, k3, k6;
402        long n;
403        double qmin, qmax, qn, qc;
404        double dq, dqa[3], phi, r;
405
406        for (k=0; k<N; k++){
407                k6 = 6*k;
408                k3 = 3*k;
409                k2 = 2*k;
410
411                qc = centroid_values[k];
412                phi = 1.0;
413
414                for (i=0; i<3; i++) {
415                    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
416                    dqa[i] = dq;                      //Save dqa for use in updating vertex values
417
418                    n = neighbours[k3+i];
419                    qn = qc;
420                    if (n >= 0)  qn = centroid_values[n]; //Neighbour's centroid value
421
422                    qmin = min(qc, qn);
423                    qmax = max(qc, qn);
424
425                    r = 1.0;
426     
427                    if (dq > 0.0) r = (qmax - qc)/dq;
428                    if (dq < 0.0) r = (qmin - qc)/dq;     
429                   
430                    phi = min( min(r*beta, 1.0), phi);   
431                   
432                }
433
434
435                //Update edge and vertex values using phi limiter
436                edge_values[k3+0] = qc + phi*dqa[0];
437                edge_values[k3+1] = qc + phi*dqa[1];
438                edge_values[k3+2] = qc + phi*dqa[2];
439               
440                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
441                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
442                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
443
444        }
445
446        return 0;
447}
448
449
450int _limit_gradient_by_neighbour(int N, double beta,
451                     double* centroid_values,
452                     double* vertex_values,
453                     double* edge_values,
454                     double* x_gradient,
455                     double* y_gradient,
456                     long*   neighbours) {
457
458        int i, k, k2, k3, k6;
459        long n;
460        double qmin, qmax, qn, qc;
461        double dq, dqa[3], phi, r;
462
463        for (k=0; k<N; k++){
464                k6 = 6*k;
465                k3 = 3*k;
466                k2 = 2*k;
467
468                qc = centroid_values[k];
469                phi = 1.0;
470
471                for (i=0; i<3; i++) {
472                    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
473                    dqa[i] = dq;                      //Save dq for use in updating vertex values
474
475                    n = neighbours[k3+i];
476                    if (n >= 0) {
477                        qn = centroid_values[n]; //Neighbour's centroid value
478
479                        qmin = min(qc, qn);
480                        qmax = max(qc, qn);
481
482                        r = 1.0;
483     
484                        if (dq > 0.0) r = (qmax - qc)/dq;
485                        if (dq < 0.0) r = (qmin - qc)/dq;     
486                   
487                        phi = min( min(r*beta, 1.0), phi);   
488                    }
489                }
490
491
492                //Update edge and vertex values using phi limiter
493                edge_values[k3+0] = qc + phi*dqa[0];
494                edge_values[k3+1] = qc + phi*dqa[1];
495                edge_values[k3+2] = qc + phi*dqa[2];
496               
497                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
498                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
499                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
500
501        }
502
503        return 0;
504}
505
506int _bound_vertices_below_by_constant(int N, double bound,
507                     double* centroid_values,
508                     double* vertex_values,
509                     double* edge_values,
510                     double* x_gradient,
511                     double* y_gradient) {
512
513        int i, k, k2, k3, k6;
514        double qmin, qc;
515        double dq, dqa[3], phi, r;
516
517        for (k=0; k<N; k++){
518                k6 = 6*k;
519                k3 = 3*k;
520                k2 = 2*k;
521
522                qc = centroid_values[k];
523                qmin = bound;
524
525
526                phi = 1.0;
527                for (i=0; i<3; i++) {   
528                    r = 1.0;
529     
530                    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
531                    dqa[i] = dq;                      //Save dq for use in updating vertex values
532     
533                    if (dq < 0.0) r = (qmin - qc)/dq;     
534 
535                   
536                    phi = min( min(r, 1.0), phi);   
537                }
538   
539
540                //Update gradient, vertex and edge values using phi limiter
541                x_gradient[k] = x_gradient[k]*phi;
542                y_gradient[k] = y_gradient[k]*phi;
543
544                vertex_values[k3+0] = qc + phi*dqa[0];
545                vertex_values[k3+1] = qc + phi*dqa[1];
546                vertex_values[k3+2] = qc + phi*dqa[2];
547               
548                edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
549                edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
550                edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
551
552        }
553
554        return 0;
555}
556
557int _bound_vertices_below_by_quantity(int N,
558                                      double* bound_vertex_values,
559                                      double* centroid_values,
560                                      double* vertex_values,
561                                      double* edge_values,
562                                      double* x_gradient,
563                                      double* y_gradient) {
564
565        int i, k, k2, k3, k6;
566        double qmin, qc;
567        double dq, dqa[3], phi, r;
568
569        for (k=0; k<N; k++){
570                k6 = 6*k;
571                k3 = 3*k;
572                k2 = 2*k;
573
574                qc = centroid_values[k];
575
576                phi = 1.0;
577                for (i=0; i<3; i++) {   
578                    r = 1.0;
579     
580                    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
581                    dqa[i] = dq;                      //Save dq for use in updating vertex values
582     
583                    qmin = bound_vertex_values[k3+i];
584                    if (dq < 0.0) r = (qmin - qc)/dq;     
585 
586                   
587                    phi = min( min(r, 1.0), phi);   
588                }
589   
590
591                //Update gradient, vertex and edge values using phi limiter
592                x_gradient[k] = x_gradient[k]*phi;
593                y_gradient[k] = y_gradient[k]*phi;
594
595                vertex_values[k3+0] = qc + phi*dqa[0];
596                vertex_values[k3+1] = qc + phi*dqa[1];
597                vertex_values[k3+2] = qc + phi*dqa[2];
598               
599                edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
600                edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
601                edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
602
603        }
604
605        return 0;
606}
607
608int _interpolate_from_vertices_to_edges(int N,
609                 double* vertex_values,
610                 double* edge_values) {
611
612        int k, k3;
613        double q0, q1, q2;
614
615
616        for (k=0; k<N; k++) {
617                k3 = 3*k;
618
619                q0 = vertex_values[k3 + 0];
620                q1 = vertex_values[k3 + 1];
621                q2 = vertex_values[k3 + 2];
622
623                edge_values[k3 + 0] = 0.5*(q1+q2);
624                edge_values[k3 + 1] = 0.5*(q0+q2);
625                edge_values[k3 + 2] = 0.5*(q0+q1);
626        }
627        return 0;
628}
629
630
631int _interpolate_from_edges_to_vertices(int N,
632                 double* vertex_values,
633                 double* edge_values) {
634
635        int k, k3;
636        double e0, e1, e2;
637
638
639        for (k=0; k<N; k++) {
640                k3 = 3*k;
641
642                e0 = edge_values[k3 + 0];
643                e1 = edge_values[k3 + 1];
644                e2 = edge_values[k3 + 2];
645
646                vertex_values[k3 + 0] = e1 + e2 - e0;
647                vertex_values[k3 + 1] = e2 + e0 - e1;
648                vertex_values[k3 + 2] = e0 + e1 - e2;
649        }
650        return 0;
651}
652
653int _backup_centroid_values(int N,
654                            double* centroid_values,
655                            double* centroid_backup_values) {
656    // Backup centroid values
657
658
659    int k;
660
661    for (k=0; k<N; k++) {
662        centroid_backup_values[k] = centroid_values[k];
663    }
664
665
666    return 0;
667}
668
669
670int _saxpy_centroid_values(int N,
671                           double a,
672                           double b,
673                           double* centroid_values,
674                           double* centroid_backup_values) {
675    // Saxby centroid values
676
677
678    int k;
679
680
681    for (k=0; k<N; k++) {
682        centroid_values[k] = a*centroid_values[k] + b*centroid_backup_values[k];
683    }
684
685
686    return 0;
687}
688
689
690int _update(int N,
691            double timestep,
692            double* centroid_values,
693            double* explicit_update,
694            double* semi_implicit_update) {
695        // Update centroid values based on values stored in
696        // explicit_update and semi_implicit_update as well as given timestep
697
698
699        int k;
700        double denominator, x;
701
702
703        // Divide semi_implicit update by conserved quantity
704        for (k=0; k<N; k++) {
705                x = centroid_values[k];
706                if (x == 0.0) {
707                        semi_implicit_update[k] = 0.0;
708                } else {
709                        semi_implicit_update[k] /= x;
710                }
711        }
712
713
714        // Semi implicit updates
715        for (k=0; k<N; k++) {
716                denominator = 1.0 - timestep*semi_implicit_update[k];
717                if (denominator == 0.0) {
718                        return -1;
719                } else {
720                        //Update conserved_quantities from semi implicit updates
721                        centroid_values[k] /= denominator;
722                }
723        }
724
725
726        // Explicit updates
727        for (k=0; k<N; k++) {
728                centroid_values[k] += timestep*explicit_update[k];
729        }
730
731
732        // MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
733        for (k=0;k<N;k++){
734                semi_implicit_update[k]=0.0;
735        }
736
737        return 0;
738}
739
740
741int _average_vertex_values(int N,
742                           long* vertex_value_indices,
743                           long* number_of_triangles_per_node,
744                           double* vertex_values, 
745                           double* A) {
746  // Average vertex values to obtain one value per node
747
748  int i, index; 
749  int k = 0; //Track triangles touching each node
750  int current_node = 0;
751  double total = 0.0;
752
753  for (i=0; i<N; i++) {
754 
755    // if (current_node == N) {
756    //   printf("Current node exceeding number of nodes (%d)", N);   
757    //   return 1;
758    // }
759   
760    index = vertex_value_indices[i];
761    k += 1;
762           
763    // volume_id = index / 3
764    // vertex_id = index % 3
765    // total += self.vertex_values[volume_id, vertex_id]
766    total += vertex_values[index];
767     
768    // printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);
769    if (number_of_triangles_per_node[current_node] == k) {
770      A[current_node] = total/k;
771               
772      // Move on to next node
773      total = 0.0;
774      k = 0;
775      current_node += 1;
776    }
777  }
778                           
779  return 0;
780}
781
782
783//-----------------------------------------------------
784// Python method Wrappers
785//-----------------------------------------------------
786
787PyObject *update(PyObject *self, PyObject *args) {
788  // FIXME (Ole): It would be great to turn this text into a Python DOC string
789
790  /*"""Update centroid values based on values stored in
791    explicit_update and semi_implicit_update as well as given timestep
792
793    Function implementing forcing terms must take on argument
794    which is the domain and they must update either explicit
795    or implicit updates, e,g,:
796
797    def gravity(domain):
798        ....
799        domain.quantities['xmomentum'].explicit_update = ...
800        domain.quantities['ymomentum'].explicit_update = ...
801
802
803
804    Explicit terms must have the form
805
806        G(q, t)
807
808    and explicit scheme is
809
810       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
811
812
813    Semi implicit forcing terms are assumed to have the form
814
815       G(q, t) = H(q, t) q
816
817    and the semi implicit scheme will then be
818
819      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
820
821  */
822
823        PyObject *quantity;
824        PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
825
826        double timestep;
827        int N, err;
828
829
830        // Convert Python arguments to C
831        if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep)) {
832          PyErr_SetString(PyExc_RuntimeError, 
833                          "quantity_ext.c: update could not parse input");
834          return NULL;
835        }
836
837        centroid_values      = get_consecutive_array(quantity, "centroid_values");
838        explicit_update      = get_consecutive_array(quantity, "explicit_update");
839        semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");
840
841        N = centroid_values -> dimensions[0];
842
843        err = _update(N, timestep,
844                      (double*) centroid_values -> data,
845                      (double*) explicit_update -> data,
846                      (double*) semi_implicit_update -> data);
847
848
849        if (err != 0) {
850          PyErr_SetString(PyExc_RuntimeError,
851                          "quantity_ext.c: update, divsion by zero in semi implicit update - call Stephen :)");
852          return NULL;
853        }
854
855        // Release and return
856        Py_DECREF(centroid_values);
857        Py_DECREF(explicit_update);
858        Py_DECREF(semi_implicit_update);
859
860        return Py_BuildValue("");
861}
862
863
864PyObject *backup_centroid_values(PyObject *self, PyObject *args) {
865
866        PyObject *quantity;
867        PyArrayObject *centroid_values, *centroid_backup_values;
868
869        int N, err;
870
871
872        // Convert Python arguments to C
873        if (!PyArg_ParseTuple(args, "O", &quantity)) {
874          PyErr_SetString(PyExc_RuntimeError, 
875                          "quantity_ext.c: backup_centroid_values could not parse input");
876          return NULL;
877        }
878
879        centroid_values        = get_consecutive_array(quantity, "centroid_values");
880        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
881
882        N = centroid_values -> dimensions[0];
883
884        err = _backup_centroid_values(N,
885                      (double*) centroid_values -> data,
886                      (double*) centroid_backup_values -> data);
887
888
889        // Release and return
890        Py_DECREF(centroid_values);
891        Py_DECREF(centroid_backup_values);
892
893        return Py_BuildValue("");
894}
895
896PyObject *saxpy_centroid_values(PyObject *self, PyObject *args) {
897
898        PyObject *quantity;
899        PyArrayObject *centroid_values, *centroid_backup_values;
900
901        double a,b;
902        int N, err;
903
904
905        // Convert Python arguments to C
906        if (!PyArg_ParseTuple(args, "Odd", &quantity, &a, &b)) {
907          PyErr_SetString(PyExc_RuntimeError, 
908                          "quantity_ext.c: saxpy_centroid_values could not parse input");
909          return NULL;
910        }
911
912        centroid_values        = get_consecutive_array(quantity, "centroid_values");
913        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
914
915        N = centroid_values -> dimensions[0];
916
917        err = _saxpy_centroid_values(N,a,b,
918                      (double*) centroid_values -> data,
919                      (double*) centroid_backup_values -> data);
920
921
922        // Release and return
923        Py_DECREF(centroid_values);
924        Py_DECREF(centroid_backup_values);
925
926        return Py_BuildValue("");
927}
928
929
930PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
931        //
932        //Compute edge values from vertex values using linear interpolation
933        //
934       
935        PyObject *quantity;
936        PyArrayObject *vertex_values, *edge_values;
937
938        int N, err;
939
940        // Convert Python arguments to C
941        if (!PyArg_ParseTuple(args, "O", &quantity)) {
942          PyErr_SetString(PyExc_RuntimeError, 
943                          "quantity_ext.c: interpolate_from_vertices_to_edges could not parse input");
944          return NULL;
945        }
946       
947        vertex_values = get_consecutive_array(quantity, "vertex_values");
948        edge_values = get_consecutive_array(quantity, "edge_values");
949
950        N = vertex_values -> dimensions[0];
951
952        err = _interpolate_from_vertices_to_edges(N,
953                           (double*) vertex_values -> data,
954                           (double*) edge_values -> data);
955
956        if (err != 0) {
957          PyErr_SetString(PyExc_RuntimeError, 
958                          "Interpolate could not be computed");
959          return NULL;
960        }
961
962        // Release and return
963        Py_DECREF(vertex_values);
964        Py_DECREF(edge_values);
965
966        return Py_BuildValue("");
967}
968
969
970PyObject *interpolate_from_edges_to_vertices(PyObject *self, PyObject *args) {
971        //
972        //Compute vertex values from edge values using linear interpolation
973        //
974       
975        PyObject *quantity;
976        PyArrayObject *vertex_values, *edge_values;
977
978        int N, err;
979
980        // Convert Python arguments to C
981        if (!PyArg_ParseTuple(args, "O", &quantity)) {
982          PyErr_SetString(PyExc_RuntimeError, 
983                          "quantity_ext.c: interpolate_from_edges_to_vertices could not parse input");
984          return NULL;
985        }
986       
987        vertex_values = get_consecutive_array(quantity, "vertex_values");
988        edge_values = get_consecutive_array(quantity, "edge_values");
989
990        N = vertex_values -> dimensions[0];
991
992        err = _interpolate_from_edges_to_vertices(N,
993                           (double*) vertex_values -> data,
994                           (double*) edge_values -> data);
995
996        if (err != 0) {
997          PyErr_SetString(PyExc_RuntimeError, 
998                          "Interpolate could not be computed");
999          return NULL;
1000        }
1001
1002        // Release and return
1003        Py_DECREF(vertex_values);
1004        Py_DECREF(edge_values);
1005
1006        return Py_BuildValue("");
1007}
1008
1009
1010PyObject *average_vertex_values(PyObject *self, PyObject *args) {
1011
1012        PyArrayObject
1013          *vertex_value_indices, 
1014          *number_of_triangles_per_node,
1015          *vertex_values,
1016          *A;
1017       
1018
1019        int N, err;
1020
1021        // Convert Python arguments to C
1022        if (!PyArg_ParseTuple(args, "OOOO",
1023                              &vertex_value_indices, 
1024                              &number_of_triangles_per_node,
1025                              &vertex_values, 
1026                              &A)) {
1027          PyErr_SetString(PyExc_RuntimeError, 
1028                          "quantity_ext.c: average_vertex_values could not parse input");
1029          return NULL;
1030        }
1031       
[6410]1032        // check that numpy array objects arrays are C contiguous memory
1033        CHECK_C_CONTIG(vertex_value_indices);
1034        CHECK_C_CONTIG(number_of_triangles_per_node);
1035        CHECK_C_CONTIG(vertex_values);
1036        CHECK_C_CONTIG(A);
1037
[5897]1038        N = vertex_value_indices -> dimensions[0];
1039        // printf("Got parameters, N=%d\n", N);
1040        err = _average_vertex_values(N,
1041                                     (long*) vertex_value_indices -> data,
1042                                     (long*) number_of_triangles_per_node -> data,
1043                                     (double*) vertex_values -> data, 
1044                                     (double*) A -> data);
1045
1046        //printf("Error %d", err);
1047        if (err != 0) {
1048          PyErr_SetString(PyExc_RuntimeError, 
1049                          "average_vertex_values could not be computed");
1050          return NULL;
1051        }
1052
1053        return Py_BuildValue("");
1054}
1055
1056
1057
1058PyObject *extrapolate_from_gradient(PyObject *self, PyObject *args) {
1059
1060        PyObject *quantity, *domain;
1061        PyArrayObject
1062            *centroids,            //Coordinates at centroids
1063            *centroid_values,      //Values at centroids
1064            *vertex_coordinates,   //Coordinates at vertices
1065            *vertex_values,        //Values at vertices
1066            *edge_values,          //Values at edges
1067            *number_of_boundaries, //Number of boundaries for each triangle
1068            *surrogate_neighbours, //True neighbours or - if one missing - self
1069            *x_gradient,           //x gradient
1070            *y_gradient;           //y gradient
1071
1072        //int N, err;
1073        //int dimensions[1];
1074        int N, err;
1075        //double *a, *b;  //Gradients
1076
1077        // Convert Python arguments to C
1078        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1079          PyErr_SetString(PyExc_RuntimeError, 
1080                          "extrapolate_gradient could not parse input");       
1081          return NULL;
1082        }
1083
1084        domain = PyObject_GetAttrString(quantity, "domain");
1085        if (!domain) {
1086          PyErr_SetString(PyExc_RuntimeError, 
1087                          "extrapolate_gradient could not obtain domain object from quantity"); 
1088          return NULL;
1089        }
1090
1091        // Get pertinent variables
1092        centroids            = get_consecutive_array(domain,   "centroid_coordinates");
1093        centroid_values      = get_consecutive_array(quantity, "centroid_values");
1094        surrogate_neighbours = get_consecutive_array(domain,   "surrogate_neighbours");
1095        number_of_boundaries = get_consecutive_array(domain,   "number_of_boundaries");
1096        vertex_coordinates   = get_consecutive_array(domain,   "vertex_coordinates");
1097        vertex_values        = get_consecutive_array(quantity, "vertex_values");
1098        edge_values          = get_consecutive_array(quantity, "edge_values");
1099        x_gradient           = get_consecutive_array(quantity, "x_gradient");
1100        y_gradient           = get_consecutive_array(quantity, "y_gradient");
1101
1102        N = centroid_values -> dimensions[0];
1103
1104        // Release
1105        Py_DECREF(domain);
1106
1107        err = _extrapolate_from_gradient(N,
1108                        (double*) centroids -> data,
1109                        (double*) centroid_values -> data,
1110                        (double*) vertex_coordinates -> data,
1111                        (double*) vertex_values -> data,
1112                        (double*) edge_values -> data,
1113                        (double*) x_gradient -> data,
1114                        (double*) y_gradient -> data);
1115
1116
1117        if (err != 0) {
1118          PyErr_SetString(PyExc_RuntimeError,
1119                          "Internal function _extrapolate failed");
1120          return NULL;
1121        }
1122
1123
1124
1125        // Release
1126        Py_DECREF(centroids);
1127        Py_DECREF(centroid_values);
1128        Py_DECREF(number_of_boundaries);
1129        Py_DECREF(surrogate_neighbours);
1130        Py_DECREF(vertex_coordinates);
1131        Py_DECREF(vertex_values);
1132        Py_DECREF(edge_values);
1133        Py_DECREF(x_gradient);
1134        Py_DECREF(y_gradient);
1135
1136        return Py_BuildValue("");
1137}
1138
1139
1140PyObject *extrapolate_second_order_and_limit_by_edge(PyObject *self, PyObject *args) {
1141    /* Compute edge values using second order approximation and limit values
1142       so that edge values are limited by the two corresponding centroid values
1143       
1144       Python Call:
1145       extrapolate_second_order_and_limit(domain,quantity,beta)
1146    */
1147
1148    PyObject *quantity, *domain;
1149       
1150    PyArrayObject
1151        *domain_centroids,              //Coordinates at centroids
1152        *domain_vertex_coordinates,     //Coordinates at vertices
1153        *domain_number_of_boundaries,   //Number of boundaries for each triangle
1154        *domain_surrogate_neighbours,   //True neighbours or - if one missing - self
1155        *domain_neighbours,             //True neighbours, or if negative a link to boundary
1156
1157        *quantity_centroid_values,      //Values at centroids   
1158        *quantity_vertex_values,        //Values at vertices
1159        *quantity_edge_values,          //Values at edges
1160        *quantity_phi,                  //limiter phi values
1161        *quantity_x_gradient,           //x gradient
1162        *quantity_y_gradient;           //y gradient
1163   
1164
1165  // Local variables
1166  int ntri;
1167  double beta;
1168  int err;
1169   
1170  // Convert Python arguments to C
1171  if (!PyArg_ParseTuple(args, "O",&quantity)) {
1172      PyErr_SetString(PyExc_RuntimeError, 
1173                      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");     
1174      return NULL;
1175  }
1176
1177  domain = PyObject_GetAttrString(quantity, "domain");
1178  if (!domain) {
1179      PyErr_SetString(PyExc_RuntimeError, 
1180                      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");       
1181      return NULL;
1182  }
1183
1184       
1185  // Get pertinent variables
1186  domain_centroids              = get_consecutive_array(domain,   "centroid_coordinates");
1187  domain_surrogate_neighbours   = get_consecutive_array(domain,   "surrogate_neighbours");
1188  domain_number_of_boundaries   = get_consecutive_array(domain,   "number_of_boundaries");
1189  domain_vertex_coordinates     = get_consecutive_array(domain,   "vertex_coordinates");
1190  domain_neighbours             = get_consecutive_array(domain,   "neighbours");
1191
1192  quantity_centroid_values      = get_consecutive_array(quantity, "centroid_values");
1193  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
1194  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
1195  quantity_phi                  = get_consecutive_array(quantity, "phi");
1196  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
1197  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
1198
1199  beta = get_python_double(quantity,"beta");
1200
1201  ntri = quantity_centroid_values -> dimensions[0];
1202
1203  err = _compute_gradients(ntri,
1204                           (double*) domain_centroids -> data,
1205                           (double*) quantity_centroid_values -> data,
1206                           (long*)   domain_number_of_boundaries -> data,
1207                           (long*)   domain_surrogate_neighbours -> data,
1208                           (double*) quantity_x_gradient -> data,
1209                           (double*) quantity_y_gradient -> data);
1210
1211  if (err != 0) {
1212      PyErr_SetString(PyExc_RuntimeError,
1213                          "quantity_ext.c: Internal function _compute_gradient failed");
1214      return NULL;
1215  }
1216
1217
1218  err = _extrapolate_from_gradient(ntri, 
1219                                   (double*) domain_centroids -> data,
1220                                   (double*) quantity_centroid_values -> data,
1221                                   (double*) domain_vertex_coordinates -> data,
1222                                   (double*) quantity_vertex_values -> data,
1223                                   (double*) quantity_edge_values -> data,
1224                                   (double*) quantity_x_gradient -> data,
1225                                   (double*) quantity_y_gradient -> data);
1226
1227  if (err != 0) {
1228      PyErr_SetString(PyExc_RuntimeError,
1229                          "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
1230      return NULL;
1231  }
1232
1233
1234  err = _limit_edges_by_all_neighbours(ntri, beta,
1235                                       (double*) quantity_centroid_values -> data,
1236                                       (double*) quantity_vertex_values -> data,
1237                                       (double*) quantity_edge_values -> data,
1238                                       (long*)   domain_neighbours -> data,
1239                                       (double*) quantity_x_gradient -> data,
1240                                       (double*) quantity_y_gradient -> data);
1241
1242  if (err != 0) {
1243      PyErr_SetString(PyExc_RuntimeError,
1244                          "quantity_ext.c: Internal function _limit_edges_by_all_neighbours failed");
1245      return NULL;
1246  }
1247
1248
1249  // Release
1250  Py_DECREF(domain_centroids);
1251  Py_DECREF(domain_surrogate_neighbours);
1252  Py_DECREF(domain_number_of_boundaries);
1253  Py_DECREF(domain_vertex_coordinates);
1254
1255  Py_DECREF(quantity_centroid_values);
1256  Py_DECREF(quantity_vertex_values);
1257  Py_DECREF(quantity_edge_values);
1258  Py_DECREF(quantity_phi);
1259  Py_DECREF(quantity_x_gradient);
1260  Py_DECREF(quantity_y_gradient);
1261
1262  return Py_BuildValue("");
1263}
1264
1265
1266PyObject *extrapolate_second_order_and_limit_by_vertex(PyObject *self, PyObject *args) {
1267    /* Compute edge values using second order approximation and limit values
1268       so that edge values are limited by the two corresponding centroid values
1269       
1270       Python Call:
1271       extrapolate_second_order_and_limit(domain,quantity,beta)
1272    */
1273
1274    PyObject *quantity, *domain;
1275       
1276    PyArrayObject
1277        *domain_centroids,              //Coordinates at centroids
1278        *domain_vertex_coordinates,     //Coordinates at vertices
1279        *domain_number_of_boundaries,   //Number of boundaries for each triangle
1280        *domain_surrogate_neighbours,   //True neighbours or - if one missing - self
1281        *domain_neighbours,             //True neighbours, or if negative a link to boundary
1282
1283        *quantity_centroid_values,      //Values at centroids   
1284        *quantity_vertex_values,        //Values at vertices
1285        *quantity_edge_values,          //Values at edges
1286        *quantity_phi,                  //limiter phi values
1287        *quantity_x_gradient,           //x gradient
1288        *quantity_y_gradient;           //y gradient
1289   
1290
1291  // Local variables
1292  int ntri;
1293  double beta;
1294  int err;
1295   
1296  // Convert Python arguments to C
1297  if (!PyArg_ParseTuple(args, "O",&quantity)) {
1298      PyErr_SetString(PyExc_RuntimeError, 
1299                      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");     
1300      return NULL;
1301  }
1302
1303  domain = PyObject_GetAttrString(quantity, "domain");
1304  if (!domain) {
1305      PyErr_SetString(PyExc_RuntimeError, 
1306                      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");       
1307      return NULL;
1308  }
1309
1310       
1311  // Get pertinent variables
1312  domain_centroids              = get_consecutive_array(domain,   "centroid_coordinates");
1313  domain_surrogate_neighbours   = get_consecutive_array(domain,   "surrogate_neighbours");
1314  domain_number_of_boundaries   = get_consecutive_array(domain,   "number_of_boundaries");
1315  domain_vertex_coordinates     = get_consecutive_array(domain,   "vertex_coordinates");
1316  domain_neighbours             = get_consecutive_array(domain,   "neighbours");
1317
1318  quantity_centroid_values      = get_consecutive_array(quantity, "centroid_values");
1319  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
1320  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
1321  quantity_phi                  = get_consecutive_array(quantity, "phi");
1322  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
1323  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
1324
1325  beta = get_python_double(quantity,"beta");
1326
1327  ntri = quantity_centroid_values -> dimensions[0];
1328
1329  err = _compute_gradients(ntri,
1330                           (double*) domain_centroids -> data,
1331                           (double*) quantity_centroid_values -> data,
1332                           (long*)   domain_number_of_boundaries -> data,
1333                           (long*)   domain_surrogate_neighbours -> data,
1334                           (double*) quantity_x_gradient -> data,
1335                           (double*) quantity_y_gradient -> data);
1336
1337  if (err != 0) {
1338      PyErr_SetString(PyExc_RuntimeError,
1339                          "quantity_ext.c: Internal function _compute_gradient failed");
1340      return NULL;
1341  }
1342
1343
1344  err = _extrapolate_from_gradient(ntri, 
1345                                   (double*) domain_centroids -> data,
1346                                   (double*) quantity_centroid_values -> data,
1347                                   (double*) domain_vertex_coordinates -> data,
1348                                   (double*) quantity_vertex_values -> data,
1349                                   (double*) quantity_edge_values -> data,
1350                                   (double*) quantity_x_gradient -> data,
1351                                   (double*) quantity_y_gradient -> data);
1352
1353  if (err != 0) {
1354      PyErr_SetString(PyExc_RuntimeError,
1355                          "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
1356      return NULL;
1357  }
1358
1359
1360  err = _limit_vertices_by_all_neighbours(ntri, beta,
1361                                       (double*) quantity_centroid_values -> data,
1362                                       (double*) quantity_vertex_values -> data,
1363                                       (double*) quantity_edge_values -> data,
1364                                       (long*)   domain_neighbours -> data,
1365                                       (double*) quantity_x_gradient -> data,
1366                                       (double*) quantity_y_gradient -> data);
1367
1368  if (err != 0) {
1369      PyErr_SetString(PyExc_RuntimeError,
1370                          "quantity_ext.c: Internal function _limit_vertices_by_all_neighbours failed");
1371      return NULL;
1372  }
1373
1374
1375  // Release
1376  Py_DECREF(domain_centroids);
1377  Py_DECREF(domain_surrogate_neighbours);
1378  Py_DECREF(domain_number_of_boundaries);
1379  Py_DECREF(domain_vertex_coordinates);
1380
1381  Py_DECREF(quantity_centroid_values);
1382  Py_DECREF(quantity_vertex_values);
1383  Py_DECREF(quantity_edge_values);
1384  Py_DECREF(quantity_phi);
1385  Py_DECREF(quantity_x_gradient);
1386  Py_DECREF(quantity_y_gradient);
1387
1388  return Py_BuildValue("");
1389}
1390
1391
1392
1393PyObject *compute_gradients(PyObject *self, PyObject *args) {
1394
1395        PyObject *quantity, *domain;
1396        PyArrayObject
1397            *centroids,            //Coordinates at centroids
1398            *centroid_values,      //Values at centroids
1399            *vertex_coordinates,   //Coordinates at vertices
1400            *vertex_values,        //Values at vertices
1401            *edge_values,          //Values at edges
1402            *number_of_boundaries, //Number of boundaries for each triangle
1403            *surrogate_neighbours, //True neighbours or - if one missing - self
1404            *x_gradient,           //x gradient
1405            *y_gradient;           //y gradient
1406
1407        //int N, err;
1408        //int dimensions[1];
1409        int N, err;
1410        //double *a, *b;  //Gradients
1411
1412        // Convert Python arguments to C
1413        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1414          PyErr_SetString(PyExc_RuntimeError, 
1415                          "compute_gradients could not parse input");   
1416          return NULL;
1417        }
1418
1419        domain = PyObject_GetAttrString(quantity, "domain");
1420        if (!domain) {
1421          PyErr_SetString(PyExc_RuntimeError, 
1422                          "compute_gradients could not obtain domain object from quantity");   
1423          return NULL;
1424        }
1425
1426        // Get pertinent variables
1427        centroids            = get_consecutive_array(domain,   "centroid_coordinates");
1428        centroid_values      = get_consecutive_array(quantity, "centroid_values");
1429        surrogate_neighbours = get_consecutive_array(domain,   "surrogate_neighbours");
1430        number_of_boundaries = get_consecutive_array(domain,   "number_of_boundaries");
1431        vertex_coordinates   = get_consecutive_array(domain,   "vertex_coordinates");
1432        vertex_values        = get_consecutive_array(quantity, "vertex_values");
1433        edge_values          = get_consecutive_array(quantity, "edge_values");
1434        x_gradient           = get_consecutive_array(quantity, "x_gradient");
1435        y_gradient           = get_consecutive_array(quantity, "y_gradient");
1436
1437        N = centroid_values -> dimensions[0];
1438
1439        // Release
1440        Py_DECREF(domain);
1441
1442
1443        err = _compute_gradients(N,
1444                        (double*) centroids -> data,
1445                        (double*) centroid_values -> data,
1446                        (long*) number_of_boundaries -> data,
1447                        (long*) surrogate_neighbours -> data,
1448                        (double*) x_gradient -> data,
1449                        (double*) y_gradient -> data);
1450
1451        if (err != 0) {
1452          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
1453          return NULL;
1454        }
1455
1456
1457
1458        // Release
1459        Py_DECREF(centroids);
1460        Py_DECREF(centroid_values);
1461        Py_DECREF(number_of_boundaries);
1462        Py_DECREF(surrogate_neighbours);
1463        Py_DECREF(vertex_coordinates);
1464        Py_DECREF(vertex_values);
1465        Py_DECREF(edge_values);
1466        Py_DECREF(x_gradient);
1467        Py_DECREF(y_gradient);
1468
1469        return Py_BuildValue("");
1470}
1471
1472
1473
1474PyObject *limit_old(PyObject *self, PyObject *args) {
1475  //Limit slopes for each volume to eliminate artificial variance
1476  //introduced by e.g. second order extrapolator
1477
1478  //This is an unsophisticated limiter as it does not take into
1479  //account dependencies among quantities.
1480
1481  //precondition:
1482  //  vertex values are estimated from gradient
1483  //postcondition:
1484  //  vertex values are updated
1485  //
1486
1487        PyObject *quantity, *domain, *Tmp;
1488        PyArrayObject
1489            *qv, //Conserved quantities at vertices
1490            *qc, //Conserved quantities at centroids
1491            *neighbours;
1492
1493        int k, i, n, N, k3;
1494        double beta_w; //Safety factor
1495        double *qmin, *qmax, qn;
1496
1497        // Convert Python arguments to C
1498        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1499          PyErr_SetString(PyExc_RuntimeError, 
1500                          "quantity_ext.c: limit_old could not parse input");
1501          return NULL;
1502        }
1503
1504        domain = PyObject_GetAttrString(quantity, "domain");
1505        if (!domain) {
1506          PyErr_SetString(PyExc_RuntimeError, 
1507                          "quantity_ext.c: limit_old could not obtain domain object from quantity");                   
1508         
1509          return NULL;
1510        }
1511
1512        //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
1513        neighbours = get_consecutive_array(domain, "neighbours");
1514
1515        // Get safety factor beta_w
1516        Tmp = PyObject_GetAttrString(domain, "beta_w");
1517        if (!Tmp) {
1518          PyErr_SetString(PyExc_RuntimeError, 
1519                          "quantity_ext.c: limit_old could not obtain beta_w object from domain");                     
1520         
1521          return NULL;
1522        }       
1523
1524        beta_w = PyFloat_AsDouble(Tmp);
1525
1526        Py_DECREF(Tmp);
1527        Py_DECREF(domain);
1528
1529
1530        qc = get_consecutive_array(quantity, "centroid_values");
1531        qv = get_consecutive_array(quantity, "vertex_values");
1532
1533
1534        N = qc -> dimensions[0];
1535
1536        // Find min and max of this and neighbour's centroid values
1537        qmin = malloc(N * sizeof(double));
1538        qmax = malloc(N * sizeof(double));
1539        for (k=0; k<N; k++) {
1540                k3=k*3;
1541
1542                qmin[k] = ((double*) qc -> data)[k];
1543                qmax[k] = qmin[k];
1544
1545                for (i=0; i<3; i++) {
1546                        n = ((long*) neighbours -> data)[k3+i];
1547                        if (n >= 0) {
1548                                qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
1549
1550                                qmin[k] = min(qmin[k], qn);
1551                                qmax[k] = max(qmax[k], qn);
1552                        }
1553                        //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]);
1554                        //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]);
1555                }
1556        }
1557
1558        // Call underlying routine
1559        _limit_old(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
1560
1561        free(qmin);
1562        free(qmax);
1563        return Py_BuildValue("");
1564}
1565
1566
1567PyObject *limit_vertices_by_all_neighbours(PyObject *self, PyObject *args) {
1568  //Limit slopes for each volume to eliminate artificial variance
1569  //introduced by e.g. second order extrapolator
1570
1571  //This is an unsophisticated limiter as it does not take into
1572  //account dependencies among quantities.
1573
1574  //precondition:
1575  //  vertex values are estimated from gradient
1576  //postcondition:
1577  //  vertex and edge values are updated
1578  //
1579
1580        PyObject *quantity, *domain, *Tmp;
1581        PyArrayObject
1582            *vertex_values,   //Conserved quantities at vertices
1583            *centroid_values, //Conserved quantities at centroids
1584            *edge_values,     //Conserved quantities at edges
1585            *neighbours,
1586            *x_gradient,
1587            *y_gradient;
1588
1589        double beta_w; //Safety factor
1590        int N, err;
1591
1592
1593        // Convert Python arguments to C
1594        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1595          PyErr_SetString(PyExc_RuntimeError, 
1596                          "quantity_ext.c: limit_by_vertex could not parse input");
1597          return NULL;
1598        }
1599
1600        domain = PyObject_GetAttrString(quantity, "domain");
1601        if (!domain) {
1602          PyErr_SetString(PyExc_RuntimeError, 
1603                          "quantity_ext.c: limit_by_vertex could not obtain domain object from quantity");                     
1604         
1605          return NULL;
1606        }
1607
1608        // Get safety factor beta_w
1609        Tmp = PyObject_GetAttrString(domain, "beta_w");
1610        if (!Tmp) {
1611          PyErr_SetString(PyExc_RuntimeError, 
1612                          "quantity_ext.c: limit_by_vertex could not obtain beta_w object from domain");                       
1613         
1614          return NULL;
1615        }       
1616
1617
1618        // Get pertinent variables
1619        neighbours       = get_consecutive_array(domain, "neighbours");
1620        centroid_values  = get_consecutive_array(quantity, "centroid_values");
1621        vertex_values    = get_consecutive_array(quantity, "vertex_values");
1622        edge_values      = get_consecutive_array(quantity, "edge_values");
1623        x_gradient       = get_consecutive_array(quantity, "x_gradient");
1624        y_gradient       = get_consecutive_array(quantity, "y_gradient");
1625        beta_w           = get_python_double(domain,"beta_w");
1626
1627
1628
1629        N = centroid_values -> dimensions[0];
1630
1631        err = _limit_vertices_by_all_neighbours(N, beta_w,
1632                                                (double*) centroid_values -> data,
1633                                                (double*) vertex_values -> data,
1634                                                (double*) edge_values -> data,
1635                                                (long*)   neighbours -> data,
1636                                                (double*) x_gradient -> data,
1637                                                (double*) y_gradient -> data);
1638
1639
1640       
1641        if (err != 0) {
1642          PyErr_SetString(PyExc_RuntimeError,
1643                          "Internal function _limit_by_vertex failed");
1644          return NULL;
1645        }       
1646
1647
1648        // Release
1649        Py_DECREF(neighbours);
1650        Py_DECREF(centroid_values);
1651        Py_DECREF(vertex_values);
1652        Py_DECREF(edge_values);
1653        Py_DECREF(x_gradient);
1654        Py_DECREF(y_gradient);
1655        Py_DECREF(Tmp);
1656
1657
1658        return Py_BuildValue("");
1659}
1660
1661
1662
1663PyObject *limit_edges_by_all_neighbours(PyObject *self, PyObject *args) {
1664  //Limit slopes for each volume to eliminate artificial variance
1665  //introduced by e.g. second order extrapolator
1666
1667  //This is an unsophisticated limiter as it does not take into
1668  //account dependencies among quantities.
1669
1670  //precondition:
1671  //  vertex values are estimated from gradient
1672  //postcondition:
1673  //  vertex and edge values are updated
1674  //
1675
1676        PyObject *quantity, *domain;
1677        PyArrayObject
1678            *vertex_values,   //Conserved quantities at vertices
1679            *centroid_values, //Conserved quantities at centroids
1680            *edge_values,     //Conserved quantities at edges
1681            *x_gradient,
1682            *y_gradient,
1683            *neighbours;
1684
1685        double beta_w; //Safety factor
1686        int N, err;
1687
1688
1689        // Convert Python arguments to C
1690        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1691          PyErr_SetString(PyExc_RuntimeError, 
1692                          "quantity_ext.c: limit_edges_by_all_neighbours could not parse input");
1693          return NULL;
1694        }
1695
1696        domain = get_python_object(quantity, "domain");
1697
1698
1699        // Get pertinent variables
1700        neighbours       = get_consecutive_array(domain, "neighbours");
1701        centroid_values  = get_consecutive_array(quantity, "centroid_values");
1702        vertex_values    = get_consecutive_array(quantity, "vertex_values");
1703        edge_values      = get_consecutive_array(quantity, "edge_values");
1704        x_gradient       = get_consecutive_array(quantity, "x_gradient");
1705        y_gradient       = get_consecutive_array(quantity, "y_gradient");
1706        beta_w           = get_python_double(domain,"beta_w");
1707
1708
1709
1710        N = centroid_values -> dimensions[0];
1711
1712        err = _limit_edges_by_all_neighbours(N, beta_w,
1713                                             (double*) centroid_values -> data,
1714                                             (double*) vertex_values -> data,
1715                                             (double*) edge_values -> data,
1716                                             (long*)   neighbours -> data,
1717                                             (double*) x_gradient -> data,
1718                                             (double*) y_gradient -> data);
1719
1720        if (err != 0) {
1721          PyErr_SetString(PyExc_RuntimeError,
1722                          "quantity_ect.c: limit_edges_by_all_neighbours internal function _limit_edges_by_all_neighbours failed");
1723          return NULL;
1724        }       
1725
1726
1727        // Release
1728        Py_DECREF(neighbours);
1729        Py_DECREF(centroid_values);
1730        Py_DECREF(vertex_values);
1731        Py_DECREF(edge_values);
1732        Py_DECREF(x_gradient);
1733        Py_DECREF(y_gradient);
1734
1735
1736
1737        return Py_BuildValue("");
1738}
1739
1740PyObject *bound_vertices_below_by_constant(PyObject *self, PyObject *args) {
1741  //Bound a quantity below by a contant (useful for ensuring positivity
1742  //precondition:
1743  //  vertex values are already calulated, gradient consistent
1744  //postcondition:
1745  //  gradient, vertex and edge values are updated
1746  //
1747
1748        PyObject *quantity, *domain;
1749        PyArrayObject
1750            *vertex_values,   //Conserved quantities at vertices
1751            *centroid_values, //Conserved quantities at centroids
1752            *edge_values,     //Conserved quantities at edges
1753            *x_gradient,
1754            *y_gradient;
1755
1756        double bound; //Safety factor
1757        int N, err;
1758
1759
1760        // Convert Python arguments to C
1761        if (!PyArg_ParseTuple(args, "Od", &quantity, &bound)) {
1762          PyErr_SetString(PyExc_RuntimeError, 
1763                          "quantity_ext.c: bound_vertices_below_by_constant could not parse input");
1764          return NULL;
1765        }
1766
1767        domain = get_python_object(quantity, "domain");
1768
1769        // Get pertinent variables
1770        centroid_values  = get_consecutive_array(quantity, "centroid_values");
1771        vertex_values    = get_consecutive_array(quantity, "vertex_values");
1772        edge_values      = get_consecutive_array(quantity, "edge_values");
1773        x_gradient       = get_consecutive_array(quantity, "x_gradient");
1774        y_gradient       = get_consecutive_array(quantity, "y_gradient");
1775
1776
1777
1778
1779        N = centroid_values -> dimensions[0];
1780
1781        err = _bound_vertices_below_by_constant(N, bound,
1782                                             (double*) centroid_values -> data,
1783                                             (double*) vertex_values -> data,
1784                                             (double*) edge_values -> data,
1785                                             (double*) x_gradient -> data,
1786                                             (double*) y_gradient -> data);
1787
1788        if (err != 0) {
1789          PyErr_SetString(PyExc_RuntimeError,
1790                          "quantity_ect.c: bound_vertices_below_by_constant internal function _bound_vertices_below_by_constant failed");
1791          return NULL;
1792        }       
1793
1794
1795        // Release
1796        Py_DECREF(centroid_values);
1797        Py_DECREF(vertex_values);
1798        Py_DECREF(edge_values);
1799        Py_DECREF(x_gradient);
1800        Py_DECREF(y_gradient);
1801
1802
1803
1804        return Py_BuildValue("");
1805}
1806
1807
1808PyObject *bound_vertices_below_by_quantity(PyObject *self, PyObject *args) {
1809  //Bound a quantity below by a contant (useful for ensuring positivity
1810  //precondition:
1811  //  vertex values are already calulated, gradient consistent
1812  //postcondition:
1813  //  gradient, vertex and edge values are updated
1814  //
1815
1816        PyObject *quantity, *bounding_quantity, *domain;
1817        PyArrayObject
1818            *vertex_values,   //Conserved quantities at vertices
1819            *centroid_values, //Conserved quantities at centroids
1820            *edge_values,     //Conserved quantities at edges
1821            *x_gradient,
1822            *y_gradient,
1823            *bound_vertex_values;
1824
1825        int N, err;
1826
1827
1828        // Convert Python arguments to C
1829        if (!PyArg_ParseTuple(args, "OO", &quantity, &bounding_quantity)) {
1830          PyErr_SetString(PyExc_RuntimeError, 
1831                          "quantity_ext.c: bound_vertices_below_by_quantity could not parse input");
1832          return NULL;
1833        }
1834
1835        domain = get_python_object(quantity, "domain");
1836
1837        // Get pertinent variables
1838        centroid_values     = get_consecutive_array(quantity, "centroid_values");
1839        vertex_values       = get_consecutive_array(quantity, "vertex_values");
1840        edge_values         = get_consecutive_array(quantity, "edge_values");
1841        x_gradient          = get_consecutive_array(quantity, "x_gradient");
1842        y_gradient          = get_consecutive_array(quantity, "y_gradient");
1843        bound_vertex_values = get_consecutive_array(bounding_quantity, "vertex_values");
1844
1845
1846
1847        Py_DECREF(domain);
1848
1849        N = centroid_values -> dimensions[0];
1850
1851        err = _bound_vertices_below_by_quantity(N, 
1852                                                (double*) bound_vertex_values -> data,
1853                                                (double*) centroid_values -> data,
1854                                                (double*) vertex_values -> data,
1855                                                (double*) edge_values -> data,
1856                                                (double*) x_gradient -> data,
1857                                                (double*) y_gradient -> data);
1858
1859        if (err != 0) {
1860          PyErr_SetString(PyExc_RuntimeError,
1861                          "quantity_ect.c: bound_vertices_below_by_quantity internal function _bound_vertices_below_by_quantity failed");
1862          return NULL;
1863        }       
1864
1865
1866        // Release
1867        Py_DECREF(centroid_values);
1868        Py_DECREF(vertex_values);
1869        Py_DECREF(edge_values);
1870        Py_DECREF(x_gradient);
1871        Py_DECREF(y_gradient);
1872        Py_DECREF(bound_vertex_values);
1873
1874
1875
1876        return Py_BuildValue("");
1877}
1878
1879
1880PyObject *limit_edges_by_neighbour(PyObject *self, PyObject *args) {
1881  //Limit slopes for each volume to eliminate artificial variance
1882  //introduced by e.g. second order extrapolator
1883
1884  //This is an unsophisticated limiter as it does not take into
1885  //account dependencies among quantities.
1886
1887  //precondition:
1888  //  vertex values are estimated from gradient
1889  //postcondition:
1890  //  vertex and edge values are updated
1891  //
1892
1893        PyObject *quantity, *domain, *Tmp;
1894        PyArrayObject
1895            *vertex_values,   //Conserved quantities at vertices
1896            *centroid_values, //Conserved quantities at centroids
1897            *edge_values,     //Conserved quantities at edges
1898            *neighbours;
1899
1900        double beta_w; //Safety factor
1901        int N, err;
1902
1903
1904        // Convert Python arguments to C
1905        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1906          PyErr_SetString(PyExc_RuntimeError, 
1907                          "quantity_ext.c: limit_edges_by_neighbour could not parse input");
1908          return NULL;
1909        }
1910
1911        domain = PyObject_GetAttrString(quantity, "domain");
1912        if (!domain) {
1913          PyErr_SetString(PyExc_RuntimeError, 
1914                          "quantity_ext.c: limit_edges_by_neighbour could not obtain domain object from quantity");                     
1915         
1916          return NULL;
1917        }
1918
1919        // Get safety factor beta_w
1920        Tmp = PyObject_GetAttrString(domain, "beta_w");
1921        if (!Tmp) {
1922          PyErr_SetString(PyExc_RuntimeError, 
1923                          "quantity_ext.c: limit_by_vertex could not obtain beta_w object from domain");                       
1924         
1925          return NULL;
1926        }       
1927
1928
1929        // Get pertinent variables
1930        neighbours       = get_consecutive_array(domain, "neighbours");
1931        centroid_values  = get_consecutive_array(quantity, "centroid_values");
1932        vertex_values    = get_consecutive_array(quantity, "vertex_values");
1933        edge_values      = get_consecutive_array(quantity, "edge_values");
1934        beta_w           = PyFloat_AsDouble(Tmp);
1935
1936
1937        N = centroid_values -> dimensions[0];
1938
1939        err = _limit_edges_by_neighbour(N, beta_w,
1940                                        (double*) centroid_values -> data,
1941                                        (double*) vertex_values -> data,
1942                                        (double*) edge_values -> data,
1943                                        (long*)   neighbours -> data);
1944       
1945        if (err != 0) {
1946          PyErr_SetString(PyExc_RuntimeError,
1947                          "Internal function _limit_by_vertex failed");
1948          return NULL;
1949        }       
1950
1951
1952        // Release
1953        Py_DECREF(domain);
1954        Py_DECREF(neighbours);
1955        Py_DECREF(centroid_values);
1956        Py_DECREF(vertex_values);
1957        Py_DECREF(edge_values);
1958        Py_DECREF(Tmp);
1959
1960
1961        return Py_BuildValue("");
1962}
1963
1964
1965PyObject *limit_gradient_by_neighbour(PyObject *self, PyObject *args) {
1966  //Limit slopes for each volume to eliminate artificial variance
1967  //introduced by e.g. second order extrapolator
1968
1969  //This is an unsophisticated limiter as it does not take into
1970  //account dependencies among quantities.
1971
1972  //precondition:
1973  //  vertex values are estimated from gradient
1974  //postcondition:
1975  //  vertex and edge values are updated
1976  //
1977
1978        PyObject *quantity, *domain, *Tmp;
1979        PyArrayObject
1980            *vertex_values,   //Conserved quantities at vertices
1981            *centroid_values, //Conserved quantities at centroids
1982            *edge_values,     //Conserved quantities at edges
1983            *x_gradient,
1984            *y_gradient,
1985            *neighbours;
1986
1987        double beta_w; //Safety factor
1988        int N, err;
1989
1990
1991        // Convert Python arguments to C
1992        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1993          PyErr_SetString(PyExc_RuntimeError, 
1994                          "quantity_ext.c: limit_gradient_by_neighbour could not parse input");
1995          return NULL;
1996        }
1997
1998        domain = PyObject_GetAttrString(quantity, "domain");
1999        if (!domain) {
2000          PyErr_SetString(PyExc_RuntimeError, 
2001                          "quantity_ext.c: limit_gradient_by_neighbour could not obtain domain object from quantity");                 
2002         
2003          return NULL;
2004        }
2005
2006        // Get safety factor beta_w
2007        Tmp = PyObject_GetAttrString(domain, "beta_w");
2008        if (!Tmp) {
2009          PyErr_SetString(PyExc_RuntimeError, 
2010                          "quantity_ext.c: limit_gradient_by_neighbour could not obtain beta_w object from domain");                   
2011         
2012          return NULL;
2013        }       
2014
2015
2016        // Get pertinent variables
2017        neighbours       = get_consecutive_array(domain, "neighbours");
2018        centroid_values  = get_consecutive_array(quantity, "centroid_values");
2019        vertex_values    = get_consecutive_array(quantity, "vertex_values");
2020        edge_values      = get_consecutive_array(quantity, "edge_values");
2021        x_gradient       = get_consecutive_array(quantity, "x_gradient");
2022        y_gradient       = get_consecutive_array(quantity, "y_gradient");
2023
2024        beta_w           = PyFloat_AsDouble(Tmp);
2025
2026
2027        N = centroid_values -> dimensions[0];
2028
2029        err = _limit_gradient_by_neighbour(N, beta_w,
2030                                        (double*) centroid_values -> data,
2031                                        (double*) vertex_values -> data,
2032                                        (double*) edge_values -> data,
2033                                        (double*) x_gradient -> data,
2034                                        (double*) y_gradient -> data,
2035                                        (long*)   neighbours -> data);
2036       
2037        if (err != 0) {
2038          PyErr_SetString(PyExc_RuntimeError,
2039                          "Internal function _limit_gradient_by_neighbour failed");
2040          return NULL;
2041        }       
2042
2043
2044        // Release
2045        Py_DECREF(neighbours);
2046        Py_DECREF(centroid_values);
2047        Py_DECREF(vertex_values);
2048        Py_DECREF(edge_values);
2049        Py_DECREF(x_gradient);
2050        Py_DECREF(y_gradient);
2051        Py_DECREF(Tmp);
2052
2053
2054        return Py_BuildValue("");
2055}
2056
2057
2058// Method table for python module
2059static struct PyMethodDef MethodTable[] = {
2060        {"limit_old", limit_old, METH_VARARGS, "Print out"},
2061        {"limit_vertices_by_all_neighbours", limit_vertices_by_all_neighbours, METH_VARARGS, "Print out"},
2062        {"limit_edges_by_all_neighbours", limit_edges_by_all_neighbours, METH_VARARGS, "Print out"},
2063        {"limit_edges_by_neighbour", limit_edges_by_neighbour, METH_VARARGS, "Print out"},
2064        {"limit_gradient_by_neighbour", limit_gradient_by_neighbour, METH_VARARGS, "Print out"},
2065        {"bound_vertices_below_by_constant", bound_vertices_below_by_constant, METH_VARARGS, "Print out"},
2066        {"bound_vertices_below_by_quantity", bound_vertices_below_by_quantity, METH_VARARGS, "Print out"},
2067        {"update", update, METH_VARARGS, "Print out"},
2068        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
2069        {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"},
2070        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
2071        {"extrapolate_from_gradient", extrapolate_from_gradient,
2072                METH_VARARGS, "Print out"},
2073        {"extrapolate_second_order_and_limit_by_edge", extrapolate_second_order_and_limit_by_edge,
2074                METH_VARARGS, "Print out"},
2075        {"extrapolate_second_order_and_limit_by_vertex", extrapolate_second_order_and_limit_by_vertex,
2076                METH_VARARGS, "Print out"},
2077        {"interpolate_from_vertices_to_edges",
2078                interpolate_from_vertices_to_edges,
2079                METH_VARARGS, "Print out"},
2080        {"interpolate_from_edges_to_vertices",
2081                interpolate_from_edges_to_vertices,
2082                METH_VARARGS, "Print out"},
2083        {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"},           
2084        {NULL, NULL, 0, NULL}   // sentinel
2085};
2086
2087// Module initialisation
2088void initquantity_ext(void){
2089  Py_InitModule("quantity_ext", MethodTable);
2090
2091  import_array(); // Necessary for handling of NumPY structures
2092}
Note: See TracBrowser for help on using the repository browser.