source: trunk/anuga_work/anuga_cuda/src/anuga_HMPP/quantity.c @ 9329

Last change on this file since 9329 was 9017, checked in by steve, 12 years ago

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