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

Last change on this file since 6654 was 6654, checked in by ole, 16 years ago

Added volumetric balance report for boundary flows and total volume.
Flows due to forcing terms have not yet been added.
File size: 55.8 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
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       
1032        N = vertex_value_indices -> dimensions[0];
1033        // printf("Got parameters, N=%d\n", N);
1034        err = _average_vertex_values(N,
1035                                     (long*) vertex_value_indices -> data,
1036                                     (long*) number_of_triangles_per_node -> data,
1037                                     (double*) vertex_values -> data, 
1038                                     (double*) A -> data);
1039
1040        //printf("Error %d", err);
1041        if (err != 0) {
1042          PyErr_SetString(PyExc_RuntimeError, 
1043                          "average_vertex_values could not be computed");
1044          return NULL;
1045        }
1046
1047        return Py_BuildValue("");
1048}
1049
1050
1051
1052PyObject *extrapolate_from_gradient(PyObject *self, PyObject *args) {
1053
1054        PyObject *quantity, *domain;
1055        PyArrayObject
1056            *centroids,            //Coordinates at centroids
1057            *centroid_values,      //Values at centroids
1058            *vertex_coordinates,   //Coordinates at vertices
1059            *vertex_values,        //Values at vertices
1060            *edge_values,          //Values at edges
1061            *number_of_boundaries, //Number of boundaries for each triangle
1062            *surrogate_neighbours, //True neighbours or - if one missing - self
1063            *x_gradient,           //x gradient
1064            *y_gradient;           //y gradient
1065
1066        //int N, err;
1067        //int dimensions[1];
1068        int N, err;
1069        //double *a, *b;  //Gradients
1070
1071        // Convert Python arguments to C
1072        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1073          PyErr_SetString(PyExc_RuntimeError, 
1074                          "extrapolate_gradient could not parse input");       
1075          return NULL;
1076        }
1077
1078        domain = PyObject_GetAttrString(quantity, "domain");
1079        if (!domain) {
1080          PyErr_SetString(PyExc_RuntimeError, 
1081                          "extrapolate_gradient could not obtain domain object from quantity"); 
1082          return NULL;
1083        }
1084
1085        // Get pertinent variables
1086        centroids            = get_consecutive_array(domain,   "centroid_coordinates");
1087        centroid_values      = get_consecutive_array(quantity, "centroid_values");
1088        surrogate_neighbours = get_consecutive_array(domain,   "surrogate_neighbours");
1089        number_of_boundaries = get_consecutive_array(domain,   "number_of_boundaries");
1090        vertex_coordinates   = get_consecutive_array(domain,   "vertex_coordinates");
1091        vertex_values        = get_consecutive_array(quantity, "vertex_values");
1092        edge_values          = get_consecutive_array(quantity, "edge_values");
1093        x_gradient           = get_consecutive_array(quantity, "x_gradient");
1094        y_gradient           = get_consecutive_array(quantity, "y_gradient");
1095
1096        N = centroid_values -> dimensions[0];
1097
1098        // Release
1099        Py_DECREF(domain);
1100
1101        err = _extrapolate_from_gradient(N,
1102                        (double*) centroids -> data,
1103                        (double*) centroid_values -> data,
1104                        (double*) vertex_coordinates -> data,
1105                        (double*) vertex_values -> data,
1106                        (double*) edge_values -> data,
1107                        (double*) x_gradient -> data,
1108                        (double*) y_gradient -> data);
1109
1110
1111        if (err != 0) {
1112          PyErr_SetString(PyExc_RuntimeError,
1113                          "Internal function _extrapolate failed");
1114          return NULL;
1115        }
1116
1117
1118
1119        // Release
1120        Py_DECREF(centroids);
1121        Py_DECREF(centroid_values);
1122        Py_DECREF(number_of_boundaries);
1123        Py_DECREF(surrogate_neighbours);
1124        Py_DECREF(vertex_coordinates);
1125        Py_DECREF(vertex_values);
1126        Py_DECREF(edge_values);
1127        Py_DECREF(x_gradient);
1128        Py_DECREF(y_gradient);
1129
1130        return Py_BuildValue("");
1131}
1132
1133
1134PyObject *extrapolate_second_order_and_limit_by_edge(PyObject *self, PyObject *args) {
1135    /* Compute edge values using second order approximation and limit values
1136       so that edge values are limited by the two corresponding centroid values
1137       
1138       Python Call:
1139       extrapolate_second_order_and_limit(domain,quantity,beta)
1140    */
1141
1142    PyObject *quantity, *domain;
1143       
1144    PyArrayObject
1145        *domain_centroids,              //Coordinates at centroids
1146        *domain_vertex_coordinates,     //Coordinates at vertices
1147        *domain_number_of_boundaries,   //Number of boundaries for each triangle
1148        *domain_surrogate_neighbours,   //True neighbours or - if one missing - self
1149        *domain_neighbours,             //True neighbours, or if negative a link to boundary
1150
1151        *quantity_centroid_values,      //Values at centroids   
1152        *quantity_vertex_values,        //Values at vertices
1153        *quantity_edge_values,          //Values at edges
1154        *quantity_phi,                  //limiter phi values
1155        *quantity_x_gradient,           //x gradient
1156        *quantity_y_gradient;           //y gradient
1157   
1158
1159  // Local variables
1160  int ntri;
1161  double beta;
1162  int err;
1163   
1164  // Convert Python arguments to C
1165  if (!PyArg_ParseTuple(args, "O",&quantity)) {
1166      PyErr_SetString(PyExc_RuntimeError, 
1167                      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");     
1168      return NULL;
1169  }
1170
1171  domain = PyObject_GetAttrString(quantity, "domain");
1172  if (!domain) {
1173      PyErr_SetString(PyExc_RuntimeError, 
1174                      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");       
1175      return NULL;
1176  }
1177
1178       
1179  // Get pertinent variables
1180  domain_centroids              = get_consecutive_array(domain,   "centroid_coordinates");
1181  domain_surrogate_neighbours   = get_consecutive_array(domain,   "surrogate_neighbours");
1182  domain_number_of_boundaries   = get_consecutive_array(domain,   "number_of_boundaries");
1183  domain_vertex_coordinates     = get_consecutive_array(domain,   "vertex_coordinates");
1184  domain_neighbours             = get_consecutive_array(domain,   "neighbours");
1185
1186  quantity_centroid_values      = get_consecutive_array(quantity, "centroid_values");
1187  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
1188  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
1189  quantity_phi                  = get_consecutive_array(quantity, "phi");
1190  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
1191  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
1192
1193  beta = get_python_double(quantity,"beta");
1194
1195  ntri = quantity_centroid_values -> dimensions[0];
1196
1197  err = _compute_gradients(ntri,
1198                           (double*) domain_centroids -> data,
1199                           (double*) quantity_centroid_values -> data,
1200                           (long*)   domain_number_of_boundaries -> data,
1201                           (long*)   domain_surrogate_neighbours -> data,
1202                           (double*) quantity_x_gradient -> data,
1203                           (double*) quantity_y_gradient -> data);
1204
1205  if (err != 0) {
1206      PyErr_SetString(PyExc_RuntimeError,
1207                          "quantity_ext.c: Internal function _compute_gradient failed");
1208      return NULL;
1209  }
1210
1211
1212  err = _extrapolate_from_gradient(ntri, 
1213                                   (double*) domain_centroids -> data,
1214                                   (double*) quantity_centroid_values -> data,
1215                                   (double*) domain_vertex_coordinates -> data,
1216                                   (double*) quantity_vertex_values -> data,
1217                                   (double*) quantity_edge_values -> data,
1218                                   (double*) quantity_x_gradient -> data,
1219                                   (double*) quantity_y_gradient -> data);
1220
1221  if (err != 0) {
1222      PyErr_SetString(PyExc_RuntimeError,
1223                          "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
1224      return NULL;
1225  }
1226
1227
1228  err = _limit_edges_by_all_neighbours(ntri, beta,
1229                                       (double*) quantity_centroid_values -> data,
1230                                       (double*) quantity_vertex_values -> data,
1231                                       (double*) quantity_edge_values -> data,
1232                                       (long*)   domain_neighbours -> data,
1233                                       (double*) quantity_x_gradient -> data,
1234                                       (double*) quantity_y_gradient -> data);
1235
1236  if (err != 0) {
1237      PyErr_SetString(PyExc_RuntimeError,
1238                          "quantity_ext.c: Internal function _limit_edges_by_all_neighbours failed");
1239      return NULL;
1240  }
1241
1242
1243  // Release
1244  Py_DECREF(domain_centroids);
1245  Py_DECREF(domain_surrogate_neighbours);
1246  Py_DECREF(domain_number_of_boundaries);
1247  Py_DECREF(domain_vertex_coordinates);
1248
1249  Py_DECREF(quantity_centroid_values);
1250  Py_DECREF(quantity_vertex_values);
1251  Py_DECREF(quantity_edge_values);
1252  Py_DECREF(quantity_phi);
1253  Py_DECREF(quantity_x_gradient);
1254  Py_DECREF(quantity_y_gradient);
1255
1256  return Py_BuildValue("");
1257}
1258
1259
1260PyObject *extrapolate_second_order_and_limit_by_vertex(PyObject *self, PyObject *args) {
1261    /* Compute edge values using second order approximation and limit values
1262       so that edge values are limited by the two corresponding centroid values
1263       
1264       Python Call:
1265       extrapolate_second_order_and_limit(domain,quantity,beta)
1266    */
1267
1268    PyObject *quantity, *domain;
1269       
1270    PyArrayObject
1271        *domain_centroids,              //Coordinates at centroids
1272        *domain_vertex_coordinates,     //Coordinates at vertices
1273        *domain_number_of_boundaries,   //Number of boundaries for each triangle
1274        *domain_surrogate_neighbours,   //True neighbours or - if one missing - self
1275        *domain_neighbours,             //True neighbours, or if negative a link to boundary
1276
1277        *quantity_centroid_values,      //Values at centroids   
1278        *quantity_vertex_values,        //Values at vertices
1279        *quantity_edge_values,          //Values at edges
1280        *quantity_phi,                  //limiter phi values
1281        *quantity_x_gradient,           //x gradient
1282        *quantity_y_gradient;           //y gradient
1283   
1284
1285  // Local variables
1286  int ntri;
1287  double beta;
1288  int err;
1289   
1290  // Convert Python arguments to C
1291  if (!PyArg_ParseTuple(args, "O",&quantity)) {
1292      PyErr_SetString(PyExc_RuntimeError, 
1293                      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");     
1294      return NULL;
1295  }
1296
1297  domain = PyObject_GetAttrString(quantity, "domain");
1298  if (!domain) {
1299      PyErr_SetString(PyExc_RuntimeError, 
1300                      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");       
1301      return NULL;
1302  }
1303
1304       
1305  // Get pertinent variables
1306  domain_centroids              = get_consecutive_array(domain,   "centroid_coordinates");
1307  domain_surrogate_neighbours   = get_consecutive_array(domain,   "surrogate_neighbours");
1308  domain_number_of_boundaries   = get_consecutive_array(domain,   "number_of_boundaries");
1309  domain_vertex_coordinates     = get_consecutive_array(domain,   "vertex_coordinates");
1310  domain_neighbours             = get_consecutive_array(domain,   "neighbours");
1311
1312  quantity_centroid_values      = get_consecutive_array(quantity, "centroid_values");
1313  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
1314  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
1315  quantity_phi                  = get_consecutive_array(quantity, "phi");
1316  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
1317  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
1318
1319  beta = get_python_double(quantity,"beta");
1320
1321  ntri = quantity_centroid_values -> dimensions[0];
1322
1323  err = _compute_gradients(ntri,
1324                           (double*) domain_centroids -> data,
1325                           (double*) quantity_centroid_values -> data,
1326                           (long*)   domain_number_of_boundaries -> data,
1327                           (long*)   domain_surrogate_neighbours -> data,
1328                           (double*) quantity_x_gradient -> data,
1329                           (double*) quantity_y_gradient -> data);
1330
1331  if (err != 0) {
1332      PyErr_SetString(PyExc_RuntimeError,
1333                          "quantity_ext.c: Internal function _compute_gradient failed");
1334      return NULL;
1335  }
1336
1337
1338  err = _extrapolate_from_gradient(ntri, 
1339                                   (double*) domain_centroids -> data,
1340                                   (double*) quantity_centroid_values -> data,
1341                                   (double*) domain_vertex_coordinates -> data,
1342                                   (double*) quantity_vertex_values -> data,
1343                                   (double*) quantity_edge_values -> data,
1344                                   (double*) quantity_x_gradient -> data,
1345                                   (double*) quantity_y_gradient -> data);
1346
1347  if (err != 0) {
1348      PyErr_SetString(PyExc_RuntimeError,
1349                          "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
1350      return NULL;
1351  }
1352
1353
1354  err = _limit_vertices_by_all_neighbours(ntri, beta,
1355                                       (double*) quantity_centroid_values -> data,
1356                                       (double*) quantity_vertex_values -> data,
1357                                       (double*) quantity_edge_values -> data,
1358                                       (long*)   domain_neighbours -> data,
1359                                       (double*) quantity_x_gradient -> data,
1360                                       (double*) quantity_y_gradient -> data);
1361
1362  if (err != 0) {
1363      PyErr_SetString(PyExc_RuntimeError,
1364                          "quantity_ext.c: Internal function _limit_vertices_by_all_neighbours failed");
1365      return NULL;
1366  }
1367
1368
1369  // Release
1370  Py_DECREF(domain_centroids);
1371  Py_DECREF(domain_surrogate_neighbours);
1372  Py_DECREF(domain_number_of_boundaries);
1373  Py_DECREF(domain_vertex_coordinates);
1374
1375  Py_DECREF(quantity_centroid_values);
1376  Py_DECREF(quantity_vertex_values);
1377  Py_DECREF(quantity_edge_values);
1378  Py_DECREF(quantity_phi);
1379  Py_DECREF(quantity_x_gradient);
1380  Py_DECREF(quantity_y_gradient);
1381
1382  return Py_BuildValue("");
1383}
1384
1385
1386
1387PyObject *compute_gradients(PyObject *self, PyObject *args) {
1388
1389        PyObject *quantity, *domain;
1390        PyArrayObject
1391            *centroids,            //Coordinates at centroids
1392            *centroid_values,      //Values at centroids
1393            *vertex_coordinates,   //Coordinates at vertices
1394            *vertex_values,        //Values at vertices
1395            *edge_values,          //Values at edges
1396            *number_of_boundaries, //Number of boundaries for each triangle
1397            *surrogate_neighbours, //True neighbours or - if one missing - self
1398            *x_gradient,           //x gradient
1399            *y_gradient;           //y gradient
1400
1401        //int N, err;
1402        //int dimensions[1];
1403        int N, err;
1404        //double *a, *b;  //Gradients
1405
1406        // Convert Python arguments to C
1407        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1408          PyErr_SetString(PyExc_RuntimeError, 
1409                          "compute_gradients could not parse input");   
1410          return NULL;
1411        }
1412
1413        domain = PyObject_GetAttrString(quantity, "domain");
1414        if (!domain) {
1415          PyErr_SetString(PyExc_RuntimeError, 
1416                          "compute_gradients could not obtain domain object from quantity");   
1417          return NULL;
1418        }
1419
1420        // Get pertinent variables
1421        centroids            = get_consecutive_array(domain,   "centroid_coordinates");
1422        centroid_values      = get_consecutive_array(quantity, "centroid_values");
1423        surrogate_neighbours = get_consecutive_array(domain,   "surrogate_neighbours");
1424        number_of_boundaries = get_consecutive_array(domain,   "number_of_boundaries");
1425        vertex_coordinates   = get_consecutive_array(domain,   "vertex_coordinates");
1426        vertex_values        = get_consecutive_array(quantity, "vertex_values");
1427        edge_values          = get_consecutive_array(quantity, "edge_values");
1428        x_gradient           = get_consecutive_array(quantity, "x_gradient");
1429        y_gradient           = get_consecutive_array(quantity, "y_gradient");
1430
1431        N = centroid_values -> dimensions[0];
1432
1433        // Release
1434        Py_DECREF(domain);
1435
1436
1437        err = _compute_gradients(N,
1438                        (double*) centroids -> data,
1439                        (double*) centroid_values -> data,
1440                        (long*) number_of_boundaries -> data,
1441                        (long*) surrogate_neighbours -> data,
1442                        (double*) x_gradient -> data,
1443                        (double*) y_gradient -> data);
1444
1445        if (err != 0) {
1446          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
1447          return NULL;
1448        }
1449
1450
1451
1452        // Release
1453        Py_DECREF(centroids);
1454        Py_DECREF(centroid_values);
1455        Py_DECREF(number_of_boundaries);
1456        Py_DECREF(surrogate_neighbours);
1457        Py_DECREF(vertex_coordinates);
1458        Py_DECREF(vertex_values);
1459        Py_DECREF(edge_values);
1460        Py_DECREF(x_gradient);
1461        Py_DECREF(y_gradient);
1462
1463        return Py_BuildValue("");
1464}
1465
1466
1467
1468PyObject *limit_old(PyObject *self, PyObject *args) {
1469  //Limit slopes for each volume to eliminate artificial variance
1470  //introduced by e.g. second order extrapolator
1471
1472  //This is an unsophisticated limiter as it does not take into
1473  //account dependencies among quantities.
1474
1475  //precondition:
1476  //  vertex values are estimated from gradient
1477  //postcondition:
1478  //  vertex values are updated
1479  //
1480
1481        PyObject *quantity, *domain, *Tmp;
1482        PyArrayObject
1483            *qv, //Conserved quantities at vertices
1484            *qc, //Conserved quantities at centroids
1485            *neighbours;
1486
1487        int k, i, n, N, k3;
1488        double beta_w; //Safety factor
1489        double *qmin, *qmax, qn;
1490
1491        // Convert Python arguments to C
1492        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1493          PyErr_SetString(PyExc_RuntimeError, 
1494                          "quantity_ext.c: limit_old could not parse input");
1495          return NULL;
1496        }
1497
1498        domain = PyObject_GetAttrString(quantity, "domain");
1499        if (!domain) {
1500          PyErr_SetString(PyExc_RuntimeError, 
1501                          "quantity_ext.c: limit_old could not obtain domain object from quantity");                   
1502         
1503          return NULL;
1504        }
1505
1506        //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
1507        neighbours = get_consecutive_array(domain, "neighbours");
1508
1509        // Get safety factor beta_w
1510        Tmp = PyObject_GetAttrString(domain, "beta_w");
1511        if (!Tmp) {
1512          PyErr_SetString(PyExc_RuntimeError, 
1513                          "quantity_ext.c: limit_old could not obtain beta_w object from domain");                     
1514         
1515          return NULL;
1516        }       
1517
1518        beta_w = PyFloat_AsDouble(Tmp);
1519
1520        Py_DECREF(Tmp);
1521        Py_DECREF(domain);
1522
1523
1524        qc = get_consecutive_array(quantity, "centroid_values");
1525        qv = get_consecutive_array(quantity, "vertex_values");
1526
1527
1528        N = qc -> dimensions[0];
1529
1530        // Find min and max of this and neighbour's centroid values
1531        qmin = malloc(N * sizeof(double));
1532        qmax = malloc(N * sizeof(double));
1533        for (k=0; k<N; k++) {
1534                k3=k*3;
1535
1536                qmin[k] = ((double*) qc -> data)[k];
1537                qmax[k] = qmin[k];
1538
1539                for (i=0; i<3; i++) {
1540                        n = ((long*) neighbours -> data)[k3+i];
1541                        if (n >= 0) {
1542                                qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
1543
1544                                qmin[k] = min(qmin[k], qn);
1545                                qmax[k] = max(qmax[k], qn);
1546                        }
1547                        //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]);
1548                        //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]);
1549                }
1550        }
1551
1552        // Call underlying routine
1553        _limit_old(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
1554
1555        free(qmin);
1556        free(qmax);
1557        return Py_BuildValue("");
1558}
1559
1560
1561PyObject *limit_vertices_by_all_neighbours(PyObject *self, PyObject *args) {
1562  //Limit slopes for each volume to eliminate artificial variance
1563  //introduced by e.g. second order extrapolator
1564
1565  //This is an unsophisticated limiter as it does not take into
1566  //account dependencies among quantities.
1567
1568  //precondition:
1569  //  vertex values are estimated from gradient
1570  //postcondition:
1571  //  vertex and edge values are updated
1572  //
1573
1574        PyObject *quantity, *domain, *Tmp;
1575        PyArrayObject
1576            *vertex_values,   //Conserved quantities at vertices
1577            *centroid_values, //Conserved quantities at centroids
1578            *edge_values,     //Conserved quantities at edges
1579            *neighbours,
1580            *x_gradient,
1581            *y_gradient;
1582
1583        double beta_w; //Safety factor
1584        int N, err;
1585
1586
1587        // Convert Python arguments to C
1588        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1589          PyErr_SetString(PyExc_RuntimeError, 
1590                          "quantity_ext.c: limit_by_vertex could not parse input");
1591          return NULL;
1592        }
1593
1594        domain = PyObject_GetAttrString(quantity, "domain");
1595        if (!domain) {
1596          PyErr_SetString(PyExc_RuntimeError, 
1597                          "quantity_ext.c: limit_by_vertex could not obtain domain object from quantity");                     
1598         
1599          return NULL;
1600        }
1601
1602        // Get safety factor beta_w
1603        Tmp = PyObject_GetAttrString(domain, "beta_w");
1604        if (!Tmp) {
1605          PyErr_SetString(PyExc_RuntimeError, 
1606                          "quantity_ext.c: limit_by_vertex could not obtain beta_w object from domain");                       
1607         
1608          return NULL;
1609        }       
1610
1611
1612        // Get pertinent variables
1613        neighbours       = get_consecutive_array(domain, "neighbours");
1614        centroid_values  = get_consecutive_array(quantity, "centroid_values");
1615        vertex_values    = get_consecutive_array(quantity, "vertex_values");
1616        edge_values      = get_consecutive_array(quantity, "edge_values");
1617        x_gradient       = get_consecutive_array(quantity, "x_gradient");
1618        y_gradient       = get_consecutive_array(quantity, "y_gradient");
1619        beta_w           = get_python_double(domain,"beta_w");
1620
1621
1622
1623        N = centroid_values -> dimensions[0];
1624
1625        err = _limit_vertices_by_all_neighbours(N, beta_w,
1626                                                (double*) centroid_values -> data,
1627                                                (double*) vertex_values -> data,
1628                                                (double*) edge_values -> data,
1629                                                (long*)   neighbours -> data,
1630                                                (double*) x_gradient -> data,
1631                                                (double*) y_gradient -> data);
1632
1633
1634       
1635        if (err != 0) {
1636          PyErr_SetString(PyExc_RuntimeError,
1637                          "Internal function _limit_by_vertex failed");
1638          return NULL;
1639        }       
1640
1641
1642        // Release
1643        Py_DECREF(neighbours);
1644        Py_DECREF(centroid_values);
1645        Py_DECREF(vertex_values);
1646        Py_DECREF(edge_values);
1647        Py_DECREF(x_gradient);
1648        Py_DECREF(y_gradient);
1649        Py_DECREF(Tmp);
1650
1651
1652        return Py_BuildValue("");
1653}
1654
1655
1656
1657PyObject *limit_edges_by_all_neighbours(PyObject *self, PyObject *args) {
1658  //Limit slopes for each volume to eliminate artificial variance
1659  //introduced by e.g. second order extrapolator
1660
1661  //This is an unsophisticated limiter as it does not take into
1662  //account dependencies among quantities.
1663
1664  //precondition:
1665  //  vertex values are estimated from gradient
1666  //postcondition:
1667  //  vertex and edge values are updated
1668  //
1669
1670        PyObject *quantity, *domain;
1671        PyArrayObject
1672            *vertex_values,   //Conserved quantities at vertices
1673            *centroid_values, //Conserved quantities at centroids
1674            *edge_values,     //Conserved quantities at edges
1675            *x_gradient,
1676            *y_gradient,
1677            *neighbours;
1678
1679        double beta_w; //Safety factor
1680        int N, err;
1681
1682
1683        // Convert Python arguments to C
1684        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1685          PyErr_SetString(PyExc_RuntimeError, 
1686                          "quantity_ext.c: limit_edges_by_all_neighbours could not parse input");
1687          return NULL;
1688        }
1689
1690        domain = get_python_object(quantity, "domain");
1691
1692
1693        // Get pertinent variables
1694        neighbours       = get_consecutive_array(domain, "neighbours");
1695        centroid_values  = get_consecutive_array(quantity, "centroid_values");
1696        vertex_values    = get_consecutive_array(quantity, "vertex_values");
1697        edge_values      = get_consecutive_array(quantity, "edge_values");
1698        x_gradient       = get_consecutive_array(quantity, "x_gradient");
1699        y_gradient       = get_consecutive_array(quantity, "y_gradient");
1700        beta_w           = get_python_double(domain,"beta_w");
1701
1702
1703
1704        N = centroid_values -> dimensions[0];
1705
1706        err = _limit_edges_by_all_neighbours(N, beta_w,
1707                                             (double*) centroid_values -> data,
1708                                             (double*) vertex_values -> data,
1709                                             (double*) edge_values -> data,
1710                                             (long*)   neighbours -> data,
1711                                             (double*) x_gradient -> data,
1712                                             (double*) y_gradient -> data);
1713
1714        if (err != 0) {
1715          PyErr_SetString(PyExc_RuntimeError,
1716                          "quantity_ect.c: limit_edges_by_all_neighbours internal function _limit_edges_by_all_neighbours failed");
1717          return NULL;
1718        }       
1719
1720
1721        // Release
1722        Py_DECREF(neighbours);
1723        Py_DECREF(centroid_values);
1724        Py_DECREF(vertex_values);
1725        Py_DECREF(edge_values);
1726        Py_DECREF(x_gradient);
1727        Py_DECREF(y_gradient);
1728
1729
1730
1731        return Py_BuildValue("");
1732}
1733
1734PyObject *bound_vertices_below_by_constant(PyObject *self, PyObject *args) {
1735  //Bound a quantity below by a contant (useful for ensuring positivity
1736  //precondition:
1737  //  vertex values are already calulated, gradient consistent
1738  //postcondition:
1739  //  gradient, vertex and edge values are updated
1740  //
1741
1742        PyObject *quantity, *domain;
1743        PyArrayObject
1744            *vertex_values,   //Conserved quantities at vertices
1745            *centroid_values, //Conserved quantities at centroids
1746            *edge_values,     //Conserved quantities at edges
1747            *x_gradient,
1748            *y_gradient;
1749
1750        double bound; //Safety factor
1751        int N, err;
1752
1753
1754        // Convert Python arguments to C
1755        if (!PyArg_ParseTuple(args, "Od", &quantity, &bound)) {
1756          PyErr_SetString(PyExc_RuntimeError, 
1757                          "quantity_ext.c: bound_vertices_below_by_constant could not parse input");
1758          return NULL;
1759        }
1760
1761        domain = get_python_object(quantity, "domain");
1762
1763        // Get pertinent variables
1764        centroid_values  = get_consecutive_array(quantity, "centroid_values");
1765        vertex_values    = get_consecutive_array(quantity, "vertex_values");
1766        edge_values      = get_consecutive_array(quantity, "edge_values");
1767        x_gradient       = get_consecutive_array(quantity, "x_gradient");
1768        y_gradient       = get_consecutive_array(quantity, "y_gradient");
1769
1770
1771
1772
1773        N = centroid_values -> dimensions[0];
1774
1775        err = _bound_vertices_below_by_constant(N, bound,
1776                                             (double*) centroid_values -> data,
1777                                             (double*) vertex_values -> data,
1778                                             (double*) edge_values -> data,
1779                                             (double*) x_gradient -> data,
1780                                             (double*) y_gradient -> data);
1781
1782        if (err != 0) {
1783          PyErr_SetString(PyExc_RuntimeError,
1784                          "quantity_ect.c: bound_vertices_below_by_constant internal function _bound_vertices_below_by_constant failed");
1785          return NULL;
1786        }       
1787
1788
1789        // Release
1790        Py_DECREF(centroid_values);
1791        Py_DECREF(vertex_values);
1792        Py_DECREF(edge_values);
1793        Py_DECREF(x_gradient);
1794        Py_DECREF(y_gradient);
1795
1796
1797
1798        return Py_BuildValue("");
1799}
1800
1801
1802PyObject *bound_vertices_below_by_quantity(PyObject *self, PyObject *args) {
1803  //Bound a quantity below by a contant (useful for ensuring positivity
1804  //precondition:
1805  //  vertex values are already calulated, gradient consistent
1806  //postcondition:
1807  //  gradient, vertex and edge values are updated
1808  //
1809
1810        PyObject *quantity, *bounding_quantity, *domain;
1811        PyArrayObject
1812            *vertex_values,   //Conserved quantities at vertices
1813            *centroid_values, //Conserved quantities at centroids
1814            *edge_values,     //Conserved quantities at edges
1815            *x_gradient,
1816            *y_gradient,
1817            *bound_vertex_values;
1818
1819        int N, err;
1820
1821
1822        // Convert Python arguments to C
1823        if (!PyArg_ParseTuple(args, "OO", &quantity, &bounding_quantity)) {
1824          PyErr_SetString(PyExc_RuntimeError, 
1825                          "quantity_ext.c: bound_vertices_below_by_quantity could not parse input");
1826          return NULL;
1827        }
1828
1829        domain = get_python_object(quantity, "domain");
1830
1831        // Get pertinent variables
1832        centroid_values     = get_consecutive_array(quantity, "centroid_values");
1833        vertex_values       = get_consecutive_array(quantity, "vertex_values");
1834        edge_values         = get_consecutive_array(quantity, "edge_values");
1835        x_gradient          = get_consecutive_array(quantity, "x_gradient");
1836        y_gradient          = get_consecutive_array(quantity, "y_gradient");
1837        bound_vertex_values = get_consecutive_array(bounding_quantity, "vertex_values");
1838
1839
1840
1841        Py_DECREF(domain);
1842
1843        N = centroid_values -> dimensions[0];
1844
1845        err = _bound_vertices_below_by_quantity(N, 
1846                                                (double*) bound_vertex_values -> data,
1847                                                (double*) centroid_values -> data,
1848                                                (double*) vertex_values -> data,
1849                                                (double*) edge_values -> data,
1850                                                (double*) x_gradient -> data,
1851                                                (double*) y_gradient -> data);
1852
1853        if (err != 0) {
1854          PyErr_SetString(PyExc_RuntimeError,
1855                          "quantity_ect.c: bound_vertices_below_by_quantity internal function _bound_vertices_below_by_quantity failed");
1856          return NULL;
1857        }       
1858
1859
1860        // Release
1861        Py_DECREF(centroid_values);
1862        Py_DECREF(vertex_values);
1863        Py_DECREF(edge_values);
1864        Py_DECREF(x_gradient);
1865        Py_DECREF(y_gradient);
1866        Py_DECREF(bound_vertex_values);
1867
1868
1869
1870        return Py_BuildValue("");
1871}
1872
1873
1874PyObject *limit_edges_by_neighbour(PyObject *self, PyObject *args) {
1875  //Limit slopes for each volume to eliminate artificial variance
1876  //introduced by e.g. second order extrapolator
1877
1878  //This is an unsophisticated limiter as it does not take into
1879  //account dependencies among quantities.
1880
1881  //precondition:
1882  //  vertex values are estimated from gradient
1883  //postcondition:
1884  //  vertex and edge values are updated
1885  //
1886
1887        PyObject *quantity, *domain, *Tmp;
1888        PyArrayObject
1889            *vertex_values,   //Conserved quantities at vertices
1890            *centroid_values, //Conserved quantities at centroids
1891            *edge_values,     //Conserved quantities at edges
1892            *neighbours;
1893
1894        double beta_w; //Safety factor
1895        int N, err;
1896
1897
1898        // Convert Python arguments to C
1899        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1900          PyErr_SetString(PyExc_RuntimeError, 
1901                          "quantity_ext.c: limit_edges_by_neighbour could not parse input");
1902          return NULL;
1903        }
1904
1905        domain = PyObject_GetAttrString(quantity, "domain");
1906        if (!domain) {
1907          PyErr_SetString(PyExc_RuntimeError, 
1908                          "quantity_ext.c: limit_edges_by_neighbour could not obtain domain object from quantity");                     
1909         
1910          return NULL;
1911        }
1912
1913        // Get safety factor beta_w
1914        Tmp = PyObject_GetAttrString(domain, "beta_w");
1915        if (!Tmp) {
1916          PyErr_SetString(PyExc_RuntimeError, 
1917                          "quantity_ext.c: limit_by_vertex could not obtain beta_w object from domain");                       
1918         
1919          return NULL;
1920        }       
1921
1922
1923        // Get pertinent variables
1924        neighbours       = get_consecutive_array(domain, "neighbours");
1925        centroid_values  = get_consecutive_array(quantity, "centroid_values");
1926        vertex_values    = get_consecutive_array(quantity, "vertex_values");
1927        edge_values      = get_consecutive_array(quantity, "edge_values");
1928        beta_w           = PyFloat_AsDouble(Tmp);
1929
1930
1931        N = centroid_values -> dimensions[0];
1932
1933        err = _limit_edges_by_neighbour(N, beta_w,
1934                                        (double*) centroid_values -> data,
1935                                        (double*) vertex_values -> data,
1936                                        (double*) edge_values -> data,
1937                                        (long*)   neighbours -> data);
1938       
1939        if (err != 0) {
1940          PyErr_SetString(PyExc_RuntimeError,
1941                          "Internal function _limit_by_vertex failed");
1942          return NULL;
1943        }       
1944
1945
1946        // Release
1947        Py_DECREF(domain);
1948        Py_DECREF(neighbours);
1949        Py_DECREF(centroid_values);
1950        Py_DECREF(vertex_values);
1951        Py_DECREF(edge_values);
1952        Py_DECREF(Tmp);
1953
1954
1955        return Py_BuildValue("");
1956}
1957
1958
1959PyObject *limit_gradient_by_neighbour(PyObject *self, PyObject *args) {
1960  //Limit slopes for each volume to eliminate artificial variance
1961  //introduced by e.g. second order extrapolator
1962
1963  //This is an unsophisticated limiter as it does not take into
1964  //account dependencies among quantities.
1965
1966  //precondition:
1967  //  vertex values are estimated from gradient
1968  //postcondition:
1969  //  vertex and edge values are updated
1970  //
1971
1972        PyObject *quantity, *domain, *Tmp;
1973        PyArrayObject
1974            *vertex_values,   //Conserved quantities at vertices
1975            *centroid_values, //Conserved quantities at centroids
1976            *edge_values,     //Conserved quantities at edges
1977            *x_gradient,
1978            *y_gradient,
1979            *neighbours;
1980
1981        double beta_w; //Safety factor
1982        int N, err;
1983
1984
1985        // Convert Python arguments to C
1986        if (!PyArg_ParseTuple(args, "O", &quantity)) {
1987          PyErr_SetString(PyExc_RuntimeError, 
1988                          "quantity_ext.c: limit_gradient_by_neighbour could not parse input");
1989          return NULL;
1990        }
1991
1992        domain = PyObject_GetAttrString(quantity, "domain");
1993        if (!domain) {
1994          PyErr_SetString(PyExc_RuntimeError, 
1995                          "quantity_ext.c: limit_gradient_by_neighbour could not obtain domain object from quantity");                 
1996         
1997          return NULL;
1998        }
1999
2000        // Get safety factor beta_w
2001        Tmp = PyObject_GetAttrString(domain, "beta_w");
2002        if (!Tmp) {
2003          PyErr_SetString(PyExc_RuntimeError, 
2004                          "quantity_ext.c: limit_gradient_by_neighbour could not obtain beta_w object from domain");                   
2005         
2006          return NULL;
2007        }       
2008
2009
2010        // Get pertinent variables
2011        neighbours       = get_consecutive_array(domain, "neighbours");
2012        centroid_values  = get_consecutive_array(quantity, "centroid_values");
2013        vertex_values    = get_consecutive_array(quantity, "vertex_values");
2014        edge_values      = get_consecutive_array(quantity, "edge_values");
2015        x_gradient       = get_consecutive_array(quantity, "x_gradient");
2016        y_gradient       = get_consecutive_array(quantity, "y_gradient");
2017
2018        beta_w           = PyFloat_AsDouble(Tmp);
2019
2020
2021        N = centroid_values -> dimensions[0];
2022
2023        err = _limit_gradient_by_neighbour(N, beta_w,
2024                                        (double*) centroid_values -> data,
2025                                        (double*) vertex_values -> data,
2026                                        (double*) edge_values -> data,
2027                                        (double*) x_gradient -> data,
2028                                        (double*) y_gradient -> data,
2029                                        (long*)   neighbours -> data);
2030       
2031        if (err != 0) {
2032          PyErr_SetString(PyExc_RuntimeError,
2033                          "Internal function _limit_gradient_by_neighbour failed");
2034          return NULL;
2035        }       
2036
2037
2038        // Release
2039        Py_DECREF(neighbours);
2040        Py_DECREF(centroid_values);
2041        Py_DECREF(vertex_values);
2042        Py_DECREF(edge_values);
2043        Py_DECREF(x_gradient);
2044        Py_DECREF(y_gradient);
2045        Py_DECREF(Tmp);
2046
2047
2048        return Py_BuildValue("");
2049}
2050
2051
2052// Method table for python module
2053static struct PyMethodDef MethodTable[] = {
2054        {"limit_old", limit_old, METH_VARARGS, "Print out"},
2055        {"limit_vertices_by_all_neighbours", limit_vertices_by_all_neighbours, METH_VARARGS, "Print out"},
2056        {"limit_edges_by_all_neighbours", limit_edges_by_all_neighbours, METH_VARARGS, "Print out"},
2057        {"limit_edges_by_neighbour", limit_edges_by_neighbour, METH_VARARGS, "Print out"},
2058        {"limit_gradient_by_neighbour", limit_gradient_by_neighbour, METH_VARARGS, "Print out"},
2059        {"bound_vertices_below_by_constant", bound_vertices_below_by_constant, METH_VARARGS, "Print out"},
2060        {"bound_vertices_below_by_quantity", bound_vertices_below_by_quantity, METH_VARARGS, "Print out"},
2061        {"update", update, METH_VARARGS, "Print out"},
2062        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
2063        {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"},
2064        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
2065        {"extrapolate_from_gradient", extrapolate_from_gradient,
2066                METH_VARARGS, "Print out"},
2067        {"extrapolate_second_order_and_limit_by_edge", extrapolate_second_order_and_limit_by_edge,
2068                METH_VARARGS, "Print out"},
2069        {"extrapolate_second_order_and_limit_by_vertex", extrapolate_second_order_and_limit_by_vertex,
2070                METH_VARARGS, "Print out"},
2071        {"interpolate_from_vertices_to_edges",
2072                interpolate_from_vertices_to_edges,
2073                METH_VARARGS, "Print out"},
2074        {"interpolate_from_edges_to_vertices",
2075                interpolate_from_edges_to_vertices,
2076                METH_VARARGS, "Print out"},
2077        {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"},           
2078        {NULL, NULL, 0, NULL}   // sentinel
2079};
2080
2081// Module initialisation
2082void initquantity_ext(void){
2083  Py_InitModule("quantity_ext", MethodTable);
2084
2085  import_array(); // Necessary for handling of NumPY structures
2086}
Note: See TracBrowser for help on using the repository browser.