use 1GP/tetra

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();

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"