signed curvature

desc/compute/_curve.py

@@ -1,6 +1,6 @@
 from interpax import interp1d
 
-from desc.backend import jnp
+from desc.backend import jnp, sign
 
 from .data_index import register_compute_fun
 from .geom_utils import rotation_matrix, rpz2xyz, rpz2xyz_vec, xyz2rpz, xyz2rpz_vec
@@ -146,6 +146,50 @@ def _Z_Curve(params, transforms, profiles, data, **kwargs):
     return data
 
 
+@register_compute_fun(
+    name="center",
+    label="\\langle\\mathbf{x}\\rangle",
+    units="m",
+    units_long="meters",
+    description="Centroid of the curve",
+    dim=3,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="s",
+    data=["x"],
+    parameterization="desc.geometry.core.Curve",
+)
+def _center_Curve(params, transforms, profiles, data, **kwargs):
+    data["center"] = jnp.mean(data["x"], axis=0) * jnp.ones_like(data["x"])
+    return data
+
+
+@register_compute_fun(
+    name="center",
+    label="\\langle\\mathbf{x}\\rangle",
+    units="m",
+    units_long="meters",
+    description="Centroid of the curve",
+    dim=3,
+    params=["center", "shift"],
+    transforms={},
+    profiles=[],
+    coordinates="s",
+    data=["x"],
+    parameterization="desc.geometry.curve.FourierPlanarCurve",
+    basis_in="{'rpz', 'xyz'}: Basis for input params vectors, Default 'xyz'",
+)
+def _center_FourierPlanarCurve(params, transforms, profiles, data, **kwargs):
+    # convert to rpz
+    if kwargs.get("basis_in", "xyz").lower() == "rpz":
+        center = params["center"]
+    else:
+        center = xyz2rpz(params["center"])
+    data["center"] = (center + xyz2rpz(params["shift"])) * jnp.ones_like(data["x"])
+    return data
+
+
 @register_compute_fun(
     name="x",
     label="\\mathbf{x}",
@@ -790,13 +834,12 @@ def _frenet_tangent(params, transforms, profiles, data, **kwargs):
     transforms={},
     profiles=[],
     coordinates="s",
-    data=["x_ss"],
+    data=["x_s", "x_ss", "frenet_tangent"],
     parameterization="desc.geometry.core.Curve",
 )
 def _frenet_normal(params, transforms, profiles, data, **kwargs):
-    data["frenet_normal"] = (
-        data["x_ss"] / jnp.linalg.norm(data["x_ss"], axis=-1)[:, None]
-    )
+    normal = cross(data["x_s"], cross(data["x_ss"], data["x_s"]))
+    data["frenet_normal"] = normal / jnp.linalg.norm(normal, axis=-1)[:, None]
     return data
 
 
@@ -826,20 +869,22 @@ def _frenet_binormal(params, transforms, profiles, data, **kwargs):
     label="\\kappa",
     units="m^{-1}",
     units_long="Inverse meters",
-    description="Scalar curvature of the curve",
+    description="Scalar curvature of the curve, positive/negative = concave/convex",
     dim=1,
     params=[],
     transforms={},
     profiles=[],
     coordinates="s",
-    data=["x_s", "x_ss"],
+    data=["center", "x", "x_s", "x_ss", "frenet_normal"],
     parameterization="desc.geometry.core.Curve",
 )
 def _curvature(params, transforms, profiles, data, **kwargs):
+    # magnitude of curvature
     dxn = jnp.linalg.norm(data["x_s"], axis=-1)[:, jnp.newaxis]
-    data["curvature"] = jnp.linalg.norm(
-        cross(data["x_s"], data["x_ss"]) / dxn**3, axis=-1
-    )
+    curvature = jnp.linalg.norm(cross(data["x_s"], data["x_ss"]) / dxn**3, axis=-1)
+    # sign of curvature (positive = "concave", negative = "convex")
+    r = data["x"] - data["center"]
+    data["curvature"] = curvature * sign(dot(r, data["frenet_normal"]))
     return data
 
 

desc/objectives/_coils.py

@@ -332,8 +332,9 @@ def compute(self, params, constants=None):
 class CoilCurvature(_CoilObjective):
     """Coil curvature.
 
-    Targets the local curvature value per grid node for each coil. A smaller curvature
-    value indicates straighter coils. All curvature values are positive.
+    Targets the local curvature at each grid node for each coil. Positive curvature
+    corresponds to "concave" curves, while negative curvature corresponds to "convex"
+    curves. Curvature values closer to 0 indicate straighter coils.
 
     Parameters
     ----------
@@ -392,7 +393,7 @@ def __init__(
         name="coil curvature",
     ):
         if target is None and bounds is None:
-            bounds = (0, 1)
+            bounds = (-1, 0)
 
         super().__init__(
             coil,
@@ -450,9 +451,8 @@ def compute(self, params, constants=None):
 class CoilTorsion(_CoilObjective):
     """Coil torsion.
 
-    Targets the local torsion value per grid node for each coil. Indicative
-    of how much the coil goes out of the poloidal plane. e.g. a torsion
-    value of 0 means the coil is completely planar.
+    Targets the local torsion value at each grid node for each coil. Indicative of how
+    non-planar the coil is (a torsion value of 0 means the coil is perfectly planar).
 
     Parameters
     ----------

tests/test_objective_funs.py

@@ -825,7 +825,7 @@ def test(coil, grid=None):
             obj = CoilCurvature(coil, grid=grid)
             obj.build()
             f = obj.compute(params=coil.params_dict)
-            np.testing.assert_allclose(f, 1 / 2, rtol=1e-8)
+            np.testing.assert_allclose(f, -1 / 2, rtol=1e-8)
             assert len(f) == obj.dim_f
 
         coil = FourierPlanarCoil(r_n=2)