source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c @ 8123

Last change on this file since 8123 was 8123, checked in by steve, 13 years ago

Added in a procedure to calculate local gradients on each triangle using
vectex values

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