Energy calculations now seem to be consistent

python/comfe/cdm.py

@@ -64,6 +64,7 @@ class CDM3D(CDMSolver):
     # fields: dict[str, df.fem.Function]
     # q_fields: dict[str, df.fem.Function]
     total_energy: float | None = None
+    distortion_expr: df.fem.Expression | None = None
     # energy_form: df.fem.FormMetaClass | None = None
 
     class Config:
@@ -80,11 +81,22 @@ def __init__(
         quadrature_rule: QuadratureRule,
         additional_output: list[str] | None = None,
         calculate_total_energy: bool = False,
+        monitor_distortion: bool = False,
     ) -> None:
         # self.del_t = None
         v = df.fem.Function(function_space, name="Velocity")
         u = df.fem.Function(function_space, name="Displacements")
         f = df.fem.Function(function_space, name="Forces")
+        fields = {"u": u, "v": v, "f": f}
+        if monitor_distortion:
+            detJ_space = df.fem.FunctionSpace(function_space.mesh, ("DG", 0))
+            detJ = df.fem.Function(detJ_space, name="detJ")
+            distortion_expr = df.fem.Expression(
+                ufl.JacobianDeterminant(detJ_space.mesh), detJ_space.element.interpolation_points()
+            )
+            fields["detJ"] = detJ
+        else:
+            distortion_expr = None
 
         model = ConstitutiveModel(
             rust_model,
@@ -131,9 +143,10 @@ def __init__(
             model=model,
             M=M,
             # nonlocal_var=nonlocal_var,
-            fields={"u": u, "v": v, "f": f},
+            fields=fields,
             q_fields=model.output,
             total_energy=total_energy,
+            distortion_expr=distortion_expr,
         )
 
     def _as_3d_tensor(self, T: ufl.core.expr.Expr):

python/comfe/helpers.py

@@ -224,6 +224,31 @@ def diagonal_mass(function_space, rho, invert=True) -> df.fem.Function:
     return M_action
 
 
+# class CriticalTimestep(mesh, rho, K, G, order=1):
+#    pass
+
+
+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_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)
+    M_form = df.fem.form(rho * ufl.inner(h_u, h_v) * ufl.dx)
+
+    h_K, h_M = (
+        df.fem.petsc.assemble_matrix(K_form),
+        df.fem.petsc.assemble_matrix(M_form),
+    )
+    h_K.assemble()
+    h_M.assemble()
+    h_M = np.array(h_M[:, :])
+    h_K = np.array(h_K[:, :])
+    max_eig = np.linalg.norm(eigvals(h_K, h_M), np.inf)
+
+    h = 2.0 / max_eig**0.5
+    return h
+
+
 def critical_timestep(l_x, l_y, G, K, rho, cell_type=df.mesh.CellType.quadrilateral, order=1):
     # todo: implement other cell_types
     # cell_type=mesh.topology.cell_type

src/gradient_jh2.rs

@@ -124,7 +124,10 @@ impl ConstitutiveModel for GradientJH23D {
                 if p_trial > p_damaged {
                     p_trial
                 } else {
-                    d_eps_vol_pl = d_eps_vol - (p_damaged - p_0) / (self.parameters.K1 * del_t);
+                    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;
                     p_damaged
                 }
             }

src/jh2.rs

@@ -124,6 +124,7 @@ impl ConstitutiveModel for JH23D {
         let f1 = del_t / 2. * 3. * d_eps_vol;
         let density_0 = input.get_scalar(Q::Density, ip);
         let density_1 = density_0 * (1. - f1) / (1. + f1);
+        assert!(density_1 > 0.0, "Negative density encountered in JH2 model: {}", density_1);
         output.set_scalar(Q::Density, ip, density_1);
 
         let mu = density_1 / self.parameters.RHO - 1.;
@@ -142,10 +143,12 @@ impl ConstitutiveModel for JH23D {
                 if p_trial > p_damaged {
                     p_trial
                 } else {
-                    d_eps_vol_pl = d_eps_vol - (p_damaged - p_0) / (self.parameters.K1 * del_t);
+                    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;
                     p_damaged
                 }
-                //(self.parameters.K1 * mu).max(-self.parameters.T * (1. - damage_1))
             }
         };
         if damage_1 > damage_0 {