Changeset 4978

Timestamp:

Jan 28, 2008, 6:14:47 PM (17 years ago)

Author:

steve

Message:

Added a C extension for the advection directory (gave up on using f2py!)

Location:

anuga_core

Files:

• compile_all.py (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/quantity_ext.c (modified) (2 diffs)
• source/anuga/advection/advection.py (modified) (3 diffs)
• source/anuga/advection/advection_ext.c (modified) (6 diffs)
• source/anuga/advection/advection_ext.pyf (deleted)
• source/anuga/advection/compile_all.py (deleted)
• source/anuga/advection/setup.py (deleted)
• source/anuga/advection/test_advection.py (modified) (1 diff)
• source/anuga/advection/test_all.py (added)
• source/anuga/compile_all.py (modified) (1 diff)
• source/anuga/shallow_water/shallow_water_ext.c (modified) (2 diffs)
• source/anuga/test_all.py (modified) (1 diff)
• source/anuga/utilities/util_ext.h (modified) (2 diffs)

anuga_core/compile_all.py

@@ -30 +30 @@
 
 os.chdir('..')
+os.chdir('advection')
+execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
+
+
+os.chdir('..')
 os.chdir('shallow_water')
 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -594 +594 @@
         }
 
-        centroid_values = get_consecutive_array(quantity, "centroid_values");
-        explicit_update = get_consecutive_array(quantity, "explicit_update");
+        centroid_values      = get_consecutive_array(quantity, "centroid_values");
+        explicit_update      = get_consecutive_array(quantity, "explicit_update");
         semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");
 
@@ -608 +608 @@
         if (err != 0) {
           PyErr_SetString(PyExc_RuntimeError,
-                          "Zero division in semi implicit update - call Stephen :)");
+                          "quantity_ext.c: update, divsion by zero in semi implicit update - call Stephen :)");
           return NULL;
         }

anuga_core/source/anuga/advection/advection.py

@@ -66 +66 @@
                                 numproc=numproc)
 
-
+        import Numeric
         if velocity is None:
-            self.velocity = [1,0]
+            self.velocity = Numeric.array([1,0],'d')
         else:
-            self.velocity = velocity
+            self.velocity = Numeric.array(velocity,'d')
 
         #Only first is implemented for advection
@@ -108 +108 @@
         return flux, max_speed
 
+
+
     def compute_fluxes(self):
-
-        
-        self.compute_fluxes_ext()
-##         try:
-##             self.compute_fluxes_ext()
-##         except:
-##             self.compute_fluxes_python()
- 
+        """Compute all fluxes and the timestep suitable for all volumes
+        in domain.
+
+        Compute total flux for each conserved quantity using "flux_function"
+
+        Fluxes across each edge are scaled by edgelengths and summed up
+        Resulting flux is then scaled by area and stored in
+        domain.explicit_update
+
+        The maximal allowable speed computed by the flux_function
+        for each volume
+        is converted to a timestep that must not be exceeded. The minimum of
+        those is computed as the next overall timestep.
+
+        Post conditions:
+        domain.explicit_update is reset to computed flux values
+        domain.timestep is set to the largest step satisfying all volumes.
+        """
+
+        import sys
+        from Numeric import zeros, Float
+        from anuga.config import max_timestep
+
+
+        huge_timestep = float(sys.maxint)
+        Stage = self.quantities['stage']
+
+        """
+        print "======================================"
+        print "BEFORE compute_fluxes"
+        print "stage_update",Stage.explicit_update
+        print "stage_edge",Stage.edge_values
+        print "stage_bdry",Stage.boundary_values
+        print "neighbours",self.neighbours
+        print "neighbour_edges",self.neighbour_edges
+        print "normals",self.normals
+        print "areas",self.areas
+        print "radii",self.radii
+        print "edgelengths",self.edgelengths
+        print "tri_full_flag",self.tri_full_flag
+        print "huge_timestep",huge_timestep
+        print "max_timestep",max_timestep
+        print "velocity",self.velocity
+        """
+
+        import advection_ext            
+        self.timestep = advection_ext.compute_fluxes(self, Stage, huge_timestep, max_timestep)
+
+
+
+    def evolve(self,
+               yieldstep = None,
+               finaltime = None,
+               duration = None,
+               skip_initial_step = False):
+
+        """Specialisation of basic evolve method from parent class
+        """
+
+        #Call basic machinery from parent class
+        for t in Generic_domain.evolve(self,
+                                       yieldstep=yieldstep,
+                                       finaltime=finaltime,
+                                       duration=duration,
+                                       skip_initial_step=skip_initial_step):
+
+            #Pass control on to outer loop for more specific actions
+            yield(t)
+
+
 
 
@@ -207 +271 @@
         self.timestep = timestep
 
