Fix NaN error

srosenbu · Oct 16, 2023 · 97c54e6 · 97c54e6
1 parent 81abf18
commit 97c54e6
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 9 deletions.
diff --git a/python/comfe/helpers.py b/python/comfe/helpers.py
@@ -229,7 +229,7 @@ def diagonal_mass(function_space, rho, invert=True) -> df.fem.Function:
 
 
 def critical_timestep_1d(l_e, E, rho, order=1):
-    h_mesh = df.mesh.create_interval(MPI.COMM_SELF, np.array([0.0, l_e]), [1], cell_type=df.mesh.CellType.interval)
+    h_mesh = df.mesh.create_interval(MPI.COMM_SELF, 1, np.array([0.0, l_e]))
     h_P1 = df.fem.FunctionSpace(h_mesh, ("CG", order))
     h_u, h_v = ufl.TrialFunction(h_P1), ufl.TestFunction(h_P1)
     K_form = df.fem.form(E * ufl.inner(ufl.grad(h_u), ufl.grad(h_v)) * ufl.dx)

diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs
@@ -124,10 +124,15 @@ impl ConstitutiveModel for GradientJH23D {
                 if p_trial > p_damaged {
                     p_trial
                 } else {
-                    let denominator = (density_1-density_0)/self.parameters.RHO;
-                    let K_pl = (p_damaged - p_0) / denominator;
-                    let K_el = (p_trial - p_0) / denominator;
-                    d_eps_vol_pl = (1. - K_pl / K_el) * d_eps_vol;
+                    let pl = p_damaged - p_0;
+                    let el = p_trial - p_0;
+                    d_eps_vol_pl = {
+                        if el != 0.0 {
+                            (1. - pl / el) * d_eps_vol
+                        } else {
+                            d_eps_vol
+                        }
+                    };
                     p_damaged
                 }
             }

diff --git a/src/jh2.rs b/src/jh2.rs
@@ -143,10 +143,16 @@ impl ConstitutiveModel for JH23D {
                 if p_trial > p_damaged {
                     p_trial
                 } else {
-                    let denominator = (density_1-density_0)/self.parameters.RHO;
-                    let K_pl = (p_damaged - p_0) / denominator;
-                    let K_el = (p_trial - p_0) / denominator;
-                    d_eps_vol_pl = (1. - K_pl / K_el) * d_eps_vol;
+                    //let denominator = (density_1-density_0)/self.parameters.RHO;
+                    let pl = p_damaged - p_0;
+                    let el = p_trial - p_0;
+                    d_eps_vol_pl = {
+                        if el != 0.0 {
+                            (1. - pl / el) * d_eps_vol
+                        } else {
+                            d_eps_vol
+                        }
+                    };
                     p_damaged
                 }
             }