Regcoil branch

CMakeLists.txt

@@ -132,6 +132,10 @@ pybind11_add_module(${PROJECT_NAME}
     src/simsoptpp/magneticfield_wireframe.cpp
     src/simsoptpp/boozerradialinterpolant.cpp
     src/simsoptpp/wireframe_optimization.cpp
+    src/simsoptpp/python_currentpotential.cpp
+    src/simsoptpp/currentpotential.cpp 
+    src/simsoptpp/currentpotentialfourier.cpp
+    src/simsoptpp/winding_surface.cpp
     )
 
 set_target_properties(${PROJECT_NAME}

examples/3_Advanced/winding_surface_advanced.py

@@ -0,0 +1,137 @@
+#!/usr/bin/env python3
+"""
+Cut Coils: Extract discrete coils from a current potential on a winding surface.
+
+This example demonstrates how to extract coil contours from a regularized current
+potential (from REGCOIL or simsopt-regcoil) and convert them into simsopt Coil
+objects. The workflow includes:
+
+1. Loading current potential data from a NetCDF file (simsopt-regcoil or legacy REGCOIL)
+2. Computing the current potential φ and |∇φ| on a (θ, ζ) grid
+3. Selecting contours by:
+   - Interactive double-click selection
+   - Specifying (θ, ζ) points that lie on desired contours
+   - Specifying contour level values
+   - Choosing N contours per field period
+4. Classifying contours as window-pane (closed), modular, or helical
+5. Computing currents for each contour
+6. Converting contours to 3D curves and Coil objects
+7. Applying stellarator symmetry to get the full coil set
+
+Usage
+-----
+From the simsopt root directory:
+
+    python winding_surface_advanced.py --surface regcoil_out.hsx.nc
+
+With custom options:
+
+    python winding_surface_advanced.py \\
+        --surface /path/to/simsopt_regcoil_out.nc \\
+        --ilambda 2 \\
+        --points "0.5,1.0" "1.0,0.5" \\
+        --output my_output_dir
+
+For interactive mode (double-click to select contours):
+
+    python winding_surface_advanced.py --surface file.nc --interactive
+
+Requirements
+------------
+- A simsopt-regcoil or REGCOIL NetCDF output file
+- For legacy REGCOIL: a surface file (nescin format) for 3D mapping
+- matplotlib for contour visualization
+"""
+
+import argparse
+from pathlib import Path
+
+# Add parent to path for standalone run
+import sys
+sys.path.insert(0, str(Path(__file__).resolve().parents[2]))
+
+from simsopt.util import run_cut_coils
+
+# Default test file
+TEST_DIR = (Path(__file__).parent / ".." / ".." / "tests" / "test_files").resolve()
+DEFAULT_REGCOIL = TEST_DIR / "regcoil_out.hsx.nc"
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Extract coils from a current potential (REGCOIL / simsopt-regcoil output).",
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        epilog=__doc__,
+    )
+    parser.add_argument(
+        "--surface",
+        type=Path,
+        default=DEFAULT_REGCOIL,
+        help="REGCOIL or simsopt-regcoil NetCDF file (relative paths are resolved against tests/test_files).",
+    )
+    parser.add_argument(
+        "--ilambda",
+        type=int,
+        default=-1,
+        help="Lambda index (0-based) to use.",
+    )
+    parser.add_argument(
+        "--points",
+        nargs="*",
+        type=str,
+        default=None,
+        help='(θ,ζ) points as "theta,zeta" for contour selection.',
+    )
+    parser.add_argument(
+        "--no-sv",
+        action="store_false",
+        help="Use full (multi-valued) current potential with net toroidal/poloidal currents.",
+    )
+    parser.add_argument(
+        "--interactive",
+        action="store_true",
+        help="Use interactive double-click contour selection.",
+    )
+    parser.add_argument(
+        "--output",
+        type=Path,
+        default=None,
+        help="Output directory for VTK and JSON (default: winding_surface_<surface_stem>/).",
+    )
+    parser.add_argument(
+        "--no-plot",
+        action="store_true",
+        help="Do not show final coil plot.",
+    )
+    parser.add_argument(
+        "--save",
+        action="store_true",
+        help="Save coils to JSON file.",
+    )
+    args = parser.parse_args()
+
+    points = None
+    if args.points:
+        points = []
+        for s in args.points:
+            parts = s.split(",")
+            if len(parts) != 2:
+                parser.error(f"Invalid point format: {s}. Use theta,zeta")
+            points.append([float(parts[0]), float(parts[1])])
+
+    surface_filename = args.surface if args.surface.is_absolute() else TEST_DIR / args.surface
+    output_path = args.output or Path(f"winding_surface_{surface_filename.stem}")
+    run_cut_coils(
+        surface_filename=surface_filename,
+        ilambda=args.ilambda,
+        points=points,
+        single_valued=not args.no_sv,
+        interactive=args.interactive,
+        show_final_coilset=not args.no_plot,
+        write_coils_to_file=args.save,
+        output_path=output_path,
+    )
+
+
+if __name__ == "__main__":
+    main()