-
-    def compute_fluxes_ext(self):
-        """Compute all fluxes and the timestep suitable for all volumes
-        in domain.
-
-        Compute total flux for each conserved quantity using "flux_function"
-
-        Fluxes across each edge are scaled by edgelengths and summed up
-        Resulting flux is then scaled by area and stored in
-        domain.explicit_update
-
-        The maximal allowable speed computed by the flux_function
-        for each volume
-        is converted to a timestep that must not be exceeded. The minimum of
-        those is computed as the next overall timestep.
-
-        Post conditions:
-        domain.explicit_update is reset to computed flux values
-        domain.timestep is set to the largest step satisfying all volumes.
-        """
-
-        import sys
-        from Numeric import zeros, Float
-        from anuga.config import max_timestep
-
-        ntri = len(self)
-
-        timestep = max_timestep
-
-        #Shortcuts
-        Stage = self.quantities['stage']
-
-        #Arrays
-        neighbours      = self.neighbours
-        neighbour_edges = self.neighbour_edges
-        normals         = self.normals
-        areas           = self.areas
-        radii           = self.radii
-        edgelengths     = self.edgelengths
-        tri_full_flag   = self.tri_full_flag
-
-        stage_edge      = Stage.edge_values
-        stage_bdry      = Stage.boundary_values
-        stage_update    = Stage.explicit_update
-
-        huge_timestep = float(sys.maxint)
-
-        v = self.velocity
-
-        nbdry = len(stage_bdry)
-
-        from advection_ext import compute_fluxes
-
-        print "stage_update",stage_update
-        print "stage_edge",stage_edge
-        print "stage_bdry",stage_bdry
-        print "neighbours",neighbours
-        print "neighbour_edges",neighbour_edges
-        print "normals",normals
-        print "areas",areas
-        print "radii",radii
-        print "edgelengths",edgelengths
-        print "tri_full_flag",tri_full_flag
-        print "huge_timestep",huge_timestep
-        print "max_timestep",max_timestep
-        print "v",v
-        print "ntri",ntri
-        print "nbdry",nbdry
-
-
-                
-        timestep = compute_fluxes(stage_update,stage_edge,stage_bdry,
-                                  neighbours,neighbour_edges,normals,
-                                  areas,radii,edgelengths,
-                                  tri_full_flag,
-                                  huge_timestep,max_timestep,v,ntri,nbdry)
-
-        print "stage_update",stage_update        
-        
-        #print 'timestep out2 =',timestep
-
-        self.timestep = timestep
-
-
-    def evolve(self,
-               yieldstep = None,
-               finaltime = None,
-               duration = None,
-               skip_initial_step = False):
-
-        """Specialisation of basic evolve method from parent class
-        """
-
-        #Call basic machinery from parent class
-        for t in Generic_domain.evolve(self,
-                                       yieldstep=yieldstep,
-                                       finaltime=finaltime,
-                                       duration=duration,
-                                       skip_initial_step=skip_initial_step):
-
-            #Pass control on to outer loop for more specific actions
-            yield(t)

anuga_core/source/anuga/advection/advection_ext.c

@@ +1 @@
+// Python - C extension module for advection.py
+//
+// To compile (Python2.X):
+// use python ../utilities/compile.py
+//
+// See the module advection.py
+//
+//
+// Steve Roberts, ANU 2008
+
+
+#include "Python.h"
+#include "Numeric/arrayobject.h"
 #include "math.h"
 #include "stdio.h"
 
