diff --git a/src/arma_numpy.cpp b/src/arma_numpy.cpp new file mode 100644 index 0000000..2761a69 --- /dev/null +++ b/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; +} diff --git a/src/arma_numpy.h b/src/arma_numpy.h new file mode 100644 index 0000000..53b586c --- /dev/null +++ b/src/arma_numpy.h @@ -0,0 +1,13 @@ +#pragma once +#ifndef ARMA_NUMPY_H +#define ARMA_NUMPY_H +#include +#include +#include "vtypes.h" +SPOILNAMESPACE +vd vd_from_npy(PyArrayObject const * const); +PyObject* npy_from_vd(const vd& armavec); + + + +#endif // ARMA_NUMPY_H diff --git a/src/solutioninstance.cpp b/src/solutioninstance.cpp index da2e80f..87c7490 100644 --- a/src/solutioninstance.cpp +++ b/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) diff --git a/src/solutioninstance.h b/src/solutioninstance.h index 40c540d..60cd1e1 100644 --- a/src/solutioninstance.h +++ b/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 + vd gpsol_to_python(const vector& gps,funptr fun){ + int size=gps.size(); + vd result(size); + for(int i=0;i + #include #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); diff --git a/timedomain.ipynb b/timedomain.ipynb index 0d476f9..b2649f8 100644 --- a/timedomain.ipynb +++ b/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": [], - "prompt_number": 1 + "outputs": [] }, { "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": [], - "prompt_number": 3 + "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ - "gp0=sol.get(0)" + "#sol.getrho()" ], "language": "python", "metadata": {}, - "outputs": [], - "prompt_number": 5 + "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ - "gp0.rho()" + "dt*gc.getfreq()" ], "language": "python", "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 8, - "text": [ - "1.2043280930847857" - ] - } - ], - "prompt_number": 8 + "outputs": [] }, { "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": [] } ],