diff --git a/cxxfem/_src/CMakeLists.txt b/cxxfem/_src/CMakeLists.txt
index f3549cc40fadd9fdb29ddb5bf1dda8023b028f5f..b8f16d368c3e917db24b3d177085cac9af30961f 100644
--- a/cxxfem/_src/CMakeLists.txt
+++ b/cxxfem/_src/CMakeLists.txt
@@ -12,6 +12,7 @@ SET_PROPERTY(TARGET femi PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
 # MACRO_DebugPostfix(femi)
 
 TARGET_INCLUDE_DIRECTORIES(femi PRIVATE 
-    ${PROJECT_SOURCE_DIR}/src
-    ${Python3_INCLUDE_DIRS})
+${PROJECT_SOURCE_DIR}/src
+${PROJECT_SOURCE_DIR}/_src
+${Python3_INCLUDE_DIRS})
 TARGET_LINK_LIBRARIES(femi PRIVATE fem ${Python3_LIBRARIES})
diff --git a/cxxfem/_src/femi.i b/cxxfem/_src/femi.i
index 41956c6d8aee6d5576b4b3205edbfaad6e8b17b0..4cf565837dd454088d3fab725ba92256c5a5a558 100644
--- a/cxxfem/_src/femi.i
+++ b/cxxfem/_src/femi.i
@@ -18,6 +18,8 @@ threads="1"
 #include "femMedium.h"
 #include "femExtractor.h"
 
+#include "wCppBuf2Py.h"
+
 #include <string>
 // #include <sstream>
 
@@ -72,6 +74,10 @@ namespace std {
    } 
 } 
 
+%warnfilter(401); //Nothing known about base class 'std::basic_streambuf< char >'
+%include "wCppBuf2Py.h"
+%warnfilter(+401);
+
 // ----------- FEM CLASSES ---------------
 %include "fem.h"
 %include "femDirichlet.h"
diff --git a/cxxfem/_src/wCppBuf2Py.cpp b/cxxfem/_src/wCppBuf2Py.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d4d4d0ee73a160f5ebac52cb9d985f8e371391f7
--- /dev/null
+++ b/cxxfem/_src/wCppBuf2Py.cpp
@@ -0,0 +1,102 @@
+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <Python.h>
+#include "wCppBuf2Py.h"
+#include <stdio.h>
+#include <algorithm>
+
+CppBuf2Py::CppBuf2Py(std::ostream &o1, bool err) : stream(o1), errstream(err)
+{
+    oldbuf = stream.rdbuf();
+    stream.rdbuf(this);
+    verb = false;
+}
+
+CppBuf2Py::~CppBuf2Py()
+{
+    stream.rdbuf(oldbuf); // restore the ostream
+}
+
+std::streamsize
+CppBuf2Py::xsputn(const char *s, std::streamsize n)
+{
+    if (verb)
+        printf("CppBuf2Py::xsputn(%ld)\n", static_cast<long>(n));
+    // PySys_WriteStdout truncates to 1000 chars
+    static const std::streamsize MAXSIZE = 1000;
+    std::streamsize written = std::min(n, MAXSIZE);
+    if (verb)
+        printf("written=%ld\n", static_cast<long>(written));
+    std::string str(s, n);
+    if (verb)
+        printf("string=%s\n", str.c_str());
+
+    PyGILState_STATE gstate = PyGILState_Ensure();
+    if (errstream)
+        PySys_WriteStderr("%s", str.c_str());
+    else
+        PySys_WriteStdout("%s", str.c_str());
+    PyGILState_Release(gstate);
+
+    return written;
+}
+
+int CppBuf2Py::overflow(int c)
+{
+    if (verb)
+        printf("CppBuf2Py::overflow(%d)\n", c);
+
+    if (c != EOF)
+    {
+        PyGILState_STATE gstate = PyGILState_Ensure();
+        if (errstream)
+            PySys_WriteStderr("%c", c);
+        else
+            PySys_WriteStdout("%c", c);
+        PyGILState_Release(gstate);
+    }
+    return c;
+}
+
+void CppBuf2Py::test()
+{
+    CppBuf2Py buf(std::cout);
+    for (int i = 0; i < 1000; ++i)
+        std::cout << "(" << i << ")"
+                  << " this is a test - ";
+}
+
+// =============================================================================
+
+StdOutErr2Py::StdOutErr2Py() : out(std::cout), err(std::cerr, true)
+{
+}
+
+void StdOutErr2Py::test()
+{
+    StdOutErr2Py redirector;
+
+    for (int i = 0; i < 50; ++i)
+    {
+        std::cout << "(" << i << ")"
+                  << " this is a test - ";
+        std::cerr << "(" << i << ")"
+                  << " error - ";
+    }
+    std::cout << '\n';
+    std::cerr << '\n';
+}
diff --git a/cxxfem/_src/wCppBuf2Py.h b/cxxfem/_src/wCppBuf2Py.h
new file mode 100644
index 0000000000000000000000000000000000000000..6ac121c5d449175e5f6bd747e627336ccb840089
--- /dev/null
+++ b/cxxfem/_src/wCppBuf2Py.h
@@ -0,0 +1,53 @@
+#ifndef CPPBUF2PY_H
+#define CPPBUF2PY_H
+
+#include "fem.h"
+#include <iostream>
+
+// http://bo-peng.blogspot.be/2004/10/how-to-re-direct-cout-to-python_05.html
+// https://docs.python.org/2/c-api/sys.html#system-functions
+
+/**
+ * @brief redirect a std::stream to python
+ */
+#ifndef SWIG
+class CppBuf2Py : public std::basic_streambuf<char>
+{
+public:
+    typedef std::char_traits<char> traits_type;
+    typedef traits_type::int_type int_type;
+
+private:
+    std::ostream &stream;
+    std::streambuf *oldbuf;
+    bool verb;
+    bool errstream;
+
+public:
+    CppBuf2Py(std::ostream &o1, bool err = false);
+    ~CppBuf2Py();
+
+    static void test();
+
+protected:
+    virtual std::streamsize xsputn(const char *s, std::streamsize n);
+    virtual int_type overflow(int c);
+};
+#endif //SWIG
+
+/**
+ * @brief redirect std::cout/std::cerr to python sys.stdout/sys.stderr
+ */
+
+class StdOutErr2Py
+{
+    CppBuf2Py out;
+    CppBuf2Py err;
+
+public:
+    StdOutErr2Py();
+
+    static void test();
+};
+
+#endif //CPPBUF2PY_H
diff --git a/run.py b/run.py
index 56af432e71f62d254f464400af03e60d2ebfb57f..1e81d706ebfc7405e3c5eb016e3538705971a239 100644
--- a/run.py
+++ b/run.py
@@ -13,6 +13,50 @@ import os, platform
 #     for p in corr_path:
 #         print(p)    
 