examples/run_serial_examples

@@ -20,6 +20,7 @@ set -ex
 ./2_Intermediate/stage_two_optimization_stochastic.py
 ./2_Intermediate/stage_two_optimization_finite_beta.py
 ./2_Intermediate/strain_optimization.py
+./2_Intermediate/winding_surface.py
 ./2_Intermediate/permanent_magnet_MUSE.py
 ./2_Intermediate/permanent_magnet_QA.py
 ./2_Intermediate/permanent_magnet_PM4Stell.py
@@ -30,3 +31,7 @@ set -ex
 ./3_Advanced/stage_two_optimization_finitebuild.py
 ./3_Advanced/coil_forces.py
 ./3_Advanced/wireframe_gsco_multistep.py
+./3_Advanced/winding_surface_advanced.py --surface regcoil_out.hsx.nc
+./3_Advanced/winding_surface_advanced.py --surface regcoil_out.w7x.nc --ilambda 0
+./3_Advanced/winding_surface_advanced.py --surface regcoil_out.near_axis.nc --points "0.5,0.5" "1.0,0.3"
+./3_Advanced/winding_surface_advanced.py --surface regcoil_out.near_axis_asym.nc --no-sv --output winding_surface_near_axis_asym

src/simsopt/field/__init__.py

@@ -7,6 +7,7 @@
 from .mgrid import *
 from .normal_field import *
 from .tracing import *
+from .currentpotential import *
 from .wireframefield import *
 from .selffield import *
 from .magnetic_axis_helpers import *
@@ -16,6 +17,7 @@
     biotsavart.__all__
     + boozermagneticfield.__all__
     + coil.__all__
+    + currentpotential.__all__
     + coilset.__all__
     + magneticfield.__all__
     + magneticfieldclasses.__all__

src/simsopt/field/coil.py

@@ -559,14 +559,15 @@ def coils_to_vtk(coils, filename, close=False, extra_data=None):
         ppl = np.asarray([c.gamma().shape[0]+1 for c in curves])
     else:
         ppl = np.asarray([c.gamma().shape[0] for c in curves])
+    ppl_cumsum = np.concatenate([[0], np.cumsum(ppl)])
 
     # get the current data, which is the same at every point on a given coil
     contig = np.ascontiguousarray
     pointData = {}
     data = np.concatenate([i*np.ones((ppl[i], )) for i in range(len(curves))])
     coil_data = np.zeros(data.shape)
     for i in range(len(currents)):
-        coil_data[i * ppl[i]: (i + 1) * ppl[i]] = currents[i]
+        coil_data[ppl_cumsum[i]:ppl_cumsum[i+1]] = currents[i]
     coil_data = np.ascontiguousarray(coil_data)
     pointData['I'] = coil_data
     pointData['I_mag'] = contig(np.abs(coil_data))
@@ -593,20 +594,20 @@ def coils_to_vtk(coils, filename, close=False, extra_data=None):
             if close:
                 coil_force_temp = np.vstack((coil_force_temp, coil_force_temp[0, :]))
                 coil_torque_temp = np.vstack((coil_torque_temp, coil_torque_temp[0, :]))
-            coil_forces[i * ppl[i]: (i + 1) * ppl[i], :] = coil_force_temp
-            coil_torques[i * ppl[i]: (i + 1) * ppl[i], :] = coil_torque_temp
+            coil_forces[ppl_cumsum[i]:ppl_cumsum[i+1], :] = coil_force_temp
+            coil_torques[ppl_cumsum[i]:ppl_cumsum[i+1], :] = coil_torque_temp
 
         # copy force and torque data over to pointwise data on a coil curve
         coil_data = np.zeros((data.shape[0], 3))
         for i in range(len(coils)):
-            coil_data[i * ppl[i]: (i + 1) * ppl[i], :] = net_forces[i, :]
+            coil_data[ppl_cumsum[i]:ppl_cumsum[i+1], :] = net_forces[i, :]
         coil_data = np.ascontiguousarray(coil_data)
         pointData['NetForces'] = (contig(coil_data[:, 0]),
                                     contig(coil_data[:, 1]),
                                     contig(coil_data[:, 2]))
         coil_data = np.zeros((data.shape[0], 3))
         for i in range(len(coils)):
-            coil_data[i * ppl[i]: (i + 1) * ppl[i], :] = net_torques[i, :]
+            coil_data[ppl_cumsum[i]:ppl_cumsum[i+1], :] = net_torques[i, :]
         coil_data = np.ascontiguousarray(coil_data)
 
         # Add pointwise force and torque data to the dictionary