From cb027ff43133596ea2d5a44e5c561a93ecb7760d Mon Sep 17 00:00:00 2001
From: Jed Brown <jed@jedbrown.org>
Date: Thu, 11 Apr 2024 15:40:11 -0600
Subject: [PATCH] samples: add volumetric: hyperelasticity example with trait
 (#6)

* samples: add volumetric: hyperelasticity example with trait

* Update samples/tests/traits/mod.rs

Co-authored-by: Jed Brown <jed@jedbrown.org>

* add small tolerance to test

---------

Co-authored-by: Manuel Drehwald <git@manuel.drehwald.info>
---
 samples/tests/mod.rs        |  5 +++--
 samples/tests/traits/mod.rs | 43 +++++++++++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+), 2 deletions(-)
 create mode 100644 samples/tests/traits/mod.rs

diff --git a/samples/tests/mod.rs b/samples/tests/mod.rs
index 457923a..3f1317a 100644
--- a/samples/tests/mod.rs
+++ b/samples/tests/mod.rs
@@ -1,7 +1,8 @@
 #![feature(autodiff)]
 
 mod forward;
-mod neohookean;
 mod reverse;
-mod higher;
 
+mod higher;
+mod neohookean;
+mod traits;
diff --git a/samples/tests/traits/mod.rs b/samples/tests/traits/mod.rs
new file mode 100644
index 0000000..65e8f5e
--- /dev/null
+++ b/samples/tests/traits/mod.rs
@@ -0,0 +1,43 @@
+samples::test! {
+    volumetric;
+    /// ANCHOR: volumetric
+    trait Volumetric {
+        /// Strain energy density
+        fn psi(&self, j: f64) -> f64;
+        /// Derivative of strain energy with respect to $J$
+        fn d_psi(&self, j: f64, b_psi: f64) -> (f64, f64);
+        /// The volumetric contribution to the second Piola-Kirchhoff stress is
+        /// a scalar multiplied by $C^{-1}$ where $C = I + 2E$ in terms of the
+        /// Green-Lagrange strain $E$. The derivative of $J$ with respect to $E$
+        /// is $J C^{-1}$. We'll call the volumetric contribution that scalar
+        /// multiple of $C^{-1}$.
+        fn stress(&self, j: f64) -> f64 {
+            let (_, d_psi) = self.d_psi(j, 1.0);
+            d_psi * j
+        }
+    }
+
+    struct Ogden {
+        k: f64,
+    }
+    impl Ogden {
+        pub fn stress_analytic(&self, j: f64) -> f64 {
+            self.k * 0.5 * (j * j - 1.0)
+        }
+    }
+    impl Volumetric for Ogden {
+        #[autodiff(d_psi, Reverse, Const, Active, Active)]
+        fn psi(&self, j: f64) -> f64 {
+            self.k * 0.25 * (j * j - 1.0 - 2.0 * j.ln())
+        }
+    }
+
+    fn main() {
+        let j = 0.8;
+        let vol = Ogden { k: 1.0 };
+        let s = vol.stress(j);
+        let s_ref = vol.stress_analytic(j);
+        assert!((s - s_ref).abs() < 1e-15, "{}", s - s_ref);
+    }
+    // ANCHOR_END: volumetric
+}