diff --git a/cxxfem/src/femProblem.cpp b/cxxfem/src/femProblem.cpp
index 6d120a18cf008cc58a34436b8583b1ee724f5695..124f4cd9fde71d5caf16908c7348c7b51dad4c18 100644
--- a/cxxfem/src/femProblem.cpp
+++ b/cxxfem/src/femProblem.cpp
@@ -30,8 +30,9 @@ Problem::~Problem()
 
 void Problem::sync()
 {
+    std::cout << "Problem::sync()\n";
     // -- entities
-    std::cout << "retrieving entities...\n";
+    std::cout << "\tretrieving entities...\n";
     gmsh::vectorpair dimTags;
     gmsh::model::getEntities(dimTags);
     for (auto &p : dimTags)
@@ -41,7 +42,7 @@ void Problem::sync()
         entities_by_dimtag[std::make_pair(e->dim, e->tag)] = e;
 
     // -- physical groups
-    std::cout << "retrieving physical groups...\n";
+    std::cout << "\tretrieving physical groups...\n";
     dimTags.clear();
     gmsh::model::getPhysicalGroups(dimTags);
     for (auto &p : dimTags)
diff --git a/cxxfem/src/femSolver.cpp b/cxxfem/src/femSolver.cpp
index 384a1b384b666dcee0f1c5d36b06cb0bbbab847b..065bc22cc3a325ec45b349925c7cc3790c61e2fc 100644
--- a/cxxfem/src/femSolver.cpp
+++ b/cxxfem/src/femSolver.cpp
@@ -33,7 +33,7 @@ Solver::Solver(Problem &_pbl) : pbl(_pbl),
     bcs.add("FREE");
     bcs.add("PRESCRIBED");
 
-    std::cout << "#threads = " << OpenMP::get_max_threads() << '\n';
+    // std::cout << "#threads = " << OpenMP::get_max_threads() << '\n';
 }
 
 Solver::~Solver()
@@ -203,7 +203,7 @@ void Solver::fill_partition()
 
     // TODO: check that all dofs are assigned?
     timers["fill_partition"].stop();
-    std::cout << "\tcpu = " << timers["fill_partition"].read() << "s\n";
+    std::cout << "\telasped = " << timers["fill_partition"].read() << "s\n";
 }
 
 /// fill right-hand side of the system of equations
@@ -528,7 +528,7 @@ void Solver::fill_rhs()
             delete Kijv[i][j];
 
     timers["fill_rhs"].stop();
-    std::cout << "\tcpu = " << timers["fill_rhs"].read() << "s\n";
+    std::cout << "\telasped = " << timers["fill_rhs"].read() << "s\n";
 }
 
 /// compute right-hand side 'f' (sources)
@@ -628,7 +628,7 @@ void Solver::fill_lhs()
     }
 
     timers["fill_lhs"].stop();
-    std::cout << "\tcpu = " << timers["fill_lhs"].read() << "s\n";
+    std::cout << "\telasped = " << timers["fill_lhs"].read() << "s\n";
 
 }
 
@@ -667,7 +667,7 @@ void Solver::apply_bc()
     }
 
     timers["apply_bc"].stop();
-    std::cout << "\tcpu = " << timers["apply_bc"].read() << "s\n";
+    std::cout << "\telasped = " << timers["apply_bc"].read() << "s\n";
 }
 
 void Solver::solve_static()
@@ -771,7 +771,7 @@ void Solver::solve_static()
         delete tmp_rhs[i];
 
     timers["solve_static"].stop();
-    std::cout << "\tcpu = " << timers["solve_static"].read() << "s\n";
+    std::cout << "\telasped = " << timers["solve_static"].read() << "s\n";
 }
 
 // strain and stresses (epsilon, sigma) per node, per element
@@ -995,7 +995,7 @@ void Solver::compute_sig_epl(std::vector<size_t> &etags_all,
             } // for (auto &elems : entity->esets)
 
     timers["compute_sig_epl"].stop();
-    std::cout << "\tcpu = " << timers["compute_sig_epl"].read() << "s\n";
+    std::cout << "\telasped = " << timers["compute_sig_epl"].read() << "s\n";
 }
 
 FEM_API std::ostream &
diff --git a/models/boneload.py b/models/boneload.py
index e69afec118487e49eb57721992925a499564b6ae..103b77324d5cd2cfe839d4c6c89d5cfcf87f733a 100644
--- a/models/boneload.py
+++ b/models/boneload.py
@@ -248,10 +248,10 @@ def identify_nodes_SaP(coords, all_no, all_coords, eps=1e-3):
     cpu_stop = time.perf_counter()
 
     # print stats
-    print(f'sweep and prune done in {cpu_stop-cpu_start:.2f} seconds')
-    print(f'({ntests} tests instead of {len(coords)*len(all_coords)})')
-    print(f'{len(coords)} to be identified')
-    print(f'{len(ntags)} successfully identified')
+    print(f'\tsweep and prune done in {cpu_stop-cpu_start:.2f} seconds')
+    print(f'\t({ntests} tests instead of {len(coords)*len(all_coords)})')
+    print(f'\t{len(coords)} to be identified')
+    print(f'\t{len(ntags)} successfully identified')
     if(len(coords) != len(ntags)):
         print('WARNING: some nodes have not been identified!')
         raise Exception(f'ERROR: in identify_nodes_SaP!')
