Working swig-numpy conversion

src/arma_numpy.cpp

@@ -0,0 +1,38 @@
+#define PY_ARRAY_UNIQUE_SYMBOL npy_array_api
+#define NO_IMPORT_ARRAY
+#include "arma_numpy.h"
+
+// Annes conversion functions. In a later stage they have to be
+// generalized for arrays of arbitrary dimensions
+
+
+
+PyObject* npy_from_vd(const vd& in){
+  long int size=in.size();
+  npy_intp dims[1]={size};
+  // This code should be improved massively?
+  if(size==0){
+    std::cout << "Size is zero!\n";
+    return NULL;
+  }
+  cout << "SFSF\n";    
+  PyArrayObject* array =  (PyArrayObject*) \
+          PyArray_SimpleNew(1, dims, NPY_DOUBLE); 
+  cout << "SFSF\n";    
+  if (array == NULL) {
+    return NULL;
+  }
+  double *pydat= (double*) PyArray_DATA(array);
+  cout << "SFSF\n";
+  mempcpy(pydat,in.memptr(),size*sizeof(double));
+  return PyArray_Return(array);
+}
+vd vd_from_npy(const PyArrayObject* const in){
+
+  double* pydata = (double*) PyArray_DATA(in);
+  vd result;
+  // long int dims=
+  // PyArra
+  //   vd res()
+  return result;
+}

src/arma_numpy.h

@@ -0,0 +1,13 @@
+#pragma once
+#ifndef ARMA_NUMPY_H
+#define ARMA_NUMPY_H
+#include <Python.h>
+#include <numpy/arrayobject.h>
+#include "vtypes.h"
+SPOILNAMESPACE
+vd vd_from_npy(PyArrayObject const * const);
+PyObject* npy_from_vd(const vd& armavec);
+
+
+
+#endif // ARMA_NUMPY_H

src/solutioninstance.cpp

@@ -14,6 +14,12 @@ namespace td{
       gp->setData(rho);
     }
   } // setrho()
