diff --git a/src/Python/CMakeLists.txt b/src/Python/CMakeLists.txt
index 6356af0046b95d4f5a61944324cdcf2827408270..598d9faf14446f4073daffccd4a9cc5f15beb7b7 100644
--- a/src/Python/CMakeLists.txt
+++ b/src/Python/CMakeLists.txt
@@ -1,8 +1,36 @@
 find_package( PythonInterp 3 )
+find_package( PythonLibs 3 )
+
+set( PYTHON_SITE_PACKAGES_DIR lib/python${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}/site-packages )
 
 if( PYTHONINTERP_FOUND )
    CONFIGURE_FILE( "__init__.py.in" "${PROJECT_BUILD_PATH}/Python/__init__.py" )
    INSTALL( FILES ${PROJECT_BUILD_PATH}/Python/__init__.py
                   LogParser.py
-            DESTINATION lib/python${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}/site-packages/TNL )
+            DESTINATION ${PYTHON_SITE_PACKAGES_DIR}/TNL )
+endif()
+
+if( PYTHONLIBS_FOUND )
+   # download and build pybind11 at configure time
+   configure_file(pybind11.cmake.in ${CMAKE_BINARY_DIR}/pybind11-download/CMakeLists.txt)
+   execute_process(COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" .
+      RESULT_VARIABLE result
+      WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/pybind11-download )
+   if(result)
+      message(FATAL_ERROR "CMake step for pybind11 failed: ${result}")
+   endif()
+   execute_process(COMMAND ${CMAKE_COMMAND} --build .
+      RESULT_VARIABLE result
+      WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/pybind11-download )
+   if(result)
+      message(FATAL_ERROR "Build step for pybind11 failed: ${result}")
+   endif()
+
+   # add the pybind11 subdirectory to provide the pybind11_add_module macro
+   add_subdirectory(${CMAKE_BINARY_DIR}/pybind11-src ${CMAKE_BINARY_DIR}/pybind11-build)
+
+   # add the subdirectory with our bindings
+   add_subdirectory(pytnl)
+else()
+   message( "The Python.h header file was not found, Python bindings will not be builg." )
 endif()
