From 170ae6efab89958aa53b89471ce993d30877cf21 Mon Sep 17 00:00:00 2001
From: Romain Boman <r.boman@uliege.be>
Date: Sun, 29 May 2022 14:06:31 +0200
Subject: [PATCH] use 1GP/tetra

---
 cxxfem/src/femProblem.h |  3 ++-
 models/bonemodel2.py    | 25 +++++++++++++------------
 2 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/cxxfem/src/femProblem.h b/cxxfem/src/femProblem.h
index fec3242..fde0369 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 6f8e08f..85c5d59 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"
-- 
GitLab