-void  print_double_array(char* name, double* array, int n, int m){
-
-    int k,i,km;
-
-    printf("%s = [",name);
-    for (k=0; k<n; k++){
-        km = m*k;
-        printf("[");
-        for (i=0; i<m ; i++){
-            printf("%g ",array[km+i]);
-        }
-        if (k==(n-1))
-            printf("]");
-        else
-            printf("]\n");
-    }
-    printf("]\n");
-}
-
-void  print_int_array(char* name, int* array, int n, int m){
-
-    int k,i,km;
-
-    printf("%s = [",name);
-    for (k=0; k<n; k++){
-        km = m*k;
-        printf("[");
-        for (i=0; i<m ; i++){
-            printf("%i ",array[km+i]);
-        }
-        if (k==(n-1))
-            printf("]");
-        else
-            printf("]\n");
-    }
-    printf("]\n");
-}
-
-
-double compute_fluxes(
-                    double* stage_update,
-                    double* stage_edge,
-                    double* stage_bdry,
-                    int* neighbours,
-                    int* neighbour_edges,
-                    double* normals,
-                    double* areas,
-                    double* radii,
-                    double* edgelengths,
-                    int*   tri_full_flag,
+// Shared code snippets
+#include "util_ext.h"
+
+
+
+double _compute_fluxes(
+                    double* quantity_update,
+                    double* quantity_edge,
+                    double* quantity_bdry,
+                    long*   domain_neighbours,
+                    long*   domain_neighbour_edges,
+                    double* domain_normals,
+                    double* domain_areas,
+                    double* domain_radii,
+                    double* domain_edgelengths,
+                    long*   domain_tri_full_flag,
+                    double* domain_velocity,
                     double  huge_timestep,
                     double  max_timestep,
-                    double* v,
                     int ntri,
                     int nbdry){
+
+ 
         //Loop
 
@@ -72 +53 @@
         timestep = max_timestep;
 
-
-        printf("======================================================\n");
-        printf("INSIDE compute_fluxes\n");
-
-
-        print_double_array( "stage_update",stage_update, ntri, 1);
-        print_double_array( "stage_edge",stage_edge, ntri, 3);
-        print_double_array( "stage_bdry",stage_bdry, nbdry, 1);
-        print_int_array( "neighbours",neighbours, ntri, 3);
-        print_int_array( "neighbour_edges",neighbour_edges, ntri, 3);
-        print_double_array( "normals",normals, ntri, 6);
-        print_double_array( "areas",areas, ntri, 1);
-        print_double_array( "radii",radii, ntri, 1);
-        print_double_array( "edgelengths",edgelengths, ntri, 3);
-        print_int_array( "tri_full_flag",tri_full_flag, ntri, 1);
-        printf("huge_timestep = %g\n",huge_timestep);
-        printf("max_timestep = %g\n",max_timestep);
-        print_double_array( "v",v, 2, 1);
-        printf("ntri = %i\n",ntri);
-        printf("nbdry = %i\n",nbdry);
-
-
         for (k=0; k<ntri; k++){
             optimal_timestep = huge_timestep;
@@ -101 +60 @@
                 k_i = k3+i;
                 //Quantities inside volume facing neighbour i
-                ql = stage_edge[k_i];
+                ql = quantity_edge[k_i];
 
 
                 //Quantities at neighbour on nearest face
-                n = neighbours[k_i];
+                n = domain_neighbours[k_i];
                 if (n < 0) {
                     m = -n-1; //Convert neg flag to index
-                    qr = stage_bdry[m];
+                    qr = quantity_bdry[m];
                 } else {
-                    m = neighbour_edges[k_i];
+                    m = domain_neighbour_edges[k_i];
                     n_m = 3*n+m;
-                    qr = stage_edge[n_m];
+                    qr = quantity_edge[n_m];
                 }
 
@@ -119 +78 @@
                 for (j=0; j<2; j++){
                     k_i_j = 6*k+2*i+j;
-                    normal[j] = normals[k_i_j];
+                    normal[j] = domain_normals[k_i_j];
                 }
 
 
                 //Flux computation using provided function
-                normal_velocity = v[0]*normal[0] + v[1]*normal[1];
+                normal_velocity = domain_velocity[0]*normal[0] + domain_velocity[1]*normal[1];
 
                 if (normal_velocity < 0) {
@@ -133 +92 @@
 
                 max_speed = fabs(normal_velocity);
-                flux = flux - edgeflux * edgelengths[k_i];
+                flux = flux - edgeflux * domain_edgelengths[k_i];
 
                 //Update optimal_timestep
-                if (tri_full_flag[k] == 1) {
+                if (domain_tri_full_flag[k] == 1) {
                     if (max_speed != 0.0) {
-                        optimal_timestep = (optimal_timestep>radii[k]/max_speed) ? radii[k]/max_speed : optimal_timestep;
+                        optimal_timestep = (optimal_timestep>domain_radii[k]/max_speed) ? 
+                          domain_radii[k]/max_speed : optimal_timestep;
                     }
                 }
@@ -146 +106 @@
             //Normalise by area and store for when all conserved
             //quantities get updated
-            stage_update[k] = flux/areas[k];
+            quantity_update[k] = flux/domain_areas[k];
 
             timestep = (timestep>optimal_timestep) ? optimal_timestep : timestep;
 
         }
-
-        //printf("timestep out = %g\n",timestep);
-
-
-
-        printf("INSIDE compute_fluxes, end \n");
-
-        print_double_array( "stage_update",stage_update, ntri, 1);
-
-        printf("FINISHED compute_fluxes\n");
-        printf("======================================================\n");
-
 
         return timestep;
 }
 
+
+//-----------------------------------------------------
+// Python method Wrappers 
+//-----------------------------------------------------
+
+
+
+PyObject *compute_fluxes(PyObject *self, PyObject *args) {
+  /*Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for the conserved quantity stage.
+
+    The maximal allowable speed computed for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Python call:
+    timestep = advection_ext.compute_fluxes(domain,quantity,huge_timestep,max_timestep)
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+
+  */
+
+  PyObject *domain, *quantity;
+  PyArrayObject 
+    * quantity_update,
+    * quantity_edge,
+    * quantity_bdry,
+    * domain_neighbours,
+    * domain_neighbour_edges,
+    * domain_normals,
+    * domain_areas,
+    * domain_radii,
+    * domain_edgelengths,
+    * domain_tri_full_flag,
+    * domain_velocity;
+                  
+
+  // Local variables
+  int ntri, nbdry;
+  double huge_timestep, max_timestep;
+  double timestep;
+
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOdd",
+                        &domain,
+                        &quantity,
+                        &huge_timestep,
+                        &max_timestep)) {
+    PyErr_SetString(PyExc_RuntimeError, "advection_ext.c: compute_fluxes could not parse input");
+    return NULL;
+  }
+
+  quantity_edge          = get_consecutive_array(quantity, "edge_values");
+  quantity_bdry          = get_consecutive_array(quantity, "boundary_values");
+  quantity_update        = get_consecutive_array(quantity, "explicit_update");
+  domain_neighbours      = get_consecutive_array(domain,   "neighbours");
+  domain_neighbour_edges = get_consecutive_array(domain,   "neighbour_edges");
+  domain_normals         = get_consecutive_array(domain,   "normals");
+  domain_areas           = get_consecutive_array(domain,   "areas");
+  domain_radii           = get_consecutive_array(domain,   "radii");
+  domain_edgelengths     = get_consecutive_array(domain,   "edgelengths");
+  domain_tri_full_flag   = get_consecutive_array(domain,   "tri_full_flag");
+  domain_velocity        = get_consecutive_array(domain,   "velocity");  
+    
+  ntri  = quantity_edge -> dimensions[0];
+  nbdry = quantity_bdry -> dimensions[0]; 
+
+  timestep = _compute_fluxes((double*) quantity_update -> data,
+                             (double*) quantity_edge -> data,
+                             (double*) quantity_bdry -> data,
+                             (long*)   domain_neighbours -> data,
+                             (long*)   domain_neighbour_edges -> data,
+                             (double*) domain_normals -> data,
+                             (double*) domain_areas ->data,
+                             (double*) domain_radii -> data,
+                             (double*) domain_edgelengths -> data,
+                             (long*)   domain_tri_full_flag -> data,
+                             (double*) domain_velocity -> data,
+                             huge_timestep,
+                             max_timestep,
+                             ntri,
+                             nbdry);
+
+  // Release and return
+  Py_DECREF(quantity_update);
+  Py_DECREF(quantity_edge);
+  Py_DECREF(quantity_bdry);
+  Py_DECREF(domain_neighbours);
+  Py_DECREF(domain_neighbour_edges);
+  Py_DECREF(domain_normals);
+  Py_DECREF(domain_areas);
+  Py_DECREF(domain_radii);
+  Py_DECREF(domain_edgelengths);
+  Py_DECREF(domain_tri_full_flag);
+  Py_DECREF(domain_velocity);
+
+
+  return Py_BuildValue("d", timestep);
+}
+
+
+
+//-------------------------------
+// Method table for python module
+//-------------------------------
+static struct PyMethodDef MethodTable[] = {
+  /* The cast of the function is necessary since PyCFunction values
+   * only take two PyObject* parameters, and rotate() takes
+   * three.
+   */
+
+  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initadvection_ext(void){
+  Py_InitModule("advection_ext", MethodTable);
+  import_array(); // Necessary for handling of NumPY structures
+}

anuga_core/source/anuga/advection/test_advection.py

@@ -60 +60 @@
         U = -domain.quantities['stage'].explicit_update
         R = -0.5/domain.areas[0]
-
-        print 'U ', U
-        print 'R ', R
-        
 
         assert U==R, '%s %s' %(U, R)

anuga_core/source/anuga/compile_all.py

@@ -7 +7 @@
 os.chdir('utilities')
 execfile('compile.py')
+
+os.chdir('..')
+os.chdir('advection')
+execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 
 os.chdir('..')

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -1521 +1521 @@
     *xmom_centroid_values,
     *ymom_centroid_values,
-        *elevation_centroid_values,
+    *elevation_centroid_values,
     *vertex_coordinates,
     *stage_vertex_values,
     *xmom_vertex_values,
     *ymom_vertex_values,
-        *elevation_vertex_values;
-  PyObject *domain, *Tmp;
+    *elevation_vertex_values;
+  
+  PyObject *domain;
 
   
@@ -1561 +1562 @@
   // This is used in the limiting process
   
-  
-  Tmp = PyObject_GetAttrString(domain, "beta_w");
-  if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain");
-    return NULL;
-  }  
-  beta_w = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp);
-  
-  Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
-  if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain");
-    return NULL;
-  }  
-  beta_w_dry = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp);
-  
-  Tmp = PyObject_GetAttrString(domain, "beta_uh");
-  if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain");
-    return NULL;
-  }  
-  beta_uh = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp);
-  
-  Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
-  if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain");
-    return NULL;
-  }  
-  beta_uh_dry = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp); 
-
-  Tmp = PyObject_GetAttrString(domain, "beta_vh");
-  if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain");
-    return NULL;
-  }  
-  beta_vh = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp);
-  
-  Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
-  if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain");
-    return NULL;
-  }  
-  beta_vh_dry = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp);
-  
-  Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
-  if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_height");
-    return NULL;
-  }  
-  minimum_allowed_height = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp);  
-
-  Tmp = PyObject_GetAttrString(domain, "epsilon");
-  if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object epsilon");
-    return NULL;
-  }  
-  epsilon = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp);  
-  
-  
+
+  beta_w      = get_double(domain,"beta_w");
+  beta_w_dry  = get_double(domain,"beta_w_dry");
+  beta_uh     = get_double(domain,"beta_uh");
+  beta_uh_dry = get_double(domain,"beta_uh_dry");
+  beta_vh     = get_double(domain,"beta_vh");
+  beta_vh_dry = get_double(domain,"beta_vh_dry");  
+
+  minimum_allowed_height = get_double(domain,"minimum_allowed_height");
+  epsilon = get_double(domain,"epsilon");
+
+  number_of_elements = stage_centroid_values -> dimensions[0];  
+
   // Call underlying computational routine