diff --git a/models/bonemodel.py b/models/bonemodel.py
index b266fd9bc8855761179e8f844f6b0cfee573cb74..d2624db217d9e0ab9ce0345e78de724ac934126c 100644
--- a/models/bonemodel.py
+++ b/models/bonemodel.py
@@ -335,8 +335,8 @@ def create_group(mesh_file_or_coord, all_no, all_coords, domain, groups, gname=N
         else:
             # one single point
             nodes = [ mesh_file_or_coord ]
-        print(f'nodes={nodes}')
-        print(f'mesh_file_or_coord={mesh_file_or_coord}')
+        print(f'\tnodes={nodes}')
+        print(f'\tmesh_file_or_coord={mesh_file_or_coord}')
         tris = []
     else:
         raise Exception(f'Data is not understood: {mesh_file_or_coord}')
@@ -356,7 +356,7 @@ def create_group(mesh_file_or_coord, all_no, all_coords, domain, groups, gname=N
     for tag in ntags:
         group.addMeshPoint(nodeset(tag))
 
-    print(f'group #{group.getNo()} ({gname}) ' \
+    print(f'\tgroup #{group.getNo()} ({gname}) ' \
     f'created from {mesh_file_or_coord} ({group.getNumberOfMeshPoints()} nodes)')
     
     return gname, nodes, tris, ntags
diff --git a/models/bonemodel2.py b/models/bonemodel2.py
index 62b2fcad8934297c984ee52ad08f8e4f940e1310..9cf1d86f6215bb3c69a2ea2cd802faba165cb8f8 100644
--- a/models/bonemodel2.py
+++ b/models/bonemodel2.py
@@ -76,7 +76,7 @@ def parms(d={}):
 def solve(p={}):
 
     import time
-    start_time = time.time()
+    start_time = time.perf_counter()
 
     p = parms(p)
 
@@ -139,20 +139,21 @@ def solve(p={}):
 
 
     # --------------------------------------------------------------------------
-    print('== ADDING MEDIUM...')
+    # print('== ADDING MEDIUM...')
+    print('')
 
     do_not_delete = []
 
     # material
     do_not_delete.append( fem.Medium(pbl, "all", E=p['Young'], nu=p['Poisson']) )
 
-    print('== ADDING DIRICHLETs...')
+    # print('== ADDING DIRICHLETs...')
     # axis of rotation
     for gname in ['axis_pt1', 'axis_pt2']:
         for d in p['fixations'][gname]:
             do_not_delete.append(fem.Dirichlet(pbl, gname, dof=d.upper(), val=0.0) )
 
-    print('== ADDING NEUMANNs...')
+    # print('== ADDING NEUMANNs...')
     # mshpts = geometry.getMesh().getPointSet()
     for name, gr in mgroups.items():
         # print(f'len(gr.ntags)={len(gr.ntags)}')
@@ -177,7 +178,7 @@ def solve(p={}):
     
     # Time integration scheme    
     # print (pbl)
-    print('SOLVE.......')
+    # print('SOLVE.......')
     pbl.quadrature = 'Gauss1' # 1 GP/tetra is enough
 
     solver = fem.Solver(pbl)
@@ -209,7 +210,7 @@ def solve(p={}):
         print(f'\t{name}_Fy = {v[1]:.3f} N')
         print(f'\t{name}_Fz = {v[2]:.3f} N')
 
-    end_time = time.time()
+    end_time = time.perf_counter()
     elapsed_time = end_time - start_time
     print(f'Total Execution time: {elapsed_time:.1f}s')
 
@@ -285,8 +286,8 @@ def create_group(mesh_file_or_coord, all_no, all_coords, groups, gname=None):
         else:
             # one single point
             nodes = [mesh_file_or_coord]
-        print(f'nodes={nodes}')
-        print(f'mesh_file_or_coord={mesh_file_or_coord}')
+        print(f'\tnodes={nodes}')
+        print(f'\tmesh_file_or_coord={mesh_file_or_coord}')
         tris = []
     else:
         raise Exception(f'Data is not understood: {mesh_file_or_coord}')
@@ -300,7 +301,7 @@ def create_group(mesh_file_or_coord, all_no, all_coords, groups, gname=None):
 
     # create a Gmsh group
     #   => a discrete entity with "points" elements
-    print(f'** adding gmsh group "{gname}"...')
+    print(f'adding gmsh group "{gname}"...')
     etag = gmsh.model.addDiscreteEntity(dim=0)
     point_type = gmsh.model.mesh.getElementType("Point", 1)
     # print(f'point_type={point_type}')
@@ -316,7 +317,7 @@ def create_group(mesh_file_or_coord, all_no, all_coords, groups, gname=None):
     # print(f'nods={nods}')
     # print(f'pos={pos}')
 
-    print(f'group {groups[gname]} ({gname}) '
+    print(f'\tgroup {groups[gname]} ({gname}) '
           f'created from {mesh_file_or_coord} ({len(nods)} nodes)')
 
     return gname, nodes, tris, ntags