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

Last change on this file since 7573 was 7573, checked in by steve, 14 years ago

Committing a version of shallow_water_balanced which passes it unit test
using a version of edge limiting which doesn't limit boundary edges. THis
is useful to allow linear functions to be reconstructed.

Had to play with the transmissive BC and use centroid values which is
set via domain set_centroid_transmissive_bc

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