+  vd SolutionInstance::getrho() const{
+    return gpsol_to_python(gps,&SolutionAtGp::rho);
+  }
+  vd SolutionInstance::getp() const{
+    return gpsol_to_python(gps,&SolutionAtGp::p);
+  }
   int SolutionInstance::save(std::ofstream& str,data d){
     TRACE(15,"SolutionInstance::save()");
     if(d==data::p)

src/solutioninstance.h

@@ -2,9 +2,21 @@
 #include "solutionatgp.h"
 #include "vtypes.h"
 
+void do_arma_to_python(const vd& armavec,double** data,int* n);
+
 namespace td{
   SPOILNAMESPACE
 
+  template<typename funptr>
+  vd gpsol_to_python(const vector<SolutionAtGp>& gps,funptr fun){
+    int size=gps.size();
+    vd result(size);
+    for(int i=0;i<size;i++)
+      result(i)=((gps[i]).*fun)();
+    return result;
+  }
+
+
   enum data{
     p,u
   };
@@ -20,6 +32,8 @@ namespace td{
     void setData(us i,SolutionAtGp& sol){ gps.at(i)=sol;}
     SolutionAtGp& get(us i){ return gps.at(i); }
 
+    vd getp() const;
+    vd getrho() const;
     void setrho(d rho);
     int save(std::ofstream&,data d);
   };

src/tube.cpp

@@ -5,7 +5,8 @@ namespace td{
 
   d Tube::pleft(d t){
     d pleft=std::sin(gc.getomg()*t);
-    // cout << "t:"<< t<<endl;
+    cout << "t:"<< t<<endl;
+    cout << "omg:" << gc.getomg()<< endl;
     // cout << "pleft:"<< pleft<<endl;
     return pleft;
 
@@ -26,6 +27,7 @@ namespace td{
     int integrations=0;
     while(integrations<n){
       Integrate(dt);
+      integrations++;
     }
   }
   void Tube::Integrate(d dt){

timedomain.i

@@ -1,7 +1,14 @@
 %module timedomaineuler
 %{
+  #define PY_ARRAY_UNIQUE_SYMBOL npy_array_api
+  /* #define NO_IMPORT_ARRAY */
+  #define SWIG_FILE_WITH_INIT
+  #include <cerrno>
+  #include <cstdlib>
   #include "vtypes.h"
+  #include "arma_numpy.h"
   #include "tube.h"
+
   #include "globalconf.h"
   SPOILNAMESPACE
     namespace td{
@@ -9,7 +16,34 @@
   }
 %}
 
+/* Call import_array() on loading module*/
+%init%{
+  cout << "Import array called\n";
+  import_array();
+%}
 
+// A better way is the following
+// Convert from numpy to vd and vice versa
+class vd{
+ public:
+  virtual ~vd();
+};
+
+%typemap(in) const vd& INPUT {
+  PyArrayObject* arr = (PyArrayObject*) PyArray_FROM_OTF($input, NPY_DOUBLE, NPY_IN_ARRAY);
+  if(arr==NULL)
+    return NULL;
+  int ndims=PyArray_DIM(arr,0);
+  if(ndims!=1){
+    PyErr_SetString(PyExc_TypeError,"Number of dimensions not equal to 1");
+    return NULL;
+  }
+  $1= vd_from_npy(arr);
+}
+%typemap(out) vd {
+  cout << "Hoi\n";
+  $result=npy_from_vd($1);
+}
 typedef double d;
 typedef unsigned us;
 
@@ -83,6 +117,8 @@ namespace td{
     SolutionInstance(int gp);
     ~SolutionInstance(){}
     d getTime() const;
+    vd getrho() const;
+    vd getp() const;
     void setTime(d t);
     void setData(us i,SolutionAtGp& sol);
     SolutionAtGp& get(us i);

timedomain.ipynb

@@ -1,7 +1,7 @@
 {
  "metadata": {
   "name": "",
-  "signature": "sha256:15ea21addc5fb3913b7c7be02aef3ba5cdc13caaf5cb403564b9e56bc36c4974"
+  "signature": "sha256:c219b2e07da88eca6f031393997d37a4b29390fb7116605d50ca9f4c6d0512b5"
  },
  "nbformat": 3,
  "nbformat_minor": 0,
@@ -17,8 +17,7 @@
      ],
      "language": "python",
      "metadata": {},
-     "outputs": [],
+     "outputs": []
-     "prompt_number": 1
     },
     {
      "cell_type": "code",
@@ -32,44 +31,13 @@
       "gp=300\n",
       "gc=cvar.gc #Reference!\n",
       "dx=L/(gp-1);                    # One left and right gp, so\n",
-      "CFL=0.95;\n",
+      "CFL=0.9;\n",
       "gc.setfreq(f)\n",
       "tube=Tube(L,gp)\n",
       "dt=min(CFL*dx/gc.c0(),T/50)"
      ],
      "language": "python",
      "metadata": {},
-     "outputs": [],
-     "prompt_number": 2
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "gc.getfreq()"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [
-      {
-       "metadata": {},
-       "output_type": "pyout",
-       "prompt_number": 3,
-       "text": [
-        "85.785"
-       ]
-      }
-     ],
-     "prompt_number": 3
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "tube.DoIntegration()"
-     ],
-     "language": "python",
-     "metadata": {},
      "outputs": []
     },
     {
@@ -80,77 +48,44 @@
      ],
      "language": "python",
      "metadata": {},
-     "outputs": [],
+     "outputs": []
-     "prompt_number": 3
     },
     {
      "cell_type": "code",
      "collapsed": false,
      "input": [
-      "gp0=sol.get(0)"
+      "#sol.getrho()"
      ],
      "language": "python",
      "metadata": {},
-     "outputs": [],
+     "outputs": []
-     "prompt_number": 5
     },
     {
      "cell_type": "code",
      "collapsed": false,
      "input": [
-      "gp0.rho()"
+      "dt*gc.getfreq()"
      ],
      "language": "python",
      "metadata": {},
-     "outputs": [
+     "outputs": []
-      {
-       "metadata": {},
-       "output_type": "pyout",
-       "prompt_number": 8,
-       "text": [
-        "1.2043280930847857"
-       ]
-      }
-     ],
-     "prompt_number": 8
     },
     {
      "cell_type": "code",
      "collapsed": false,
      "input": [
-      "gp0.setData(1.222)"
+      "tube.DoIntegration(dt,10)\n",
+      "sol=tube.getSol()\n",
+      "rho=tube.getSol().getrho()\n",
+      "p=tube.getSol().getp()\n",
+      "figure(figsize=(12,4))\n",
+      "subplot(121)\n",
+      "plot(p)\n",
+      "subplot(122)\n",
+      "plot(rho)"
      ],
      "language": "python",
      "metadata": {},
-     "outputs": [],
-     "prompt_number": 11
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "gp0.rho()"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [
-      {
-       "metadata": {},
-       "output_type": "pyout",
-       "prompt_number": 12,
-       "text": [
-        "1.222"
-       ]
-      }
-     ],
-     "prompt_number": 12
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [],
-     "language": "python",
-     "metadata": {},
      "outputs": []
     }
    ],