+# Duplicate console output to file
+
+class DupStream:
+    def __init__(self, stream1, stream2):
+        self.stream1 = stream1
+        self.stream2 = stream2
+
+    def write(self, data):
+        self.stream1.write(data)
+        self.stream2.write(data)
+
+    def flush(self):
+        self.stream1.flush()
+        self.stream2.flush()
+
+class Tee:
+    def __init__(self, name):
+        import sys
+        self.file = open(name, 'w')
+        self.stdoutbak = sys.stdout
+        self.stderrbak = sys.stderr
+        sys.stdout = DupStream(sys.stdout, self.file)
+        sys.stderr = DupStream(sys.stderr, self.file)
+
+    def __del__(self):
+        import sys
+        sys.stdout = self.stdoutbak
+        sys.stderr = self.stderrbak
+        self.file.close()
+
+
+def parseargs():
+    """
+    parses command line arguments
+    """
+    import argparse
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-v", "--verb", help="increase output verbosity", action="count", default=0)
+    parser.add_argument("--nogui", help="disable any graphical output", action="store_true")
+    parser.add_argument("-k", help="nb of threads", type=int, default=1)
+    parser.add_argument('file', help='python files')
+    args = parser.parse_args()
+    return args
+
 if __name__ == "__main__":
     import sys
     import os
@@ -30,6 +74,11 @@ if __name__ == "__main__":
     sys.path.append(os.path.join(thisdir, 'cxxfem',
                     'build', 'bin', 'Release'))  # msvc
 
+
+    args = parseargs()
+    os.environ['OMP_NUM_THREADS'] = str(args.k)
+
+
     # display env variables
     try:
         print(f"OMP_NUM_THREADS={os.environ['OMP_NUM_THREADS']}")
@@ -37,29 +86,35 @@ if __name__ == "__main__":
         print(f"OMP_NUM_THREADS=[not set]")
 
     # parse args
-    import argparse
-    parser = argparse.ArgumentParser()
-    parser.add_argument(
-        "-v", "--verb", help="increase output verbosity", action="count", default=0)
-    parser.add_argument("--nogui", help="disable any graphical output",
-                        action="store_true")
-    parser.add_argument('file', nargs='*', help='python files')
-    args = parser.parse_args()
+    # import argparse
+    # parser = argparse.ArgumentParser()
+    # parser.add_argument(
+    #     "-v", "--verb", help="increase output verbosity", action="count", default=0)
+    # parser.add_argument("--nogui", help="disable any graphical output",
+    #                     action="store_true")
+    # parser.add_argument('file', nargs='*', help='python files')
+    # args = parser.parse_args()
+
+    # run test file
+    testname = os.path.abspath(args.file)
+    testname = os.path.normcase(testname)  # F:/ => f:/ on Windows
+    # create workspace
+    common = os.path.commonprefix((testname, thisdir + os.sep))
+    resdir = testname[len(common):].replace(os.sep, "_")
+    resdir, ext = os.path.splitext(resdir)
+    wdir = os.path.join('workspace', resdir)
+    print('workspace=', wdir)
+    if not os.path.isdir(wdir):
+        os.makedirs(wdir)
+    os.chdir(wdir)
+
+    # redirect C++ streams to Python
+    import cxxfem
+    redirect = cxxfem.StdOutErr2Py()
+
+    tee = Tee('stdout.txt') # split streams
 
-    # run all tests sequentially
-    for testname in args.file:
-        testname = os.path.abspath(testname)
-        testname = os.path.normcase(testname)  # F:/ => f:/ on Windows
-        # create workspace
-        common = os.path.commonprefix((testname, thisdir + os.sep))
-        resdir = testname[len(common):].replace(os.sep, "_")
-        resdir, ext = os.path.splitext(resdir)
-        wdir = os.path.join('workspace', resdir)
-        print('workspace=', wdir)
-        if not os.path.isdir(wdir):
-            os.makedirs(wdir)
-        os.chdir(wdir)
-        if ext == '.py':
-            # python script
-            __file__ = testname
-            exec(open(testname, encoding='utf-8').read())
+    if ext == '.py':
+        # run python script
+        __file__ = testname
+        exec(open(testname, encoding='utf-8').read())