diff --git a/src/Python/pybind11.cmake.in b/src/Python/pybind11.cmake.in
new file mode 100644
index 0000000000000000000000000000000000000000..18f1aeac607f719a067a8a61c1c3449cfd64687b
--- /dev/null
+++ b/src/Python/pybind11.cmake.in
@@ -0,0 +1,12 @@
+cmake_minimum_required(VERSION 2.8.2)
+
+project(pybind11-download)
+include(ExternalProject)
+
+ExternalProject_Add(pybind11
+  GIT_REPOSITORY    https://github.com/pybind/pybind11.git
+  GIT_TAG           master
+  SOURCE_DIR        "${CMAKE_BINARY_DIR}/pybind11-src"
+  BINARY_DIR        "${CMAKE_BINARY_DIR}/pybind11-build"
+  CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${CMAKE_INSTALL_PREFIX} -DPYBIND11_TEST=FALSE
+)
diff --git a/src/Python/pytnl/CMakeLists.txt b/src/Python/pytnl/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..555f8fc5e233653eab12064bfd88c8e791a6eb08
--- /dev/null
+++ b/src/Python/pytnl/CMakeLists.txt
@@ -0,0 +1,13 @@
+add_subdirectory( tnl )
+
+set( headers
+         exceptions.h
+         RawIterator.h
+         tnl_conversions.h
+         tnl_indexing.h
+         tnl_str_conversion.h
+         tnl_tuple_conversion.h
+         typedefs.h
+)
+
+install( FILES ${headers} DESTINATION "pytnl" )
diff --git a/src/Python/pytnl/RawIterator.h b/src/Python/pytnl/RawIterator.h
new file mode 100644
index 0000000000000000000000000000000000000000..6e9aa7b41cfb3f06c4630675062a5b9cbe77c848
--- /dev/null
+++ b/src/Python/pytnl/RawIterator.h
@@ -0,0 +1,51 @@
+#pragma once
+
+#include <iterator>
+
+template< typename DataType >
+class RawIterator : public std::iterator<std::random_access_iterator_tag,
+                                           DataType,
+                                           ptrdiff_t,
+                                           DataType*,
+                                           DataType&>
+{
+protected:
+    DataType*               m_ptr;
+
+public:
+    RawIterator( DataType* ptr = nullptr ) { m_ptr = ptr; }
+    RawIterator( const RawIterator<DataType> & rawIterator ) = default;
+    ~RawIterator(){}
+
+    RawIterator<DataType>&  operator=( const RawIterator<DataType> & rawIterator ) = default;
+    RawIterator<DataType>&  operator=( DataType* ptr ) { m_ptr = ptr; return (*this); }
+
+    operator                bool() const
+    {
+        if(m_ptr)
+            return true;
+        else
+            return false;
+    }
+
+    bool                    operator==( const RawIterator<DataType> & rawIterator ) const { return ( m_ptr == rawIterator.getConstPtr() ); }
+    bool                    operator!=( const RawIterator<DataType> & rawIterator ) const { return ( m_ptr != rawIterator.getConstPtr() ); }
+
+    RawIterator<DataType>&  operator+=( const ptrdiff_t & movement ){ m_ptr += movement; return (*this); }
+    RawIterator<DataType>&  operator-=( const ptrdiff_t & movement ){ m_ptr -= movement; return (*this); }
+    RawIterator<DataType>&  operator++() { ++m_ptr; return (*this); }
+    RawIterator<DataType>&  operator--() { --m_ptr; return (*this); }
+    RawIterator<DataType>   operator++( int ) { auto temp(*this); ++m_ptr; return temp; }
+    RawIterator<DataType>   operator--( int ) { auto temp(*this); --m_ptr; return temp; }
+    RawIterator<DataType>   operator+( const ptrdiff_t & movement ) { auto oldPtr = m_ptr; m_ptr+=movement; auto temp(*this); m_ptr = oldPtr; return temp; }
+    RawIterator<DataType>   operator-( const ptrdiff_t & movement ) { auto oldPtr = m_ptr; m_ptr-=movement; auto temp(*this); m_ptr = oldPtr; return temp; }
+
+    ptrdiff_t               operator-( const RawIterator<DataType>& rawIterator ) { return std::distance(rawIterator.getPtr(), this->getPtr()); }
+
+    DataType&               operator*() { return *m_ptr; }
+    const DataType&         operator*() const { return *m_ptr; }
+    DataType*               operator->() { return m_ptr; }
+
+    DataType*               getPtr() const { return m_ptr; }
+    const DataType*         getConstPtr() const { return m_ptr; }
+};
diff --git a/src/Python/pytnl/exceptions.h b/src/Python/pytnl/exceptions.h
new file mode 100644
index 0000000000000000000000000000000000000000..669fb3214c9e6bfad69bd2a80cde79ede532f781
--- /dev/null
+++ b/src/Python/pytnl/exceptions.h
@@ -0,0 +1,35 @@
+#pragma once
+
+#include <stdexcept>
+
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+
+#include <TNL/Assert.h>
+
+struct NotImplementedError
+   : public std::runtime_error
+{
+   NotImplementedError( const char* msg )
+   : std::runtime_error( msg )
+   {}
+};
+
+static void register_exceptions( py::module & m )
+{
+    py::register_exception_translator(
+        [](std::exception_ptr p) {
+            try {
+                if (p) std::rethrow_exception(p);
+            }
+            // translate exceptions used in the bindings
+            catch (const NotImplementedError & e) {
+                PyErr_SetString(PyExc_NotImplementedError, e.what());
+            }
+            // translate TNL::Assert::AssertionError
+            catch (const TNL::Assert::AssertionError & e) {
+                PyErr_SetString(PyExc_AssertionError, e.what());
+            }
+        }
+    );
+}
diff --git a/src/Python/pytnl/tnl/Array.h b/src/Python/pytnl/tnl/Array.h
new file mode 100644
index 0000000000000000000000000000000000000000..0d74acdccef35773616fb6b7e695caa32c767361
--- /dev/null
+++ b/src/Python/pytnl/tnl/Array.h
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+#include <pybind11/operators.h>
+namespace py = pybind11;
+
+#include "../tnl_indexing.h"
+
+#include <TNL/Containers/Array.h>
+
+template< typename ArrayType >
+void export_Array(py::module & m, const char* name)
+{
+    auto array = py::class_<ArrayType, TNL::Object>(m, name, py::buffer_protocol())
+        .def(py::init<>())
+        .def(py::init<int>())
+        .def_static("getType",              &ArrayType::getType)
+        .def("getTypeVirtual",              &ArrayType::getTypeVirtual)
+        .def_static("getSerializationType", &ArrayType::getSerializationType)
+        .def("getSerializationTypeVirtual", &ArrayType::getSerializationTypeVirtual)
+        .def("setSize", &ArrayType::setSize)
+        .def("setLike", &ArrayType::template setLike<ArrayType>)
+        .def("swap", &ArrayType::swap)
+        .def("reset", &ArrayType::reset)
+        .def("getSize", &ArrayType::getSize)
+        .def("setElement", &ArrayType::setElement)
+        .def("getElement", &ArrayType::getElement)
+        // operator=
+        .def("assign", []( ArrayType& array, const ArrayType& other ) -> ArrayType& {
+                return array = other;
+            })
+        .def(py::self == py::self)
+        .def(py::self != py::self)
+        .def("setValue", &ArrayType::setValue)
+
+        .def("__str__", []( ArrayType & a ) {
+                std::stringstream ss;
+                ss << a;
+                return ss.str();
+            } )
+
+        // Python buffer protocol: http://pybind11.readthedocs.io/en/master/advanced/pycpp/numpy.html
+        .def_buffer( [](ArrayType & a) -> py::buffer_info {
+            return py::buffer_info(
+                // Pointer to buffer
+                a.getData(),
+                // Size of one scalar
+                sizeof( typename ArrayType::ElementType ),
+                // Python struct-style format descriptor
+                py::format_descriptor< typename ArrayType::ElementType >::format(),
+                // Number of dimensions
+                1,
+                // Buffer dimensions
+                { a.getSize() },
+                // Strides (in bytes) for each index
+                { sizeof( typename ArrayType::ElementType ) }
+            );
+        })
+    ;
+
+    tnl_indexing< ArrayType >( array );
+    tnl_slice_indexing< ArrayType >( array );
+}
diff --git a/src/Python/pytnl/tnl/CMakeLists.txt b/src/Python/pytnl/tnl/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..6a284a02bb8b631ba1c649fe5b3a5cfa7aeca34f
--- /dev/null
+++ b/src/Python/pytnl/tnl/CMakeLists.txt
@@ -0,0 +1,26 @@
+set( sources
+      Grid1D.cpp
+      Grid2D.cpp
+      Grid3D.cpp
+#      Mesh.cpp
+      Object.cpp
+      SparseMatrix.cpp
+      String.cpp
+      tnl.cpp
+)
+pybind11_add_module( pytnl ${sources} )
+
+# link against tnl.so
+target_link_libraries( pytnl PRIVATE tnl )
+
+# rename the shared library to tnl.cpython-XXm-x86_64-linux-gnu.so
+set_target_properties( pytnl PROPERTIES LIBRARY_OUTPUT_NAME tnl )
+
+# We have bindings for unsafe objects (e.g. Array::operator[]) where assertion
+# is the only safeguard, so we need to translate the TNL::AssertionError to
+# Python's AssertionError.
+# NDEBUG is defined in the global CMAKE_CXX_FLAGS and cannot be easily removed
+# per-target, so we need to undefine it by passing -U NDEBUG.
+target_compile_options( pytnl PRIVATE -U NDEBUG -D TNL_THROW_ASSERTION_ERROR )
+
+install( TARGETS pytnl DESTINATION ${PYTHON_SITE_PACKAGES_DIR} )
diff --git a/src/Python/pytnl/tnl/EntityTypes.h b/src/Python/pytnl/tnl/EntityTypes.h
new file mode 100644
index 0000000000000000000000000000000000000000..1f10e2827dd2cc24d006c88f509e6b8d5a5cbf90
--- /dev/null
+++ b/src/Python/pytnl/tnl/EntityTypes.h
@@ -0,0 +1,36 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+
+enum class EntityTypes { Cell, Face, Vertex };
+
+inline void
+export_EntityTypes( py::module & m )
+{
+    // avoid duplicate conversion -> export only once
+    static bool exported = false;
+    if( ! exported ) {
+        // TODO: should be nested types instead
+        py::enum_< EntityTypes >( m, "EntityTypes" )
+            .value("Cell", EntityTypes::Cell)
+            .value("Face", EntityTypes::Face)
+            .value("Vertex", EntityTypes::Vertex)
+        ;
+        exported = true;
+    }
+}
+
+template< typename Mesh >
+typename Mesh::GlobalIndexType
+mesh_getEntitiesCount( const Mesh & self, const EntityTypes & entity )
+{
+    if( entity == EntityTypes::Cell )
+        return self.template getEntitiesCount< typename Mesh::Cell >();
+    else if( entity == EntityTypes::Face )
+        return self.template getEntitiesCount< typename Mesh::Face >();
+    else if( entity == EntityTypes::Vertex )
+        return self.template getEntitiesCount< typename Mesh::Vertex >();
+    else
+        throw py::value_error("The entity parameter must be either Cell, Face or Vertex.");
+}
diff --git a/src/Python/pytnl/tnl/Grid.h b/src/Python/pytnl/tnl/Grid.h
new file mode 100644
index 0000000000000000000000000000000000000000..afc5b39749362a08248befb7716c8f446e888dad
--- /dev/null
+++ b/src/Python/pytnl/tnl/Grid.h
@@ -0,0 +1,98 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+
+#include "StaticVector.h"
+#include "Grid_getSpaceStepsProducts.h"
+#include "EntityTypes.h"
+
+#include <type_traits>
+
+template< typename GridEntity, typename PyGrid >
+void export_GridEntity( PyGrid & scope, const char* name )
+{
+    typename GridEntity::CoordinatesType const & (GridEntity::* _getCoordinates1)(void) const = &GridEntity::getCoordinates;
+    typename GridEntity::CoordinatesType       & (GridEntity::* _getCoordinates2)(void) = &GridEntity::getCoordinates;
+
+    py::class_< GridEntity >( scope, name )
+        .def( py::init< typename GridEntity::GridType >(),
+              py::return_value_policy::reference_internal )
+        .def( py::init< typename GridEntity::GridType,
+                        typename GridEntity::CoordinatesType,
+                        typename GridEntity::EntityOrientationType,
+                        typename GridEntity::EntityBasisType >(),
+              py::return_value_policy::reference_internal )
+        .def("getEntityDimension", &GridEntity::getEntityDimension)
+        .def("getMeshDimension", &GridEntity::getMeshDimension)
+        // TODO: constructors
+        .def("getCoordinates", _getCoordinates1, py::return_value_policy::reference_internal)
+        .def("getCoordinates", _getCoordinates2, py::return_value_policy::reference_internal)
+        .def("setCoordinates", &GridEntity::setCoordinates)
+        .def("refresh", &GridEntity::refresh)
+        .def("getIndex", &GridEntity::getIndex)
+        // FIXME: some templates return reference, others value
+//        .def("getOrientation", &GridEntity::getOrientation, py::return_internal_reference<>(), py::return_value_policy<py::copy_const_reference>())
+        .def("setOrientation", &GridEntity::setOrientation)
+        // FIXME: some templates return reference, others value
+//        .def("getBasis", &GridEntity::getBasis, py::return_internal_reference<>(), py::return_value_policy<py::copy_const_reference>())
+        .def("setBasis", &GridEntity::setBasis)
+        // TODO: getNeighbourEntities
+        .def("isBoundaryEntity", &GridEntity::isBoundaryEntity)
+        .def("getCenter", &GridEntity::getCenter)
+        // FIXME: some templates return reference, others value
+//        .def("getMeasure", &GridEntity::getMeasure, py::return_value_policy<py::copy_const_reference>())
+        .def("getMesh", &GridEntity::getMesh, py::return_value_policy::reference_internal)
+    ;
+}
+
+template< typename Grid >
+void export_Grid( py::module & m, const char* name )
+{
+    // function pointers for overloaded methods
+// FIXME: number of parameters depends on the grid dimension
+//    void (Grid::* _setDimensions1)(const IndexType) = &Grid::setDimensions;
+    void (Grid::* _setDimensions2)(const typename Grid::CoordinatesType &) = &Grid::setDimensions;
+
+    export_EntityTypes(m);
+
+    auto grid = py::class_<Grid, TNL::Object>( m, name )
+        .def(py::init<>())
+        .def_static("getMeshDimension", &Grid::getMeshDimension)
+        .def_static("getType",              &Grid::getType)
+        .def("getTypeVirtual",              &Grid::getTypeVirtual)
+        .def_static("getSerializationType", &Grid::getSerializationType)
+        .def("getSerializationTypeVirtual", &Grid::getSerializationTypeVirtual)
+        // FIXME: number of parameters depends on the grid dimension
+//        .def("setDimensions", _setDimensions1)
+        .def("setDimensions", _setDimensions2)
+        .def("getDimensions", &Grid::getDimensions, py::return_value_policy::reference_internal)
+        .def("setDomain", &Grid::setDomain)
+        .def("getOrigin", &Grid::getOrigin, py::return_value_policy::reference_internal)
+        .def("getProportions", &Grid::getProportions, py::return_value_policy::reference_internal)
+        .def("getEntitiesCount", &mesh_getEntitiesCount< Grid >)
+        // TODO: if combined, the return type would depend on the runtime parameter (entity)
+        .def("getEntity_cell", &Grid::template getEntity<typename Grid::Cell>)
+        .def("getEntity_face", &Grid::template getEntity<typename Grid::Face>)
+        .def("getEntity_vertex", &Grid::template getEntity<typename Grid::Vertex>)
+        .def("getEntityIndex", &Grid::template getEntityIndex<typename Grid::Cell>)
+        .def("getEntityIndex", &Grid::template getEntityIndex<typename Grid::Face>)
+        .def("getEntityIndex", &Grid::template getEntityIndex<typename Grid::Vertex>)
+        .def("getCellMeasure", &Grid::getCellMeasure, py::return_value_policy::reference_internal)
+        .def("getSpaceSteps", &Grid::getSpaceSteps, py::return_value_policy::reference_internal)
+        .def("getSmallestSpaceStep", &Grid::getSmallestSpaceStep)
+        // TODO: writeProlog()
+    ;
+
+    // complicated methods
+    SpaceStepsProductsGetter< Grid >::export_getSpaceSteps( grid );
+
+    // nested types
+    export_StaticVector< typename Grid::CoordinatesType >( grid, "CoordinatesType" );
+    export_StaticVector< typename Grid::PointType >( grid, "PointType" );
+    export_GridEntity< typename Grid::Cell >( grid, "Cell" );
+    export_GridEntity< typename Grid::Face >( grid, "Face" );
+    // avoid duplicate conversion if the type is the same
+    if( ! std::is_same< typename Grid::Face, typename Grid::Vertex >::value )
+        export_GridEntity< typename Grid::Vertex >( grid, "Vertex" );
+}
diff --git a/src/Python/pytnl/tnl/Grid1D.cpp b/src/Python/pytnl/tnl/Grid1D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d10952b713e89a8deab65c6fe321bfeea5a3b16f
--- /dev/null
+++ b/src/Python/pytnl/tnl/Grid1D.cpp
@@ -0,0 +1,9 @@
+// conversions have to be registered for each object file
+#include "../tnl_conversions.h"
+
+#include "Grid.h"
+
+void export_Grid1D( py::module & m )
+{
+    export_Grid< Grid1D >( m, "Grid1D" );
+}
diff --git a/src/Python/pytnl/tnl/Grid2D.cpp b/src/Python/pytnl/tnl/Grid2D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..69b8683c10453d44ac415da95e83067ba4139376
--- /dev/null
+++ b/src/Python/pytnl/tnl/Grid2D.cpp
@@ -0,0 +1,9 @@
+// conversions have to be registered for each object file
+#include "../tnl_conversions.h"
+
+#include "Grid.h"
+
+void export_Grid2D( py::module & m )
+{
+    export_Grid< Grid2D >( m, "Grid2D" );
+}
diff --git a/src/Python/pytnl/tnl/Grid3D.cpp b/src/Python/pytnl/tnl/Grid3D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..15d7dd48071ac310c78d7cf057f5e7d700f72ad2
--- /dev/null
+++ b/src/Python/pytnl/tnl/Grid3D.cpp
@@ -0,0 +1,9 @@
+// conversions have to be registered for each object file
+#include "../tnl_conversions.h"
+
+#include "Grid.h"
+
+void export_Grid3D( py::module & m )
+{
+    export_Grid< Grid3D >( m, "Grid3D" );
+}
diff --git a/src/Python/pytnl/tnl/Grid_getSpaceStepsProducts.h b/src/Python/pytnl/tnl/Grid_getSpaceStepsProducts.h
new file mode 100644
index 0000000000000000000000000000000000000000..f54b9eebc5a6dde5a373a00f57733aa7a15c1885
--- /dev/null
+++ b/src/Python/pytnl/tnl/Grid_getSpaceStepsProducts.h
@@ -0,0 +1,82 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+
+#include "../typedefs.h"
+
+template< typename Grid >
+struct SpaceStepsProductsGetter
+{};
+
+template<>
+struct SpaceStepsProductsGetter< Grid1D >
+{
+    static inline
+    typename Grid1D::RealType
+    get( const Grid1D & grid, const int & xPow )
+    {
+        if( xPow == -2 ) return grid.template getSpaceStepsProducts< -2 >();
+        if( xPow == -1 ) return grid.template getSpaceStepsProducts< -1 >();
+        if( xPow ==  0 ) return grid.template getSpaceStepsProducts<  0 >();
+        if( xPow ==  1 ) return grid.template getSpaceStepsProducts<  1 >();
+        if( xPow ==  2 ) return grid.template getSpaceStepsProducts<  2 >();
+        const auto hx = grid.template getSpaceStepsProducts< 1 >();
+        return pow( hx, xPow );
+    }
+
+    template< typename PyGrid >
+    static void
+    export_getSpaceSteps( PyGrid & scope, const char* name = "getSpaceSteps" )
+    {
+        scope.def("getSpaceStepsProducts", get,   py::arg("xPow") );
+    }
+};
+
+template<>
+struct SpaceStepsProductsGetter< Grid2D >
+{
+    static inline
+    typename Grid2D::RealType
+    get( const Grid2D & grid, const int & xPow, const int & yPow = 0 )
+    {
+        if( xPow == 0 && yPow == 0 ) return grid.template getSpaceStepsProducts< 0, 0 >();
+        auto hx = grid.template getSpaceStepsProducts< 1, 0 >();
+        auto hy = grid.template getSpaceStepsProducts< 0, 1 >();
+        if( xPow != 1 ) hx = pow( hx, xPow );
+        if( yPow != 1 ) hy = pow( hy, yPow );
+        return hx * hy;
+    }
+
+    template< typename PyGrid >
+    static void
+    export_getSpaceSteps( PyGrid & scope, const char* name = "getSpaceSteps" )
+    {
+        scope.def("getSpaceStepsProducts", get,   py::arg("xPow"), py::arg("yPow")=0 );
+    }
+};
+
+template<>
+struct SpaceStepsProductsGetter< Grid3D >
+{
+    static inline
+    typename Grid3D::RealType
+    get( const Grid3D & grid, const int & xPow, const int & yPow = 0, const int & zPow = 0 )
+    {
+        if( xPow == 0 && yPow == 0 && zPow == 0 ) return grid.template getSpaceStepsProducts< 0, 0, 0 >();
+        auto hx = grid.template getSpaceStepsProducts< 1, 0, 0 >();
+        auto hy = grid.template getSpaceStepsProducts< 0, 1, 0 >();
+        auto hz = grid.template getSpaceStepsProducts< 0, 0, 1 >();
+        if( xPow != 1 ) hx = pow( hx, xPow );
+        if( yPow != 1 ) hy = pow( hy, yPow );
+        if( zPow != 1 ) hz = pow( hz, zPow );
+        return hx * hy * hz;
+    }
+
+    template< typename PyGrid >
+    static void
+    export_getSpaceSteps( PyGrid & scope, const char* name = "getSpaceSteps" )
+    {
+        scope.def("getSpaceStepsProducts", get,   py::arg("xPow"), py::arg("yPow")=0, py::arg("zPow")=0 );
+    }
+};
diff --git a/src/Python/pytnl/tnl/Mesh.cpp b/src/Python/pytnl/tnl/Mesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6a3bedeedd12e7985074173bcb82a490897e1a9a
--- /dev/null
+++ b/src/Python/pytnl/tnl/Mesh.cpp
@@ -0,0 +1,30 @@
+// conversions have to be registered for each object file
+#include "../tnl_conversions.h"
+
+#include "Mesh.h"
+#include <TNL/Meshes/Readers/VTKReader.h>
+
+void export_Meshes( py::module & m )
+{
+    export_Mesh< MeshOfEdges >( m, "MeshOfEdges" );
+    export_Mesh< MeshOfTriangles >( m, "MeshOfTriangles" );
+    export_Mesh< MeshOfTetrahedrons >( m, "MeshOfTetrahedrons" );
+
+    using Reader = TNL::Meshes::Readers::VTKReader;
+
+    py::scope entity = py::class_< Reader >( m, "VTKReader" )
+        .def(py::init<>())
+        .def("readMesh", &Reader::template readMesh< MeshOfEdges >)
+        .def("readMesh", &Reader::template readMesh< MeshOfTriangles >)
+        .def("readMesh", &Reader::template readMesh< MeshOfTetrahedrons >)
+//        .def("readMesh", []( Reader& reader, const std::string& name, MeshOfEdges & mesh ) {
+//                return reader.readMesh( name.c_str(), mesh );
+//            } )
+//        .def("readMesh", []( Reader& reader, const std::string& name, MeshOfTriangles & mesh ) {
+//                return reader.readMesh( name.c_str(), mesh );
+//            } )
+//        .def("readMesh", []( Reader& reader, const std::string& name, MeshOfTetrahedrons & mesh ) {
+//                return reader.readMesh( name.c_str(), mesh );
+//            } )
+    ;
+}
diff --git a/src/Python/pytnl/tnl/Mesh.h b/src/Python/pytnl/tnl/Mesh.h
new file mode 100644
index 0000000000000000000000000000000000000000..7ff7d61b0f24f30ca617bf17fae01333838c8c26
--- /dev/null
+++ b/src/Python/pytnl/tnl/Mesh.h
@@ -0,0 +1,144 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+
+#include "../typedefs.h"
+#include "StaticVector.h"
+#include "EntityTypes.h"
+
+#include <TNL/String.h>
+#include "../../lib_general/mesh_helpers.h"
+
+#include <type_traits>
+
+template< typename MeshEntity >
+typename MeshEntity::MeshTraitsType::GlobalIndexType
+getIndex( const MeshEntity& entity )
+{
+    return entity.getIndex();
+};
+
+struct _general {};
+struct _special : _general {};
+
+template< typename MeshEntity,
+          int Superdimension,
+          typename Scope,
+          typename = typename std::enable_if< Superdimension <= MeshEntity::MeshTraitsType::meshDimension >::type >
+//                                           && MeshEntity::template SuperentityTraits< Superdimension >::storageEnabled >::type >
+void export_getSuperentityIndex( Scope & m, _special )
+{
+    m.def("getSuperentityIndex", []( const MeshEntity& entity, const typename MeshEntity::LocalIndexType& i ) {
+                                        return entity.template getSuperentityIndex< Superdimension >( i );
+                                    }
+            );
+}
+
+template< typename MeshEntity,
+          int Superdimension,
+          typename Scope >
+void export_getSuperentityIndex( Scope &, _general )
+{
+}
+
+
+template< typename MeshEntity,
+          int Subdimension,
+          typename Scope,
+          typename = typename std::enable_if< Subdimension <= MeshEntity::MeshTraitsType::meshDimension
+                                              && (Subdimension < MeshEntity::getEntityDimension()) >::type >
+void export_getSubentityIndex( Scope & m, const char* name, _special )
+{
+    m.def(name, []( const MeshEntity& entity, const typename MeshEntity::LocalIndexType& i ) {
+                    return entity.template getSubentityIndex< Subdimension >( i );
+                }
+            );
+}
+
+template< typename MeshEntity,
+          int Superdimension,
+          typename Scope >
+void export_getSubentityIndex( Scope &, const char*, _general )
+{
+}
+
+
+template< typename MeshEntity,
+          typename Scope,
+          typename = typename std::enable_if< MeshEntity::getEntityDimension() == 0 >::type >
+void export_getPoint( Scope & scope, _special )
+{
+    scope.def("getPoint", []( const MeshEntity& entity ) {
+                            return entity.getPoint();
+                        }
+            );
+}
+
+template< typename MeshEntity,
+          typename Scope >
+void export_getPoint( Scope &, _general )
+{
+}
+
+
+template< typename MeshEntity, typename Scope >
+void export_MeshEntity( Scope & scope, const char* name )
+{
+    auto entity = py::class_< MeshEntity >( scope, name )
+        .def_static("getEntityDimension", &MeshEntity::getEntityDimension)
+        // FIXME: boost chokes on this
+//        .def("getIndex", &MeshEntity::getIndex, py::return_internal_reference<>())
+        .def("getIndex", getIndex< MeshEntity >)
+        // TODO
+    ;
+
+    export_getSuperentityIndex< MeshEntity, MeshEntity::getEntityDimension() + 1 >( entity, _special() );
+    export_getSubentityIndex< MeshEntity, 0 >( entity, "getSubvertexIndex", _special() );
+    export_getPoint< MeshEntity >( entity, _special() );
+}
+
+template< typename Mesh >
+void export_Mesh( py::module & m, const char* name )
+{
+    // there are two templates - const and non-const - take only the const
+    auto (Mesh::* getEntity_cell)(const typename Mesh::GlobalIndexType&) const = &Mesh::template getEntity<typename Mesh::Cell>;
+    auto (Mesh::* getEntity_face)(const typename Mesh::GlobalIndexType&) const = &Mesh::template getEntity<typename Mesh::Face>;
+    auto (Mesh::* getEntity_vertex)(const typename Mesh::GlobalIndexType&) const = &Mesh::template getEntity<typename Mesh::Vertex>;
+
+    export_EntityTypes(m);
+
+    auto mesh = py::class_< Mesh, TNL::Object >( m, name )
+        .def(py::init<>())
+        .def_static("getMeshDimension", &Mesh::getMeshDimension)
+        .def_static("getType",              &Mesh::getType)
+        .def("getTypeVirtual",              &Mesh::getTypeVirtual)
+        .def_static("getSerializationType", &Mesh::getSerializationType)
+        .def("getSerializationTypeVirtual", &Mesh::getSerializationTypeVirtual)
+        .def("getEntitiesCount", &mesh_getEntitiesCount< Mesh >)
+        // TODO: if combined, the return type would depend on the runtime parameter (entity)
+        .def("getEntity_cell", getEntity_cell, py::return_value_policy::reference_internal)
+        .def("getEntity_face", getEntity_face, py::return_value_policy::reference_internal)
+        .def("getEntity_vertex", getEntity_vertex, py::return_value_policy::reference_internal)
+        .def("getEntityCenter", []( const Mesh& mesh, const typename Mesh::Cell& cell ){ return getEntityCenter( mesh, cell ); } )
+        .def("getEntityCenter", []( const Mesh& mesh, const typename Mesh::Face& face ){ return getEntityCenter( mesh, face ); } )
+        .def("getEntityCenter", []( const Mesh& mesh, const typename Mesh::Vertex& vertex ){ return getEntityCenter( mesh, vertex ); } )
+        .def("getEntityMeasure", []( const Mesh& mesh, const typename Mesh::Cell& cell ){ return getEntityMeasure( mesh, cell ); } )
+        .def("getEntityMeasure", []( const Mesh& mesh, const typename Mesh::Face& face ){ return getEntityMeasure( mesh, face ); } )
+        .def("getEntityMeasure", []( const Mesh& mesh, const typename Mesh::Vertex& vertex ){ return getEntityMeasure( mesh, vertex ); } )
+        .def("isBoundaryEntity", []( const Mesh& mesh, const typename Mesh::Cell& cell ){
+                                       return mesh.template isBoundaryEntity< Mesh::Cell::getEntityDimension() >( cell.getIndex() ); } )
+        .def("isBoundaryEntity", []( const Mesh& mesh, const typename Mesh::Face& face ){
+                                       return mesh.template isBoundaryEntity< Mesh::Face::getEntityDimension() >( face.getIndex() ); } )
+        .def("isBoundaryEntity", []( const Mesh& mesh, const typename Mesh::Vertex& vertex ){
+                                        return mesh.template isBoundaryEntity< Mesh::Vertex::getEntityDimension() >( vertex.getIndex() ); } )
+        // TODO: more?
+    ;
+
+    // nested types
+    export_MeshEntity< typename Mesh::Cell >( mesh, "Cell" );
+    export_MeshEntity< typename Mesh::Face >( mesh, "Face" );
+    // avoid duplicate conversion if the type is the same
+    if( ! std::is_same< typename Mesh::Face, typename Mesh::Vertex >::value )
+        export_MeshEntity< typename Mesh::Vertex >( mesh, "Vertex" );
+}
diff --git a/src/Python/pytnl/tnl/Object.cpp b/src/Python/pytnl/tnl/Object.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cd592a2c6ac8eb0e3cd24ce09c5e27110e92f8a2
--- /dev/null
+++ b/src/Python/pytnl/tnl/Object.cpp
@@ -0,0 +1,22 @@
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+
+// conversions have to be registered for each object file
+#include "../tnl_conversions.h"
+
+#include <TNL/Object.h>
+
+void export_Object( py::module & m )
+{
+    py::class_< TNL::Object >( m, "Object" )
+        // TODO: make it abstract class in Python
+        .def("save", (bool (TNL::Object::*)(const TNL::String &) const) &TNL::Object::save)
+        .def("load", (bool (TNL::Object::*)(const TNL::String &)) &TNL::Object::load)
+        .def("boundLoad", (bool (TNL::Object::*)(const TNL::String &)) &TNL::Object::boundLoad)
+        // FIXME: why does it not work?
+//        .def("save", py::overload_cast<TNL::File>(&TNL::Object::save, py::const_))
+//        .def("load", py::overload_cast<TNL::File>(&TNL::Object::load))
+        .def("save", (bool (TNL::Object::*)(TNL::File &) const) &TNL::Object::save)
+        .def("load", (bool (TNL::Object::*)(TNL::File &)) &TNL::Object::load)
+    ;
+}
diff --git a/src/Python/pytnl/tnl/SparseMatrix.cpp b/src/Python/pytnl/tnl/SparseMatrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e6584998313fa9e3c1314c6f67b99267815cf0a8
--- /dev/null
+++ b/src/Python/pytnl/tnl/SparseMatrix.cpp
@@ -0,0 +1,29 @@
+// conversions have to be registered for each object file
+#include "../tnl_conversions.h"
+
+#include "SparseMatrix.h"
+
+#include <TNL/Matrices/CSR.h>
+#include <TNL/Matrices/Ellpack.h>
+#include <TNL/Matrices/SlicedEllpack.h>
+
+using CSR_host = TNL::Matrices::CSR< double, TNL::Devices::Host, int >;
+using CSR_cuda = TNL::Matrices::CSR< double, TNL::Devices::Cuda, int >;
+using E_host = TNL::Matrices::Ellpack< double, TNL::Devices::Host, int >;
+using E_cuda = TNL::Matrices::Ellpack< double, TNL::Devices::Cuda, int >;
+using SE_host = TNL::Matrices::SlicedEllpack< double, TNL::Devices::Host, int >;
+using SE_cuda = TNL::Matrices::SlicedEllpack< double, TNL::Devices::Cuda, int >;
+
+void export_SparseMatrices( py::module & m )
+{
+    export_Matrix< CSR_host >( m, "CSR" );
+    export_Matrix< E_host   >( m, "Ellpack" );
+    export_Matrix< SE_host  >( m, "SlicedEllpack" );
+
+    m.def("copySparseMatrix", &TNL::Matrices::copySparseMatrix< CSR_host, E_host >);
+    m.def("copySparseMatrix", &TNL::Matrices::copySparseMatrix< E_host, CSR_host >);
+    m.def("copySparseMatrix", &TNL::Matrices::copySparseMatrix< CSR_host, SE_host >);
+    m.def("copySparseMatrix", &TNL::Matrices::copySparseMatrix< SE_host, CSR_host >);
+    m.def("copySparseMatrix", &TNL::Matrices::copySparseMatrix< E_host, SE_host >);
+    m.def("copySparseMatrix", &TNL::Matrices::copySparseMatrix< SE_host, E_host >);
+}
diff --git a/src/Python/pytnl/tnl/SparseMatrix.h b/src/Python/pytnl/tnl/SparseMatrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..5d278ae8a17fc204f51b39652cc5dfe735c28dca
--- /dev/null
+++ b/src/Python/pytnl/tnl/SparseMatrix.h
@@ -0,0 +1,113 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+
+#include <TNL/String.h>
+#include <TNL/Containers/Vector.h>
+#include <TNL/Matrices/CSR.h>
+
+template< typename Matrix >
+struct SpecificExports
+{
+    template< typename Scope >
+    static void exec( Scope & s ) {}
+};
+
+template< typename Real, typename Device, typename Index >
+struct SpecificExports< TNL::Matrices::CSR< Real, Device, Index > >
+{
+    template< typename Scope >
+    static void exec( Scope & s )
+    {
+        using Matrix = TNL::Matrices::CSR< Real, Device, Index >;
+
+        s.def("getRowPointers",   py::overload_cast<>(&Matrix::getRowPointers),   py::return_value_policy::reference_internal);
+        s.def("getColumnIndexes", py::overload_cast<>(&Matrix::getColumnIndexes), py::return_value_policy::reference_internal);
+        s.def("getValues",        py::overload_cast<>(&Matrix::getValues),        py::return_value_policy::reference_internal);
+    }
+};
+
+template< typename MatrixRow >
+void export_MatrixRow( py::module & m, const char* name )
+{
+    // guard against duplicate to-Python converters for the same type
+    static bool defined = false;
+    if( ! defined ) {
+        py::class_< MatrixRow >( m, name )
+            .def("setElement", &MatrixRow::setElement)
+            .def("getElementColumn", &MatrixRow::getElementColumn, py::return_value_policy::reference_internal)
+            .def("getElementValue", &MatrixRow::getElementValue, py::return_value_policy::reference_internal)
+//            .def(py::self_ns::str(py::self_ns::self))
+        ;
+        defined = true;
+    }
+}
+
+template< typename Matrix >
+void export_Matrix( py::module & m, const char* name )
+{
+    typename Matrix::MatrixRow (Matrix::* _getRow)(typename Matrix::IndexType) = &Matrix::getRow;
+
+    using VectorType = TNL::Containers::Vector< typename Matrix::RealType, typename Matrix::DeviceType, typename Matrix::IndexType >;
+
+    auto matrix = py::class_< Matrix, TNL::Object >( m, name )
+        .def(py::init<>())
+        // overloads (defined in Object)
+        .def_static("getType",              &Matrix::getType)
+        .def("getTypeVirtual",              &Matrix::getTypeVirtual)
+        .def_static("getSerializationType", &Matrix::getSerializationType)
+        .def("getSerializationTypeVirtual", &Matrix::getSerializationTypeVirtual)
+        .def("print", &Matrix::print)
+        .def("__str__", []( Matrix & m ) {
+                std::stringstream ss;
+                ss << m;
+                return ss.str();
+            } )
+
+        // Matrix
+        .def("setDimensions",           &Matrix::setDimensions)
+        .def("setCompressedRowLengths", &Matrix::setCompressedRowLengths)
+        .def("getRowLength",            &Matrix::getRowLength)
+        .def("getCompressedRowLengths", &Matrix::getCompressedRowLengths)
+        // TODO: export for more types
+        .def("setLike",                 &Matrix::template setLike< typename Matrix::RealType, typename Matrix::DeviceType, typename Matrix::IndexType >)
+        .def("getNumberOfMatrixElements", &Matrix::getNumberOfMatrixElements)
+        .def("getNumberOfNonzeroMatrixElements", &Matrix::getNumberOfNonzeroMatrixElements)
+        .def("reset",                   &Matrix::reset)
+        .def("getRows",                 &Matrix::getRows)
+        .def("getColumns",              &Matrix::getColumns)
+        .def("setElement",              &Matrix::setElement)
+        .def("addElement",              &Matrix::addElement)
+        // setRow and addRow operate on pointers
+        //.def("setRow",                  &Matrix::setRow)
+        //.def("addRow",                  &Matrix::addRow)
+        .def("getElement",              &Matrix::getElement)
+        // TODO: operator== and operator!= are general and very slow
+        
+        // Sparse
+        .def("getMaxRowLength",     &Matrix::getMaxRowLength)
+        .def("getPaddingIndex",     &Matrix::getPaddingIndex)
+        // TODO: this one is empty in the C++ code
+//        .def("printStructure",      &Matrix::printStructure)
+
+        // specific to each format, but with common interface
+        .def("getRow",              _getRow)
+        // TODO: export for more types
+        .def("rowVectorProduct",    &Matrix::template rowVectorProduct< VectorType >)
+        .def("vectorProduct",       &Matrix::template vectorProduct< VectorType, VectorType >)
+        // TODO: these two don't work
+        //.def("addMatrix",           &Matrix::addMatrix)
+        //.def("getTransposition",    &Matrix::getTransposition)
+        .def("performSORIteration", &Matrix::template performSORIteration< VectorType >)
+//        .def("assign",              &Matrix::operator=)
+        .def("assign", []( Matrix& matrix, const Matrix& other ) -> Matrix& {
+                return matrix = other;
+            })
+    ;
+
+    // export format-specific methods
+    SpecificExports< Matrix >::exec( matrix );
+
+    export_MatrixRow< typename Matrix::MatrixRow >( m, "MatrixRow" );
+}
diff --git a/src/Python/pytnl/tnl/StaticVector.h b/src/Python/pytnl/tnl/StaticVector.h
new file mode 100644
index 0000000000000000000000000000000000000000..53197af56195fa0c22dceba6e05bcc7e3b6b9319
--- /dev/null
+++ b/src/Python/pytnl/tnl/StaticVector.h
@@ -0,0 +1,43 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+#include <pybind11/operators.h>
+namespace py = pybind11;
+
+#include "../tnl_indexing.h"
+
+// needed for discovery of operator<< for tnlStaticArray
+#include <TNL/Containers/StaticArray.h>
+
+template< typename VectorType, typename Scope >
+void export_StaticVector( Scope & scope, const char* name )
+{
+    using RealType = typename VectorType::RealType;
+
+    auto vector = py::class_<VectorType>(scope, name)
+        .def(py::init< typename VectorType::RealType >())
+        .def(py::init< VectorType >())
+        .def_static("getType", &VectorType::getType)
+        .def("getSize", &VectorType::getSize)
+        .def("assign", &VectorType::operator=)
+        .def(py::self == py::self)
+        .def(py::self != py::self)
+        .def("setValue", &VectorType::setValue)
+        // TODO: pybind11
+        // explicit namespace resolution is necessary, see http://stackoverflow.com/a/3084341/4180822
+//        .def(py::self_ns::str(py::self))
+        .def(py::self += py::self)
+        .def(py::self -= py::self)
+        .def(py::self *= typename VectorType::RealType())
+        .def(py::self + py::self)
+        .def(py::self - py::self)
+        .def(py::self * typename VectorType::RealType())
+        .def(py::self * py::self)
+        .def(py::self < py::self)
+        .def(py::self > py::self)
+        .def(py::self <= py::self)
+        .def(py::self >= py::self)
+    ;
+
+    tnl_indexing< VectorType >( vector );
+}
diff --git a/src/Python/pytnl/tnl/String.cpp b/src/Python/pytnl/tnl/String.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..309a6fe25fdafae45c672cd321ad851a55f57056
--- /dev/null
+++ b/src/Python/pytnl/tnl/String.cpp
@@ -0,0 +1,41 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/operators.h>
+namespace py = pybind11;
+
+// conversions have to be registered for each object file
+#include "../tnl_conversions.h"
+
+#include <TNL/String.h>
+#include <TNL/File.h>
+
+void export_String( py::module & m )
+{
+    // function pointers for overloaded methods
+//    const char* (TNL::String::* _String_getString)() const = &TNL::String::getString;
+    bool (TNL::String::* _String_save)(TNL::File &) const = &TNL::String::save;
+    bool (TNL::String::* _String_load)(TNL::File &)       = &TNL::String::load;
+
+    py::class_<TNL::String>(m, "String")
+        .def(py::init<const char*>())
+        .def(py::init<const char*, int>())
+        .def(py::init<const char*, int, int>())
+        .def(py::init<int>())
+        .def(py::init<double>())
+        .def_static("getType", &TNL::String::getType)
+        // __str__ (uses operator<<)
+        // explicit namespace resolution is necessary, see http://stackoverflow.com/a/3084341/4180822
+//        .def(py::self_ns::str(py::self_ns::self))
+        // TODO: operator[]
+//        .def(vector_indexing_suite<TNL::String>())
+        .def(py::self + py::self)
+        .def(py::self += py::self)
+        .def(py::self == py::self)
+        .def(py::self != py::self)
+        .def("getLength", &TNL::String::getLength)
+        .def("__len__", &TNL::String::getLength)
+        // FIXME
+//        .def("replace", &TNL::String::replace)
+        .def("save", _String_save)
+        .def("load", _String_load)
+    ;
+}
diff --git a/src/Python/pytnl/tnl/Vector.h b/src/Python/pytnl/tnl/Vector.h
new file mode 100644
index 0000000000000000000000000000000000000000..7f4232e1bc48954750bacf7c47f4da5b23301dee
--- /dev/null
+++ b/src/Python/pytnl/tnl/Vector.h
@@ -0,0 +1,52 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+#include <pybind11/operators.h>
+namespace py = pybind11;
+
+#include <TNL/Containers/Vector.h>
+
+template< typename ArrayType, typename VectorType >
+void export_Vector(py::module & m, const char* name)
+{
+    // function pointers for overloaded methods
+    void (VectorType::* _addElement1)(const typename VectorType::IndexType,
+                                      const typename VectorType::RealType &)
+        = &VectorType::addElement;
+    void (VectorType::* _addElement2)(const typename VectorType::IndexType,
+                                      const typename VectorType::RealType &,
+                                      const typename VectorType::RealType &)
+        = &VectorType::addElement;
+
+    py::class_<VectorType, ArrayType>(m, name)
+        .def(py::init<>())
+        .def(py::init<int>())
+        .def_static("getType",              &VectorType::getType)
+        .def("getTypeVirtual",              &VectorType::getTypeVirtual)
+        .def_static("getSerializationType", &VectorType::getSerializationType)
+        .def("getSerializationTypeVirtual", &VectorType::getSerializationTypeVirtual)
+        .def("addElement", _addElement1)
+        .def("addElement", _addElement2)
+        .def(py::self == py::self)
+        .def(py::self != py::self)
+        .def(py::self += py::self)
+        .def(py::self -= py::self)
+        .def(py::self *= typename VectorType::RealType())
+        .def(py::self /= typename VectorType::RealType())
+        .def("max", &VectorType::max)
+        .def("min", &VectorType::min)
+        .def("absMax", &VectorType::absMax)
+        .def("absMin", &VectorType::absMin)
+        .def("lpNorm", &VectorType::template lpNorm<double, double>)
+        .def("sum", &VectorType::template sum<double>)
+        .def("differenceMax", &VectorType::template differenceMax<VectorType>)
+        .def("differenceMin", &VectorType::template differenceMin<VectorType>)
+        .def("differenceAbsMax", &VectorType::template differenceAbsMax<VectorType>)
+        .def("differenceAbsMin", &VectorType::template differenceAbsMin<VectorType>)
+        .def("differenceLpNorm", &VectorType::template differenceLpNorm<double, VectorType, double>)
+        .def("differenceSum", &VectorType::template differenceSum<double, VectorType>)
+        .def("scalarProduct", &VectorType::template scalarProduct<VectorType>)
+        .def("addVector", &VectorType::template addVector<VectorType>)
+        .def("addVectors", &VectorType::template addVectors<VectorType>)
+    ;
+}
diff --git a/src/Python/pytnl/tnl/tnl.cpp b/src/Python/pytnl/tnl/tnl.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c71a67fa4ec44f69724d2bd2968b93f5cef7c169
--- /dev/null
+++ b/src/Python/pytnl/tnl/tnl.cpp
@@ -0,0 +1,47 @@
+#include "../exceptions.h"
+#include "../typedefs.h"
+
+// conversions have to be registered for each object file
+#include "../tnl_conversions.h"
+
+#include "Array.h"
+#include "Vector.h"
+
+// external functions
+void export_Object( py::module & m );
+void export_String( py::module & m );
+void export_Grid1D( py::module & m );
+void export_Grid2D( py::module & m );
+void export_Grid3D( py::module & m );
+//void export_Meshes( py::module & m );
+void export_SparseMatrices( py::module & m );
+
+template< typename T >
+using _array = TNL::Containers::Array< T, TNL::Devices::Host, IndexType >;
+
+template< typename T >
+using _vector = TNL::Containers::Vector< T, TNL::Devices::Host, IndexType >;
+
+// Python module definition
+PYBIND11_MODULE(tnl, m)
+{
+    register_exceptions(m);
+
+    export_Object(m);
+    // TODO: TNL::File
+    export_String(m);
+
+    export_Array< _array<double> >(m, "Array");
+    export_Vector< _array<double>, _vector<double> >(m, "Vector");
+    export_Array< _array<int> >(m, "Array_int");
+    export_Vector< _array<int>, _vector<int> >(m, "Vector_int");
+    export_Array< _array<bool> >(m, "Array_bool");
+
+    export_Grid1D(m);
+    export_Grid2D(m);
+    export_Grid3D(m);
+
+//    export_Meshes(m);
+
+    export_SparseMatrices(m);
+}
diff --git a/src/Python/pytnl/tnl_conversions.h b/src/Python/pytnl/tnl_conversions.h
new file mode 100644
index 0000000000000000000000000000000000000000..602d1cffd1ccd32660b72630c023f2d913f1da19
--- /dev/null
+++ b/src/Python/pytnl/tnl_conversions.h
@@ -0,0 +1,3 @@
+// conversion has to be registered for each object file
+#include "tnl_str_conversion.h"
+#include "tnl_tuple_conversion.h"
diff --git a/src/Python/pytnl/tnl_indexing.h b/src/Python/pytnl/tnl_indexing.h
new file mode 100644
index 0000000000000000000000000000000000000000..29663d66b553d8095ca0a67f980b346fac04e712
--- /dev/null
+++ b/src/Python/pytnl/tnl_indexing.h
@@ -0,0 +1,81 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+
+#include "RawIterator.h"
+
+template< typename Array, typename Scope >
+void tnl_indexing( Scope & scope )
+{
+    using Index = typename Array::IndexType;
+    using Element = typename Array::ElementType;
+
+    scope.def("__len__", &Array::getSize);
+
+    scope.def("__iter__",
+        []( Array& array ) {
+            return py::make_iterator(
+                        RawIterator<Element>(array.getData()),
+                        RawIterator<Element>(array.getData() + array.getSize()) );
+        },
+        py::keep_alive<0, 1>()  // keep array alive while iterator is used
+    );
+
+    scope.def("__getitem__",
+        [](Array &a, Index i) {
+            if (i >= a.getSize())
+                throw py::index_error();
+            return a[i];
+        }
+    );
+
+    scope.def("__setitem__",
+        [](Array &a, Index i, const Element& e) {
+            if (i >= a.getSize())
+                throw py::index_error();
+            a[i] = e;
+        }
+    );
+}
+
+template< typename Array, typename Scope >
+void tnl_slice_indexing( Scope & scope )
+{
+    /// Slicing protocol
+    scope.def("__getitem__",
+        [](const Array& a, py::slice slice) -> Array* {
+            size_t start, stop, step, slicelength;
+
+            if (!slice.compute(a.getSize(), &start, &stop, &step, &slicelength))
+                throw py::error_already_set();
+
+            Array* seq = new Array();
+            seq->setSize(slicelength);
+
+            for (size_t i = 0; i < slicelength; ++i) {
+                seq->operator[](i) = a[start];
+                start += step;
+            }
+            return seq;
+        },
+        "Retrieve list elements using a slice object"
+    );
+
+    scope.def("__setitem__",
+        [](Array& a, py::slice slice,  const Array& value) {
+            size_t start, stop, step, slicelength;
+            if (!slice.compute(a.getSize(), &start, &stop, &step, &slicelength))
+                throw py::error_already_set();
+
+            if (slicelength != value.getSize())
+                throw std::runtime_error("Left and right hand size of slice assignment have different sizes!");
+
+            for (size_t i = 0; i < slicelength; ++i) {
+                a[start] = value[i];
+                start += step;
+            }
+        },
+        "Assign list elements using a slice object"
+    );
+}
diff --git a/src/Python/pytnl/tnl_str_conversion.h b/src/Python/pytnl/tnl_str_conversion.h
new file mode 100644
index 0000000000000000000000000000000000000000..b62886e58a201d30b8d52abba5e3306d203bbd31
--- /dev/null
+++ b/src/Python/pytnl/tnl_str_conversion.h
@@ -0,0 +1,51 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+
+#include <TNL/String.h>
+
+namespace pybind11 { namespace detail {
+
+    template <>
+    struct type_caster<TNL::String>
+    {
+        using StdStringCaster = type_caster<std::string>;
+        StdStringCaster _caster;
+
+    public:
+        /**
+         * This macro establishes the name 'TNL::String' in
+         * function signatures and declares a local variable
+         * 'value' of type TNL::String
+         */
+        PYBIND11_TYPE_CASTER(TNL::String, _("TNL::String"));
+
+        /**
+         * Conversion part 1 (Python->C++): convert a PyObject into a TNL::String
+         * instance or return false upon failure. The second argument
+         * indicates whether implicit conversions should be applied.
+         */
+        bool load(handle src, bool implicit)
+        {
+            if( ! _caster.load(src, implicit) )
+                return false;
+            const std::string& str = (std::string&) _caster;
+            value = TNL::String(str.c_str());
+            return true;
+        }
+
+        /**
+         * Conversion part 2 (C++ -> Python): convert an TNL::String instance into
+         * a Python object. The second and third arguments are used to
+         * indicate the return value policy and parent object (for
+         * ``return_value_policy::reference_internal``) and are generally
+         * ignored by implicit casters.
+         */
+        static handle cast(TNL::String src, return_value_policy policy, handle parent)
+        {
+            return StdStringCaster::cast( src.getString(), policy, parent );
+        }
+    };
+
+}} // namespace pybind11::detail
diff --git a/src/Python/pytnl/tnl_tuple_conversion.h b/src/Python/pytnl/tnl_tuple_conversion.h
new file mode 100644
index 0000000000000000000000000000000000000000..f4ca8045ea18bfd7e9d40a75d3201533628e662f
--- /dev/null
+++ b/src/Python/pytnl/tnl_tuple_conversion.h
@@ -0,0 +1,79 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+namespace py = pybind11;
+
+#include <TNL/Containers/StaticArray.h>
+#include <TNL/Containers/StaticVector.h>
+
+namespace pybind11 { namespace detail {
+
+    template< typename ArrayType >
+    struct _tnl_tuple_caster
+    {
+        using Value = typename std::remove_reference< decltype(ArrayType()[0]) >::type;
+        using StdArray = std::array< Value, ArrayType::size >;
+        using StdArrayCaster = type_caster< StdArray >;
+//        StdArrayCaster _caster;
+		using value_conv = make_caster<Value>;
+
+    public:
+//        PYBIND11_TYPE_CASTER(ArrayType, StdArrayCaster::name);
+        PYBIND11_TYPE_CASTER(ArrayType, _("Tuple[") + value_conv::name + _<false>(_(""), _("[") + _<ArrayType::size>() + _("]")) + _("]"));
+
+        /**
+         * Conversion part 1 (Python -> C++): convert a PyObject into an ArrayType
+         * instance or return false upon failure. The second argument indicates
+         * whether implicit conversions should be applied.
+         */
+        bool load(handle src, bool implicit)
+        {
+            // we don't use StdArrayCaster here because we want to convert Python tuples, not lists
+//            if( ! _caster.load(src, implicit) )
+//                return false;
+//            const StdArray& arr = (StdArray&) _caster;
+//            for( int i = 0; i < ArrayType::size; i++ )
+//                value[ i ] = arr[ i ];
+//            return true;
+
+			if (!isinstance<tuple>(src))
+			    return false;
+			auto t = reinterpret_borrow<tuple>(src);
+			if (t.size() != ArrayType::size)
+			    return false;
+			size_t ctr = 0;
+			for (auto it : t) {
+			    value_conv conv;
+			    if (!conv.load(it, implicit))
+			        return false;
+			    value[ctr++] = cast_op<Value &&>(std::move(conv));
+			}
+			return true;
+        }
+
+        /**
+         * Conversion part 2 (C++ -> Python): convert an ArrayType instance into
+         * a Python object. The second and third arguments are used to
+         * indicate the return value policy and parent object (for
+         * ``return_value_policy::reference_internal``) and are generally
+         * ignored by implicit casters.
+         */
+        static handle cast(const ArrayType& src, return_value_policy policy, handle parent)
+        {
+            StdArray arr;
+            for( int i = 0; i < ArrayType::size; i++ )
+                arr[ i ] = src[ i ];
+            return StdArrayCaster::cast( arr, policy, parent );
+        }
+    };
+
+    template< typename T, int Size >
+    struct type_caster< TNL::Containers::StaticArray< Size, T > >
+        : _tnl_tuple_caster< TNL::Containers::StaticArray< Size, T > > {};
+
+    template< typename T, int Size >
+    struct type_caster< TNL::Containers::StaticVector< Size, T > >
+        : _tnl_tuple_caster< TNL::Containers::StaticVector< Size, T > > {};
+
+}} // namespace pybind11::detail
diff --git a/src/Python/pytnl/typedefs.h b/src/Python/pytnl/typedefs.h
new file mode 100644
index 0000000000000000000000000000000000000000..6644a8351dd601e049add5c6bf4a28c62ac3a904
--- /dev/null
+++ b/src/Python/pytnl/typedefs.h
@@ -0,0 +1,42 @@
+#pragma once
+
+#include <TNL/Meshes/Grid.h>
+//#include <TNL/Meshes/Mesh.h>
+//#include <TNL/Meshes/DefaultConfig.h>
+//#include <TNL/Meshes/Topologies/Edge.h>
+//#include <TNL/Meshes/Topologies/Triangle.h>
+//#include <TNL/Meshes/Topologies/Tetrahedron.h>
+
+using RealType = double;
+using DeviceType = TNL::Devices::Host;
+using IndexType = int;
+
+using Grid1D = TNL::Meshes::Grid<1, RealType, DeviceType, IndexType>;
+using Grid2D = TNL::Meshes::Grid<2, RealType, DeviceType, IndexType>;
+using Grid3D = TNL::Meshes::Grid<3, RealType, DeviceType, IndexType>;
+
+//using LocalIndexType = short int;
+//using EdgeTopology = TNL::Meshes::Topologies::Edge;
+//using TriangleTopology = TNL::Meshes::Topologies::Triangle;
+//using TetrahedronTopology = TNL::Meshes::Topologies::Tetrahedron;
+//using MeshOfEdges = TNL::Meshes::Mesh< TNL::Meshes::DefaultConfig<
+//                            EdgeTopology,
+//                            EdgeTopology::dimension,
+//                            RealType,
+//                            IndexType,
+//                            LocalIndexType,
+//                            IndexType > >;
+//using MeshOfTriangles = TNL::Meshes::Mesh< TNL::Meshes::DefaultConfig<
+//                            TriangleTopology,
+//                            TriangleTopology::dimension,
+//                            RealType,
+//                            IndexType,
+//                            LocalIndexType,
+//                            IndexType > >;
+//using MeshOfTetrahedrons = TNL::Meshes::Mesh< TNL::Meshes::DefaultConfig<
+//                            TetrahedronTopology,
+//                            TetrahedronTopology::dimension,
+//                            RealType,
+//                            IndexType,
+//                            LocalIndexType,
+//                            IndexType > >;