diff --git a/cxxfem/src/femProblem.h b/cxxfem/src/femProblem.h
index fec3242b5ff48f485a8d3340cb5b25c7343edf8f..fde0369c08cfe720fb23f3c6af1da69935426d13 100644
--- a/cxxfem/src/femProblem.h
+++ b/cxxfem/src/femProblem.h
@@ -28,9 +28,10 @@ public:
     std::vector<Medium *> media;
     std::vector<Dirichlet *> dirichlets;
     std::vector<Neumann *> neumanns;
+#endif
 
     std::string quadrature;
-#endif
+
 public:
     Problem();
     ~Problem();
diff --git a/models/bonemodel2.py b/models/bonemodel2.py
index 6f8e08fd56c0a8af11b4c43c84e3fdcf8998ac3c..85c5d596480f749827fdb9587ecec6b3dd014343 100644
--- a/models/bonemodel2.py
+++ b/models/bonemodel2.py
@@ -48,26 +48,26 @@ def parms(d={}):
         'axis_pt2': ['x', 'y', 'z']
     }
     # rigid cylinder (used if p['food']=='cylinder')
-    R = 20
-    p['cylinder'] = {
-        'origin': [-R-8, 0, -150],
-        'radius': R,
-        'width': 100,
-        'friction': 0.1,
-        'penalty': 1e3
-    }
+    # R = 20
+    # p['cylinder'] = {
+    #     'origin': [-R-8, 0, -150],
+    #     'radius': R,
+    #     'width': 100,
+    #     'friction': 0.1,
+    #     'penalty': 1e3
+    # }
 
     # material properties
-    p['density'] = 1.850e-9  # [T/mm³] - bone: 1.850 kg/l
+    # p['density'] = 1.850e-9  # [T/mm³] - bone: 1.850 kg/l
     p['Young'] = 17000.      # [MPa] elastic modulus - bone: 17-20 GPa
     p['Poisson'] = 0.3       # [-] Poisson's ratio
 
     # numerical parameters
-    p['tolNR'] = 1e-6        # [-] equilibrium tolerance
-    p['dt0'] = 1.0           # [s] time step size
+    # p['tolNR'] = 1e-6        # [-] equilibrium tolerance
+    # p['dt0'] = 1.0           # [s] time step size
 
     # gmsh toolbox
-    p['use_gmshOld'] = False  # use old gmsh interface
+    # p['use_gmshOld'] = False  # use old gmsh interface
 
     p.update(d)
     return p
@@ -175,6 +175,7 @@ def solve(p={}):
     # Time integration scheme    
     # print (pbl)
     print('SOLVE.......')
+    pbl.quadrature = 'Gauss1' # 1 GP/tetra is enough
 
     solver = fem.Solver(pbl)
     solver.linear_solver = "pardiso"