Commit 7cd385f7 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

Add python interface

parent 31d88f1e
......@@ -4,7 +4,7 @@ set(PYTHON_CXX_EXTENSION_HDRS scalfmm_interface.hpp scalfmm_python.hpp scalfmm_p
SET_SOURCE_FILES_PROPERTIES(${PYTHON_CXX_EXTENSION_SRCS} PROPERTIES LANGUAGE CXX )
set(CMAKE_CXX_FLAGS)
add_library(python SHARED ${PYTHON_CXX_EXTENSION_SRCS} )
add_library(python SHARED ${PYTHON_CXX_EXTENSION_SRCS} ${PYTHON_CXX_EXTENSION_HDRS} )
target_compile_options(python PRIVATE "-Wold-style-cast")
add_dependencies(python scalfmm)
target_link_libraries(python PUBLIC ${CMAKE_PROJECT_NAME} )
......
......@@ -5,6 +5,7 @@
#include <Python.h>
#include "numpy/ndarrayobject.h"
#include "scalfmm_interface.hpp"
/*!
* Methods exported by the extension.
*/
......
......@@ -6,8 +6,8 @@
#include <Python.h>
#include <numpy/ndarrayobject.h>
#include "scalfmm_python.hpp"
#include "Utils/FPoint.hpp"
#
#include "Utils/FPoint.hpp"
#include "Adaptive/FBox.hpp"
/////////////////////////////////////////////////////////////////////////////////////////////////////
......@@ -18,42 +18,47 @@ static char docstring_sorted[] =
"\n"
"Parameters\n"
" points 2-darray (nbELts, dim) of points (array of float)\n"
" level (int) levelto build the morton index if value=0 then \n"
" Centre 1-darray (dim) of centre of the cubic box\n"
" length (double) the size of the cubic box\n"
" level (int) level to build the morton index if value=0 then \n"
" we consider the maximal level (20)\n"
" \n"
"\n"
"Returns\n"
" None";
static PyObject *ScalFMMsorted(PyObject *self, PyObject *args) {
PyArrayObject *py_points=nullptr, *py_centre=nullptr,*py_length=nullptr;
PyArrayObject *py_points = nullptr, *py_centre = nullptr;
double length ;
int level ;
if (! PyArg_ParseTuple( args, "OOOi",&py_points, &py_centre,&py_length, &level) ) {
if (! PyArg_ParseTuple( args, "OOdi",&py_points, &py_centre,&length, &level) ) {
return nullptr;
}
FPoint<double>* points_data = static_cast<FPoint<double>*>(PyArray_DATA(py_points));
FPoint<double> centre = *(static_cast<FPoint<double>*>(PyArray_DATA(py_centre)));
double* length = static_cast<double*>(PyArray_DATA(py_length));
const int nbElts = static_cast<int>(PyArray_DIM(py_points, 0));
const int dim = static_cast<int>(PyArray_DIM(py_points, 1));
std::cout << "dim " << dim << std::endl;
std::cout << "nbElts " << nbElts << std::endl;
std::cout << "dim " << dim << std::endl;
std::cout << "centre " << centre <<std::endl;
std::cout << "length " << length[0] << " " << length[1] << " " << length[2] <<std::endl;
std::cout << "length " << length <<std::endl;
std::cout << "level " << level <<std::endl;
//
// Build the Box
// Allocate the morton
npy_intp dims[1]={nbElts};
std::vector< FPoint<double>> vec_points(points_data,points_data+nbElts); //COPY HERE !!!!
PyArrayObject* py_morton = reinterpret_cast<PyArrayObject*>((PyArray_SimpleNew(1,dims,NPY_INT64)));
MortonIndex* morton_data = static_cast<MortonIndex*>(PyArray_DATA(py_morton));
PyArrayObject* py_perm = reinterpret_cast<PyArrayObject*>((PyArray_SimpleNew(1,dims,NPY_INT)));
std::vector<std::int64_t> vec_morton(nbElts,-1);
std::vector<std::int32_t> permutation(nbElts,-1);
//
const FBox<FPoint<double>> box{length[0],centre};
const FBox<FPoint<double>> box{length, centre};
//
//
std::vector< FPoint<double>> vec_points(points_data,points_data+nbElts); //COPY HERE !!!!
scalfmm::python::sortArrayWithMortonIndex(box, vec_points , vec_morton, permutation) ;
//
std::move(vec_morton.begin(), vec_morton.end(), morton_data);
......@@ -61,10 +66,7 @@ static PyObject *ScalFMMsorted(PyObject *self, PyObject *args) {
int* py_perm_data = static_cast<int*>(PyArray_DATA(py_perm));
std::move(permutation.begin(), permutation.end(), py_perm_data);
std::cout << "\n\n"<<std::endl;
for(int i = 0; i < 5 ; ++i){
std::cout << i << " "<< vec_points[i] << " " << morton_data[i] <<std::endl;
}
Py_INCREF(py_perm);
Py_INCREF(py_morton);
return PyTuple_Pack(2, py_morton, py_perm) ;
......
......@@ -5,10 +5,10 @@ namespace scalfmm {
template<class BOX_T, class POINTS_T, typename MORTON_T>
void sortPoints(const BOX_T &box,POINTS_T& points, MORTON_T& mortonIdx) {
const int nbParticles = mortonIdx.size() ;
const std::size_t max_level = 20 ; //sizeof(MortonIndex) * 8 / dim;
std::cout << "Nb elements: " << mortonIdx.size() << std::endl;
std::cout << "first points: " <<points[0] << std::endl;
// const auto nbParticles = mortonIdx.size() ;
// const std::size_t max_level = 20 ; //sizeof(MortonIndex) * 8 / dim;
// std::cout << "Nb elements: " << mortonIdx.size() << std::endl;
// std::cout << "first points: " <<points[0] << std::endl;
// std::vector<PARTICLE_T> vect(nbParticles);
// for(int i = 0; i < nbParticles ; ++i){
// vect[i].val = array[i] ;
......
# # Example # import sys import numpy as np sys.path.append("../Build/pyscalfmm/c_extension") import pyscalfmm # # Generate data distribution # Build tree with data --> linearTree (thepoints are sorted according # to their Morton Index ) # Select the method (Kernel, Interpolation points, order) # Build the real tree and compute the potential and the forces # Save the result # # Other methods # - compute the direct interactions by # results = directInteraction(Kernel, points) # results contains the potential or the forces or both. # CoinMin=np.amin(points, axis=0) # CoinMax=np.amax(points, axis=0) # CLength =(CoinMax-CoinMin) # CCenter =0.5*(CoinMax+CoinMin) # # dim = 3 # dimension of the space N = 100 # number of points points = np.random.uniform(-1.0,1.0,N*dim) #points.reshape([dim,N]) # x_1, ...x_N, y_1, ... y_N, z_1, ..., z_N points=np.reshape(points,[N,dim]) # x_1,y_1,z_1, ...x_N, y_N, z_N print(points.shape) print("points 0 ",points[0,:]) Centre=np.array([0.0,0.0,0.0]) L=np.array([2.0,2.0,2.0]) morton,perm = pyscalfmm.sorted(points,Centre,L,4) print("\n points 0 ",points[0,:]) print("\n Morton ",morton[0:5]) print("\n Permutation ",perm[0:5]) print(type(morton),type(morton[0]),morton.shape) print(type(perm),type(perm[0]),perm.shape) charges = np.random.uniform(-1,1,N) # compute the physical values method="treeHeigh=5, points=equispaced, order=3" forces = pyscalfmm.computeforce(points,charges,method) print("Force: ",forces.shape) print("Points: ",points.shape) pyscalfmm.write("result.fma",points,charges,forces) print("--- END ---")
\ No newline at end of file
# # Example # import sys import numpy as np sys.path.append("../Build/pyscalfmm/c_extension") import pyscalfmm # # Generate data distribution # Build tree with data --> linearTree (thepoints are sorted according # to their Morton Index ) # Select the method (Kernel, Interpolation points, order) # Build the real tree and compute the potential and the forces # Save the result # # Other methods # - compute the direct interactions by # results = directInteraction(Kernel, points) # results contains the potential or the forces or both. # CoinMin=np.amin(points, axis=0) # CoinMax=np.amax(points, axis=0) # CLength =(CoinMax-CoinMin) # CCenter =0.5*(CoinMax+CoinMin) # # dim = 3 # dimension of the space N = 100 # number of points points = np.random.uniform(-1.0,1.0,N*dim) #points.reshape([dim,N]) # x_1, ...x_N, y_1, ... y_N, z_1, ..., z_N points=np.reshape(points,[N,dim]) # x_1,y_1,z_1, ...x_N, y_N, z_N print(points.shape) print("points 0 ",points[0,:]) Centre=np.array([0.0,0.0,0.0]) L=np.array([2.0,2.0,2.0]) morton,perm = pyscalfmm.sorted(points,Centre,L[0],8) print("\n points 0 ",points[0,:]) print("\n Morton ",morton[0:5]) print("\n Permutation ",perm[0:5]) print(type(morton),type(morton[0]),morton.shape) print(type(perm),type(perm[0]),perm.shape) charges = np.random.uniform(-1,1,N) # compute the physical values method="treeHeigh=5, points=equispaced, order=3" forces = pyscalfmm.computeforce(points,charges,method) print("Force: ",forces.shape) print("Points: ",points.shape) pyscalfmm.write("result.fma",points,charges,forces) print("--- END ---")
\ No newline at end of file
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment