Tidy up docs for strain optimization. Renamed StrainOpt -> CoilStrain.

examples/2_Intermediate/strain_optimization.py

@@ -13,7 +13,7 @@
 
 import numpy as np
 from scipy.optimize import minimize
-from simsopt.geo import StrainOpt, LPTorsionalStrainPenalty, LPBinormalCurvatureStrainPenalty
+from simsopt.geo import CoilStrain, LPTorsionalStrainPenalty, LPBinormalCurvatureStrainPenalty
 from simsopt.geo import FrameRotation, FramedCurveFrenet, CurveXYZFourier
 from simsopt.configs import get_hsx_data
 from simsopt.util import in_github_actions
@@ -40,7 +40,7 @@
 Jbin = LPBinormalCurvatureStrainPenalty(
     framedcurve, p=2, threshold=cur_threshold)
 
-strain = StrainOpt(framedcurve, width)
+strain = CoilStrain(framedcurve, width)
 JF = Jtor + Jbin
 
 

src/simsopt/geo/__init__.py

@@ -22,7 +22,6 @@
 from .surfacexyzfourier import *
 from .surfacexyztensorfourier import *
 from .strain_optimization import *
-from .framedcurve import *
 
 from .permanent_magnet_grid import *
 

src/simsopt/geo/finitebuild.py

@@ -1,12 +1,5 @@
 import numpy as np
-import jax.numpy as jnp
-from jax import vjp, jvp
-
-import simsoptpp as sopp
-from .._core.optimizable import Optimizable
-from .._core.derivative import Derivative
-from .curve import Curve
-from .jit import jit
+
 from .framedcurve import FramedCurve, FrameRotation, ZeroRotation, FramedCurveCentroid, FramedCurveFrenet
 
 """
@@ -50,7 +43,6 @@ def recompute_bell(self, parent=None):
     def gamma_impl(self, gamma, quadpoints):
         assert quadpoints.shape[0] == self.curve.quadpoints.shape[0]
         assert np.linalg.norm(quadpoints - self.curve.quadpoints) < 1e-15
-        c = self.curve
         t, n, b = self.framedcurve.rotated_frame()
         gamma[:] = self.curve.gamma() + self.dn * n + self.db * b
 

src/simsopt/geo/framedcurve.py

@@ -261,7 +261,6 @@ def frame_torsion(self):
         gamma = self.curve.gamma()
         d1gamma = self.curve.gammadash()
         d2gamma = self.curve.gammadashdash()
-        d3gamma = self.curve.gammadashdashdash()
         alpha = self.rotation.alpha(self.curve.quadpoints)
         alphadash = self.rotation.alphadash(self.curve.quadpoints)
         return self.torsion(gamma, d1gamma, d2gamma, alpha, alphadash)
@@ -270,7 +269,6 @@ def frame_binormal_curvature(self):
         gamma = self.curve.gamma()
         d1gamma = self.curve.gammadash()
         d2gamma = self.curve.gammadashdash()
-        d3gamma = self.curve.gammadashdashdash()
         alpha = self.rotation.alpha(self.curve.quadpoints)
         alphadash = self.rotation.alphadash(self.curve.quadpoints)
         return self.binorm(gamma, d1gamma, d2gamma, alpha, alphadash)
@@ -422,7 +420,7 @@ def __init__(self, quadpoints):
 
         .. code-block:: python
 
-            rot = FilamentRotation(...)
+            rot = FrameRotation(...)
             rot.fix_all()
 
         """
@@ -604,12 +602,6 @@ def inner(a, b):
     return np.sum(a*b, axis=1)
 
 
-torsion2vjp0 = jit(lambda ndash, b, v: vjp(
-    lambda nd: torsion_pure(nd, b), ndash)[1](v)[0])
-torsion2vjp1 = jit(lambda ndash, b, v: vjp(
-    lambda bi: torsion_pure(ndash, bi), b)[1](v)[0])
-
-
 def torsion_pure_frenet(gamma, gammadash, gammadashdash, gammadashdashdash,
                         alpha, alphadash):
     """Torsion function for export/evaulate coil sets"""

src/simsopt/geo/strain_optimization.py

@@ -6,7 +6,7 @@
 from simsopt.geo.curveobjectives import Lp_torsion_pure
 
 __all__ = ['LPBinormalCurvatureStrainPenalty',
-           'LPTorsionalStrainPenalty', 'StrainOpt']
+           'LPTorsionalStrainPenalty', 'CoilStrain']
 
 
 class LPBinormalCurvatureStrainPenalty(Optimizable):
@@ -15,14 +15,21 @@ class LPBinormalCurvatureStrainPenalty(Optimizable):
     of the binormal curvature strain, and penalizes where the local strain exceeds a threshold
 
     .. math::