-  number_of_elements = stage_centroid_values -> dimensions[0];  
   e = _extrapolate_second_order_sw(number_of_elements,
                                    epsilon,

anuga_core/source/anuga/test_all.py

@@ -24 +24 @@
 
 # Exclude test_advection on all platforms for the time being. See ticket:205
-exclude_files.append('test_advection.py') #Weave doesn't work on Linux
+#exclude_files.append('test_advection.py') #Weave doesn't work on Linux
 
 

anuga_core/source/anuga/utilities/util_ext.h

@@ -166 +166 @@
 
 
-
+void  print_double_array(char* name, double* array, int n, int m){
+
+    int k,i,km;
+
+    printf("%s = [",name);
+    for (k=0; k<n; k++){
+        km = m*k;
+        printf("[");
+        for (i=0; i<m ; i++){
+            printf("%g ",array[km+i]);
+        }
+        if (k==(n-1))
+            printf("]");
+        else
+            printf("]\n");
+    }
+    printf("]\n");
+}
+
+void  print_int_array(char* name, int* array, int n, int m){
+
+    int k,i,km;
+
+    printf("%s = [",name);
+    for (k=0; k<n; k++){
+        km = m*k;
+        printf("[");
+        for (i=0; i<m ; i++){
+            printf("%i ",array[km+i]);
+        }
+        if (k==(n-1))
+            printf("]");
+        else
+            printf("]\n");
+    }
+    printf("]\n");
+}
+
+
+void  print_long_array(char* name, long* array, int n, int m){
+
+    int k,i,km;
+
+    printf("%s = [",name);
+    for (k=0; k<n; k++){
+        km = m*k;
+        printf("[");
+        for (i=0; i<m ; i++){
+          printf("%i ",(int) array[km+i]);
+        }
+        if (k==(n-1))
+            printf("]");
+        else
+            printf("]\n");
+    }
+    printf("]\n");
+}
 
 void print_numeric_array(PyArrayObject *x) {  
@@ -217 +273 @@
 }
 
+double get_double(PyObject *O, char *name) {
+  PyObject *TObject;
+  double tmp;
+  
+
+  //Get double from attribute
+  TObject = PyObject_GetAttrString(O, name);
+  if (!TObject) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_double could not obtain double from object");
+    return 0.0;
+  }  
+  
+  tmp = PyFloat_AsDouble(TObject);
+  
+  Py_DECREF(TObject);
+  
+  return tmp;
+}
+
+
+