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

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

Updating to relocated repository

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