-        J = \frac{1}{p} \int_{\text{curve}} \text{max}(\epsilon_{\text{bend}} - \epsilon_0, 0)^p ~dl
+        J = \frac{1}{p} \int_{\text{curve}} \text{max}(\epsilon_{\text{bend}} - \epsilon_0, 0)^p ~dl,
 
-    where :math:`\epsilon_0` is a threshold strain, given by the argument ``threshold``.
+    where
+
+    .. math::
+        \epsilon_{\text{bend}} = \frac{w |\hat{\textbf{b}} \cdot \boldsymbol{\kappa}|}{2},
+
+    :math:`w` is the width of the tape, :math:`\hat{\textbf{b}}` is the 
+    frame binormal vector, :math:`\boldsymbol{\kappa}` is the curvature vector of the 
+    filamentary coil, and :math:`\epsilon_0` is a threshold strain, given by the argument ``threshold``.
     """
 
     def __init__(self, framedcurve, width=1e-3, p=2, threshold=0):
         self.framedcurve = framedcurve
-        self.strain = StrainOpt(framedcurve, width)
+        self.strain = CoilStrain(framedcurve, width)
         self.width = width
         self.p = p
         self.threshold = threshold
@@ -65,12 +72,18 @@ class LPTorsionalStrainPenalty(Optimizable):
     .. math::
         J = \frac{1}{p} \int_{\text{curve}} \text{max}(\epsilon_{\text{tor}} - \epsilon_0, 0)^p ~dl
 
-    where :math:`\epsilon_0` is a threshold strain, given by the argument ``threshold``.
+    where
+
+    .. math::
+        \epsilon_{\text{tor}} = \frac{\tau^2 w^2}{12},
+
+    :math:`\tau` is the torsion of the tape frame, :math:`w` is the width of the tape,
+    and :math:`\epsilon_0` is a threshold strain, given by the argument ``threshold``.
     """
 
     def __init__(self, framedcurve, width=1e-3, p=2, threshold=0):
         self.framedcurve = framedcurve
-        self.strain = StrainOpt(framedcurve, width)
+        self.strain = CoilStrain(framedcurve, width)
         self.width = width
         self.p = p
         self.threshold = threshold
@@ -105,7 +118,7 @@ def dJ(self):
     return_fn_map = {'J': J, 'dJ': dJ}
 
 
-class StrainOpt(Optimizable):
+class CoilStrain(Optimizable):
     r"""
     This class evaluates the torsional and binormal curvature strains on HTS, based on
     a filamentary model of the coil and the orientation of the HTS tape. 
@@ -119,13 +132,18 @@ class StrainOpt(Optimizable):
     the expressions for the strains are: 
 
     .. math::
-        \epsilon_{\text{tor}} = \frac{\tau^2 w^2}{12} 
+        \epsilon_{\text{tor}} = \frac{\tau^2 w^2}{12}
+
         \epsilon_{\text{bend}} = \frac{w |\hat{\textbf{b}} \cdot \boldsymbol{\kappa}|}{2},
 
     where :math:`\tau` is the torsion of the tape frame, :math:`\hat{\textbf{b}}` is the 
-    frame binormal vector, and :math:`\boldsymbol{\kappa}` is the curvature vector of the 
-    filamentary coil. 
+    frame binormal vector, :math:`\boldsymbol{\kappa}` is the curvature vector of the 
+    filamentary coil, and :math:`w` is the width of the tape.
 
+    This class is not intended to be used as an objective function inside
+    optimization. For that purpose you should instead use
+    :obj:`LPBinormalCurvatureStrainPenalty` or :obj:`LPTorsionalStrainPenalty`.
+    Those classes also compute gradients whereas this class does not.
     """
 
     def __init__(self, framedcurve, width=1e-3):
@@ -143,14 +161,14 @@ def __init__(self, framedcurve, width=1e-3):
         super().__init__(depends_on=[framedcurve])
 
     def torsional_strain(self):
-        """
+        r"""
         Returns the value of the torsional strain, :math:`\epsilon_{\text{tor}}`, along 
         the quadpoints defining the filamentary coil. 
         """
         return self.torstrain_jax(self.framedcurve.frame_torsion(), self.width)
 
     def binormal_curvature_strain(self):
-        """
+        r"""
         Returns the value of the torsional strain, :math:`\epsilon_{\text{bend}}`, along 
         the quadpoints defining the filamentary coil. 
         """

tests/geo/test_strainopt.py

@@ -7,7 +7,7 @@
 from scipy.optimize import minimize
 
 
-class StrainOptTesting(unittest.TestCase):
+class CoilStrainTesting(unittest.TestCase):
 
     def test_strain_opt(self):
         """