diff --git a/cxxfem/__init__.py b/cxxfem/__init__.py
index 3a3577e00eb47acd97f625cc658016a766d30094..0e5c73c937259b88e7603c30e3abc9b93d779a4a 100644
--- a/cxxfem/__init__.py
+++ b/cxxfem/__init__.py
@@ -2,3 +2,4 @@
 # fem MODULE initialization file
 
 from femi import *
+from .utils import *
diff --git a/cxxfem/src/femSolver.cpp b/cxxfem/src/femSolver.cpp
index ec1058bb7f8ffec5816d563a9b0bdaf5737c9d60..384a1b384b666dcee0f1c5d36b06cb0bbbab847b 100644
--- a/cxxfem/src/femSolver.cpp
+++ b/cxxfem/src/femSolver.cpp
@@ -274,8 +274,8 @@ void Solver::fill_rhs()
         Kijv[i].resize(partition->parts.size());
         for (size_t j = 0; j < partition->parts.size(); ++j)
         {
-            std::cout << "expecting " << ntriplets[i][j] 
-                      << " triplets for submatrix " << i << "," << j << '\n';
+            // std::cout << "expecting " << ntriplets[i][j] 
+            //           << " triplets for submatrix " << i << "," << j << '\n';
             Kijv[i][j]->reserve(ntriplets[i][j]);
         }
     }
@@ -763,7 +763,7 @@ void Solver::solve_static()
 
     // energy: work of ext forces (using their "final" value - thus twice the int energy)
     double workF = (*f[MP]).dot(*X[MP]) +  (*f[MF]).dot(*X[MF]);
-    std::cout << "workF = " << workF << '\n';
+    //std::cout << "workF = " << workF << '\n';
     int_energy = workF/2;
 
     // free vectors of tmp_rhs
@@ -805,7 +805,7 @@ void Solver::compute_sig_epl(std::vector<size_t> &etags_all,
                     }
             }
         }
-    std::cout << "\texpecting " << nelems << " elements\n";
+    // std::cout << "\texpecting " << nelems << " elements\n";
     etags_all.reserve(nelems);
     epsilon_all.reserve(nelems);
     sigma_all.reserve(nelems);
diff --git a/cxxfem/utils.py b/cxxfem/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..c3ab0dd7768f7880d48853c844994ba1cec0ae02
--- /dev/null
+++ b/cxxfem/utils.py
@@ -0,0 +1,49 @@
+# -*- coding: utf-8 -*-
+# utilities
+
+# Duplicate console output to file
+
+import sys
+
+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):
+        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):
+        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
diff --git a/fossils.py b/fossils.py
index 1e81d706ebfc7405e3c5eb016e3538705971a239..3b0546f6dfa98a181a9a8b2e796f3787bc82b9e0 100644
--- a/fossils.py
+++ b/fossils.py
@@ -2,7 +2,8 @@
 # -*- coding: utf-8 -*-
 # runs a test as if it was installed
 
-import os, platform
+import os
+import platform
 
 # if 'Windows' in platform.uname():
 #     # remove tbb from the PATH
@@ -11,51 +12,8 @@ import os, platform
 #     corr_path = [ p for p in paths if not 'tbb' in p]
 #     os.environ['PATH'] = ';'.join(corr_path)
 #     for p in corr_path:
-#         print(p)    
+#         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
@@ -74,27 +32,20 @@ if __name__ == "__main__":
     sys.path.append(os.path.join(thisdir, 'cxxfem',
                     'build', 'bin', 'Release'))  # msvc
 
+    # redirect C++ streams to Python
+    import cxxfem
+    redirect = cxxfem.StdOutErr2Py()
 
-    args = parseargs()
+    args = cxxfem.parseargs()
+    # set number of threads from the "-k" argument
     os.environ['OMP_NUM_THREADS'] = str(args.k)
 
-
     # display env variables
     try:
         print(f"OMP_NUM_THREADS={os.environ['OMP_NUM_THREADS']}")
     except:
         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()
-
     # run test file
     testname = os.path.abspath(args.file)
     testname = os.path.normcase(testname)  # F:/ => f:/ on Windows
@@ -108,11 +59,7 @@ if __name__ == "__main__":
         os.makedirs(wdir)
     os.chdir(wdir)
 
-    # redirect C++ streams to Python
-    import cxxfem
-    redirect = cxxfem.StdOutErr2Py()
-
-    tee = Tee('stdout.txt') # split streams
+    tee = cxxfem.Tee('stdout.txt')  # split streams
 
     if ext == '.py':
         # run python script
diff --git a/models/bonemodel2.py b/models/bonemodel2.py
index 85c5d596480f749827fdb9587ecec6b3dd014343..62b2fcad8934297c984ee52ad08f8e4f940e1310 100644
--- a/models/bonemodel2.py
+++ b/models/bonemodel2.py
@@ -75,9 +75,12 @@ def parms(d={}):
 
 def solve(p={}):
 
+    import time
+    start_time = time.time()
+
     p = parms(p)
 
-    import json
+    # import json
     # print(f'parameters={json.dumps(p, indent=4, sort_keys=True)}')
 
     pbl = fem.Problem()
@@ -206,7 +209,13 @@ def solve(p={}):
         print(f'\t{name}_Fy = {v[1]:.3f} N')
         print(f'\t{name}_Fz = {v[2]:.3f} N')
 
-    post.view()
+    end_time = time.time()
+    elapsed_time = end_time - start_time
+    print(f'Total Execution time: {elapsed_time:.1f}s')
+
+    args = fem.parseargs()
+    if not args.nogui:
+        post.view()
 
 
 # ------------------------------------------------------------------------------