From 31c76a6ae98ea8ce36a09aeb044a3f304251a52e Mon Sep 17 00:00:00 2001
From: Mark Messner <messner@anl.gov>
Date: Fri, 24 Jan 2025 15:02:38 -0600
Subject: [PATCH] Fix the Bree cylinder example

---
 examples/bree.py |  4 ++--
 neml/axisym.py   | 15 ++++++++++++---
 2 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/examples/bree.py b/examples/bree.py
index a0568252..d0a189d6 100755
--- a/examples/bree.py
+++ b/examples/bree.py
@@ -35,7 +35,7 @@ def Pfn(t):
       nhalf, 0, ncycles, delay = delay, ndelay = nhalf)
   dts = np.diff(times)
 
-  problem = axisym.BreeProblem(rs, mats, ns, Tfn, Pfn)
+  problem = axisym.BreeProblem(rs, mats, ns, Tfn, Pfn, specify_hoop_stress = True)
 
   for t in times[1:]:
     problem.step(t)
@@ -159,7 +159,7 @@ def grid_classical(E, nu, sY, alpha,
   nu = 0.25
   sY = 100.0
   alpha = 1.0e-5
-  n = 20
+  n = 12
   
   points, conditions, dx, dy, nx, ny = grid_classical(
       E, nu, sY, alpha, max_x = 0.9999, max_y = 5.0, nx = 15, ny = 75,
diff --git a/neml/axisym.py b/neml/axisym.py
index 21592f64..01a779fc 100644
--- a/neml/axisym.py
+++ b/neml/axisym.py
@@ -283,7 +283,7 @@ class BreeProblem(VesselSectionProblem):
     pressure.
   """
   def __init__(self, rs, mats, ns, T, p, rtol = 1.0e-6, atol = 1.0e-8,
-      itype = "gauss", p_ext = lambda t: 0.0):
+      itype = "gauss", p_ext = lambda t: 0.0, specify_hoop_stress = False):
     """
       Parameters:
         rs      radii deliminating each region
@@ -298,9 +298,12 @@ def __init__(self, rs, mats, ns, T, p, rtol = 1.0e-6, atol = 1.0e-8,
         bias    bias the mesh towards the interfaces
         itype   "rectangle" or "gauss", defaults to "gauss"
         p_ext   external pressure as a function of t
+        specify_hoop_stress if True, directly specify a hoop stress rather than calculate it
     """
     super(BreeProblem, self).__init__(rs, mats, ns, T, p,
         rtol = rtol, atol = atol, p_ext = p_ext)
+
+    self.specify_hoop_stress = specify_hoop_stress
     
     self.regions = np.diff(self.rs)
     self.dlength = [r for r,n in zip(self.regions, self.ns) for i in range(n)]
@@ -320,7 +323,10 @@ def __init__(self, rs, mats, ns, T, p, rtol = 1.0e-6, atol = 1.0e-8,
 
     self.npoints = len(self.ipoints)
     
-    self.P = lambda t: (p(t) * self.r_inner - p_ext(t) * self.r_outer) / self.t
+    if self.specify_hoop_stress:
+        self.P = p
+    else:
+        self.P = lambda t: (p(t) * self.r_inner - p_ext(t) * self.r_outer) / self.t
 
     # Setup the n bar model
     self.barmodel = arbbar.BarModel()
@@ -360,7 +366,10 @@ def update_loading(self, new_p, new_T, new_p_ext = lambda t: 0.0):
     """
       Impose new loading functions
     """
-    self.P = lambda t: (new_p(t) * self.r_inner - new_p_ext(t) * self.r_outer) / self.t
+    if self.specify_hoop_stress:
+        self.P = new_p
+    else:
+        self.P = lambda t: (new_p(t) * self.r_inner - new_p_ext(t) * self.r_outer) / self.t
     self.barmodel.nodes[2]['force bc'] = lambda t: self.P(t) * np.sum(self.weights)
 
     self.T = new_T