Add sim files

5 files changed, 632 insertions, 4 deletions

diff --git a/coil_parasitics.py b/coil_parasitics.py
new file mode 100644
index 0000000..b57e1f2
--- /dev/null
+++ b/coil_parasitics.py
@@ -0,0 +1,149 @@
+#!/usr/bin/env python3
+
+from pathlib import Path
+import multiprocessing
+import re
+import tempfile
+import subprocess
+import fnmatch
+import shutil
+import numpy as np
+
+from pyelmer import elmer
+import click
+from scipy import constants
+
+def enumerate_mesh_bodies(msh_file):
+    with open(msh_file, 'r') as f:
+        for line in f:
+            if line.startswith('$PhysicalNames'):
+                break
+        else:
+            raise ValueError('No physcial bodies found in mesh file.')
+
+        _num_names = next(f)
+
+        for line in f:
+            if line.startswith('$EndPhysicalNames'):
+                break
+            
+            dim, _, line = line.strip().partition(' ')
+            tag, _, name = line.partition(' ')
+            yield name.strip().strip('"'), (int(dim), int(tag))
+
+INPUT_EXT_MAP = {
+        '.grd': 1,
+        '.mesh*': 2,
+        '.ep': 3,
+        '.ansys': 4,
+        '.inp': 5,
+        '.fil': 6,
+        '.FDNEUT': 7,
+        '.unv': 8,
+        '.mphtxt': 9,
+        '.dat': 10,
+        '.node': 11,
+        '.ele': 11,
+        '.mesh': 12,
+        '.msh': 14,
+        '.ep.i': 15,
+        '.2dm': 16}
+
+OUTPUT_EXT_MAP = {
+        '.grd': 1,
+        '.mesh*': 2,
+        '.ep': 3,
+        '.msh': 4,
+        '.vtu': 5}
+
+def elmer_grid(infile, outfile=None, intype=None, outtype=None, cwd=None, **kwargs):
+    infile = Path(infile)
+    if outfile is not None:
+        outfile = Path(outfile)
+
+    if intype is None:
+        intype = str(INPUT_EXT_MAP[infile.suffix])
+
+    if outtype is None:
+        if outfile is not None and outfile.suffix:
+            outtype = str(OUTPUT_EXT_MAP[outfile.suffix])
+        else:
+            outtype = '2'
+
+    if outfile is not None:
+        kwargs['out'] = str(outfile)
+
+    args = ['ElmerGrid', intype, outtype, infile]
+    for key, value in kwargs.items():
+        args.append(f'-{key}')
+        if isinstance(value, (tuple, list)):
+            args.extend(str(v) for v in value)
+        else:
+            args.append(str(value))
+    subprocess.run(args, cwd=cwd)
+
+def elmer_solver(cwd):
+    subprocess.run(['ElmerSolver'], cwd=cwd)
+
+
+@click.command()
+@click.option('-d', '--sim-dir', type=click.Path(dir_okay=True, file_okay=False, path_type=Path))
+@click.argument('mesh_file', type=click.Path(dir_okay=False, path_type=Path))
+def run_simulation(mesh_file, sim_dir):
+    physical = dict(enumerate_mesh_bodies(mesh_file))
+    if sim_dir is not None:
+        sim_dir = Path(sim_dir)
+        sim_dir.mkdir(exist_ok=True)
+
+    sim = elmer.load_simulation('3D_steady', 'coil_parasitics_sim.yml')
+    mesh_dir = '.'
+    mesh_fn = 'mesh'
+    sim.header['Mesh DB'] = f'"{mesh_dir}" "{mesh_fn}"'
+    sim.constants.update({
+        'Permittivity of Vacuum': str(constants.epsilon_0),
+        'Gravity(4)': f'0 -1 0 {constants.g}',
+        'Boltzmann Constant': str(constants.Boltzmann),
+        'Unit Charge': str(constants.elementary_charge)})
+
+    air = elmer.load_material('air', sim, 'coil_parasitics_materials.yml')
+    ro4003c = elmer.load_material('ro4003c', sim, 'coil_parasitics_materials.yml')
+
+    solver_electrostatic = elmer.load_solver('Electrostatics_Capacitance', sim, 'coil_parasitics_solvers.yml')
+    solver_electrostatic.data['Potential Difference'] = '1.0'
+    eqn = elmer.Equation(sim, 'main', [solver_electrostatic])
+
+    bdy_sub = elmer.Body(sim, 'substrate', [physical['substrate'][1]])
+    bdy_sub.material = ro4003c
+    bdy_sub.equation = eqn
+
+    bdy_ab = elmer.Body(sim, 'airbox', [physical['airbox'][1]])
+    bdy_ab.material = air
+    bdy_ab.equation = eqn
+
+    # boundaries
+    for name, identity in physical.items():
+        if (m := re.fullmatch(r'trace([0-9]+)', name)):
+            num = int(m.group(1))
+
+            bndry_m2 = elmer.Boundary(sim, name, [identity[1]])
+            bndry_m2.data['Capacitance Body'] = str(num)
+
+    boundary_airbox = elmer.Boundary(sim, 'FarField', [physical['airbox_surface'][1]])
+    boundary_airbox.data['Electric Infinity BC'] = 'True'
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        if sim_dir:
+            tmpdir = str(sim_dir)
+
+        sim.write_startinfo(tmpdir)
+        sim.write_sif(tmpdir)
+        # Convert mesh from gmsh to elemer formats. Also scale it from 1 unit = 1 mm to 1 unit = 1 m (SI units)
+        elmer_grid(mesh_file.name, 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3])
+        elmer_solver(tmpdir)
+        
+        capacitance_matrix = np.loadtxt(tmpdir / 'capacitance.txt')
+
+
+
+if __name__ == '__main__':
+    run_simulation()
diff --git a/coil_parasitics_materials.yml b/coil_parasitics_materials.yml
new file mode 100644
index 0000000..ecb49b7
--- /dev/null
+++ b/coil_parasitics_materials.yml
@@ -0,0 +1,76 @@
+air:
+  Density: 1.1885  # 20°C
+  Electric Conductivity: 0.0 
+  Heat Capacity: 1006.4  # 20°C
+  Heat Conductivity: 0.025873  # 20°C
+  Relative Permeability: 1
+  Relative Permittivity: 1
+ro4003c:
+  Density: 1790  # 23°C
+  Relative Permeability: 1
+  Relative Permittivity: 3.55  
+ideal:
+  Relative Permittivity: 1 
+copper_inductor:
+  Density: 8960.0  # 20°C
+  Electric Conductivity: 0.0  # necessary for 2D
+  Emissivity: 0.012  # 327°C
+  Heat Capacity: 384.4  # interpolated for 20°C
+  Heat Conductivity: 401.0
+  Relative Permeability: 1
+  Relative Permittivity: 1
+copper:
+  Density: 8960.0  # 0°C
+  Electric Conductivity: 32300000  # 200°C
+  Emissivity: 0.012  # 327°C
+  Heat Capacity: 415.0  # 200°C
+  Heat Conductivity: 401.0  # 0°C
+  Relative Permeability: 1
+  Relative Permittivity: 1
+graphite_CZ3-R6300:  # crucible
+  Density: 1730.0
+  Electric Conductivity: 58800
+  Emissivity: 0.81  # 205°C
+  Heat Capacity: 1237.0
+  Heat Conductivity: 65  # 20°C
+  Relative Permeability: 1
+  Relative Permittivity: 1
+graphite_FU8957:  # heater
+  Density: 1750.0
+  Emissivity: 0.81  # 250°C
+  Heat Capacity: 1237.0
+  Heat Conductivity: 105  # averaged over different given values
+  Relative Permeability: 1
+  Relative Permittivity: 1
+steel_1.4541:
+  Density: 7900.0  # 20°C
+  Electric Conductivity: 1370
+  Emissivity: 0.111  # 200°C
+  Heat Capacity:  470.0 # 20°C
+  Heat Conductivity: 15.0 # 20°C
+  Relative Permeability: 1
+  Relative Permittivity: 1
+tin_liquid:
+  Density: 6980.0
+  Electric Conductivity: 2080000
+  Emissivity: 0.064  # set equal to solid
+  Heat Capacity: 252.7
+  Heat Conductivity: 29.0
+  Relative Permeability: 1
+  Relative Permittivity: 1
+  Liquid: 'Logical True'
+tin_solid:
+  Density: 7179.0
+  Electric Conductivity: 4380000
+  Emissivity: 0.064
+  Heat Capacity: 244.0
+  Heat Conductivity: 60.0
+  Relative Permeability: 1
+  Relative Permittivity: 1
+  Solid: 'Logical True'
+  Melting Point: 505
+  Latent Heat: 59600
+water:
+  Density: 1000.0
+  Heat Capacity: 4182.0
+  Heat Conductivity: 0.6
diff --git a/coil_parasitics_sim.yml b/coil_parasitics_sim.yml
new file mode 100644
index 0000000..f65dd35
--- /dev/null
+++ b/coil_parasitics_sim.yml
@@ -0,0 +1,50 @@
+2D_steady:
+  Max Output Level: 4
+  Coordinate System: Cartesian 2D
+  Simulation Type: Steady state
+  Steady State Max Iterations: 10
+  
+3D_steady:
+  Max Output Level: 5
+  Coordinate System: Cartesian
+  Coordinate Mapping(3): 1 2 3
+  Simulation Type: Steady state
+  Steady State Max Iterations: 1
+  Output Intervals: 1
+  Timestepping Method: BDF
+  BDF Order: 1
+  Solver Input File: case.sif
+  Post File: case.vtu
+  Output File: case.result
+  
+2D_transient:
+  Max Output Level: 4
+  Coordinate System: Cartesian 2D
+  Simulation Type: Transient
+  Steady State Max Iterations: 10
+  Output File: case.result
+  Output Intervals: 10
+  Timestep Sizes: 0.1
+  Timestep Intervals: 100
+  Timestepping Method: BDF
+  BDF Order: 1
+
+
+axi-symmetric_steady:
+  Max Output Level: 4
+  Coordinate System: Axi Symmetric
+  Simulation Type: Steady state
+  Steady State Max Iterations: 10
+  Output File: case.result
+
+axi-symmetric_transient:
+  Max Output Level: 4
+  Coordinate System: Axi Symmetric
+  Simulation Type: Transient
+  Steady State Max Iterations: 10
+  Output File: case.result
+  Output Intervals: 10
+  Timestep Sizes: 0.1
+  Timestep Intervals: 100
+  Timestepping Method: BDF
+  BDF Order: 1
diff --git a/coil_parasitics_solvers.yml b/coil_parasitics_solvers.yml
new file mode 100644
index 0000000..93a935c
--- /dev/null
+++ b/coil_parasitics_solvers.yml
@@ -0,0 +1,203 @@
+Electrostatics_Capacitance:
+  Equation: Electrostatics
+  Calculate Electric Field: True
+  Calculate Capacitance Matrix: True
+  Capacitance Matrix Filename: capacitance.txt
+  Procedure: '"StatElecSolve" "StatElecSolver"'
+  Variable: Potential
+  Calculate Electric Energy: True
+  Exec Solver: Always
+  Stabilize: True
+  Bubbles: False
+  Lumped Mass Matrix: False
+  Optimize Bandwidth: True
+  Steady State Convergence Tolerance: 1.0e-5
+  Nonlinear System Convergence Tolerance: 1.0e-7
+  Nonlinear System Max Iterations: 20
+  Nonlinear System Newton After Iterations: 3
+  Nonlinear System Newton After Tolerance: 1.0e-3
+  Nonlinear System Relaxation Factor: 1
+  Linear System Solver: Iterative
+  Linear System Iterative Method: BiCGStab
+  Linear System Max Iterations: 500
+  Linear System Convergence Tolerance: 1.0e-10
+  BiCGstabl polynomial degree: 2
+  Linear System Preconditioning: ILU0
+  Linear System ILUT Tolerance: 1.0e-3
+  Linear System Abort Not Converged: False
+  Linear System Residual Output: 10
+  Linear System Precondition Recompute: 1
+Electrostatics:
+  Equation: Electrostatics
+  Calculate Electric Field: True
+  Procedure: '"StatElecSolve" "StatElecSolver"'
+  Variable: Potential
+  Calculate Electric Energy: True
+  Exec Solver: Always
+  Stabilize: True
+  Bubbles: False
+  Lumped Mass Matrix: False
+  Optimize Bandwidth: True
+  Steady State Convergence Tolerance: 1.0e-5
+  Nonlinear System Convergence Tolerance: 1.0e-7
+  Nonlinear System Max Iterations: 20
+  Nonlinear System Newton After Iterations: 3
+  Nonlinear System Newton After Tolerance: 1.0e-3
+  Nonlinear System Relaxation Factor: 1
+  Linear System Solver: Iterative
+  Linear System Iterative Method: BiCGStab
+  Linear System Max Iterations: 500
+  Linear System Convergence Tolerance: 1.0e-10
+  BiCGstabl polynomial degree: 2
+  Linear System Preconditioning: ILU0
+  Linear System ILUT Tolerance: 1.0e-3
+  Linear System Abort Not Converged: False
+  Linear System Residual Output: 10
+  Linear System Precondition Recompute: 1
+ThermoElectricSolver:
+  Equation: ThermoElectric
+  Procedure: '"ThermoElectricSolver" "ThermoElectricSolver"'
+  Variable: POT[Temperature:1 Potential:1]
+  Element: '"p:1"'
+  Calculate Loads: True
+  Exec Solver: Always
+  Nonlinear System Convergence Tolerance: 1.0e-6
+  Nonlinear System Max Iterations: 100
+  Nonlinear System Newton After Iterations : 1
+  Nonlinear System Newton After Tolerance: 1e-9
+  Linear System Solver: '"Iterative"'
+  Linear System Iterative Method: BicgstabL
+  Bicgstabl Polynomial Degree: 2
+  Linear System Max Iterations: 200
+  Linear System Residual Output: 40
+  Linear System Preconditioning: Ilu
+  Linear System Convergence Tolerance: 1e-8
+  Steady State Convergence Tolerance: 1e-6
+HeatSolver:
+  Equation: HeatSolver
+  Procedure: '"HeatSolve" "HeatSolver"'
+  Variable: '"Temperature"'
+  Variable Dofs: 1
+  Calculate Loads: True
+  Exec Solver: Always
+  Nonlinear System Convergence Tolerance: 1.0e-6
+  Nonlinear System Max Iterations: 1000
+  Nonlinear System Relaxation Factor: 0.7
+  Steady State Convergence Tolerance: 1.0e-6
+  Stabilize: True  # Necessary in convection-dominated systems
+  Optimize Bandwidth: True
+  Linear System Solver: Iterative  
+  Linear System Iterative Method: BiCGStab
+  Linear System Max Iterations: 1000
+  Linear System Preconditioning: ILU
+  Linear System Precondition Recompute: 1
+  Linear System Convergence Tolerance: 1.0e-8
+  Linear System Abort Not Converged: True
+  Linear System Residual Output: 1
+  Smart Heater Control After Tolerance: 1.0e-4
+MagnetoDynamics2DHarmonic:
+  Equation: MagnetoDynamics2DHarmonic
+  Procedure: '"MagnetoDynamics2D" "MagnetoDynamics2DHarmonic"'
+  Variable: 'Potential[Potential Re:1 Potential Im:1]'
+  Variable Dofs: 2
+  Exec Solver: Always
+  Nonlinear System Convergence Tolerance: 1.0e-6
+  Nonlinear System Max Iterations: 1000
+  Nonlinear System Relaxation Factor: 0.7
+  Steady State Convergence Tolerance: 1.0e-6
+  Optimize Bandwidth: True
+  Linear System Solver: Iterative  
+  Linear System Iterative Method: BiCGStab
+  Linear System Max Iterations: 1000
+  Linear System Preconditioning: ILU
+  Linear System Precondition Recompute: 1
+  Linear System Convergence Tolerance: 1.0e-8
+  Linear System Abort Not Converged: True
+  Linear System Residual Output: 1
+MagnetoDynamicsCalcFields:
+  Equation: MagnetoDynamicsCalcFields
+  Procedure: '"MagnetoDynamics" "MagnetoDynamicsCalcFields"'
+  Potential Variable: Potential
+  Angular Frequency: 8.48e4
+  Calculate Joule Heating: True
+  Calculate Magnetic Field Strength: True
+  Calculate Electric Field: True
+  Exec Solver: Always
+  Calculate Nodal Fields: Logical False
+  Calculate Elemental Fields: Logical True
+StatMagSolver:
+  Equation: StatMagSolver
+  Procedure: '"StatMagSolve" "StatMagSolver"'
+  Variable: Potential
+  Variable DOFs: 2
+  Calculate Joule Heating: 'Logical True'
+  Calculate Magnetic Flux: 'Logical True'
+  Exec Solver: Always
+  Nonlinear System Convergence Tolerance: 1.0e-6
+  Nonlinear System Max Iterations: 1000
+  Nonlinear System Relaxation Factor: 0.7
+  Steady State Convergence Tolerance: 1.0e-6
+  Optimize Bandwidth: True
+  Linear System Solver: Iterative  
+  Linear System Iterative Method: BiCGStab
+  Linear System Max Iterations: 1000
+  Linear System Preconditioning: ILU
+  Linear System Precondition Recompute: 1
+  Linear System Convergence Tolerance: 1.0e-8
+  Linear System Abort Not Converged: True
+  Linear System Residual Output: 1
+SaveMaterials: 
+  Exec Solver: 'After timestep'
+  Procedure: 'File "SaveData" "SaveMaterials"'
+  Parameter 1: 'String "Heat Conductivity"'
+ResultOutputSolver:
+  Exec Solver: 'After timestep'
+  Equation: ResultOutputSolver
+  Procedure: '"ResultOutputSolve" "ResultOutputSolver"'
+  VTU Format: True
+  Save Geometry Ids: 'Logical True'
+FluxSolver:
+  Exec Solver: 'After timestep'
+  Equation: 'Flux Solver'
+  Procedure: '"FluxSolver" "FluxSolver"'
+  Calculate Grad: 'Logical True'
+  Calculate Flux: 'Logical True'
+  Target Variable: 'String "Temperature"'
+  Flux Coefficient: 'String "Heat Conductivity"'
+  Exec Solver: Always
+  Nonlinear System Convergence Tolerance: 1.0e-6
+  Nonlinear System Max Iterations: 1000
+  Nonlinear System Relaxation Factor: 0.7
+  Steady State Convergence Tolerance: 1.0e-6
+  Optimize Bandwidth: True
+  Linear System Solver: Iterative  
+  Linear System Iterative Method: BiCGStab
+  Linear System Max Iterations: 1000
+  Linear System Preconditioning: ILU
+  Linear System Precondition Recompute: 1
+  Linear System Convergence Tolerance: 1.0e-8
+  Linear System Abort Not Converged: True
+  Linear System Residual Output: 1
+SaveScalars:
+  Exec Solver: 'After timestep'
+  Equation: SaveScalars
+  Procedure: '"SaveData" "SaveScalars"'
+  Filename: '"boundary_scalars.dat"'
+  Output Directory: './results'
+  Operator 1: 'boundary sum'
+  Variable 1: 'Temperature Loads'
+  Operator 2: 'diffusive flux'
+  Variable 2: Temperature
+  Coefficient 2: 'Heat Conductivity'
+SaveLine:
+  Exec Solver: 'After timestep'
+  Equation: '"SaveLine"'
+  Procedure: '"SaveData" "SaveLine"'
+  Filename: '"boundary_lines.dat"'
+  Output Directory: './results'
+  Variable 1: Temperature Loads
+SteadyPhaseChange:
+  Equation: SteadyPhaseChange
+  Variable: '"Phase Surface"'
+  Procedure: '"SteadyPhaseChange" "SteadyPhaseChange"'
+  Internal Mesh Movement: 'Logical True'
diff --git a/twisted_coil_gen_twolayer.py b/twisted_coil_gen_twolayer.py
index 77fea7a..475fe16 100644
--- a/twisted_coil_gen_twolayer.py
+++ b/twisted_coil_gen_twolayer.py
@@ -7,7 +7,10 @@ import os
 from math import *
 from pathlib import Path
 from itertools import cycle
+
 from scipy.constants import mu_0
+import numpy as np
+import click
 import matplotlib as mpl
 
 from gerbonara.cad.kicad import pcb as kicad_pcb
@@ -16,7 +19,7 @@ from gerbonara.cad.kicad import graphical_primitives as kicad_gr
 from gerbonara.cad.kicad import primitives as kicad_pr
 from gerbonara.utils import Tag
 from gerbonara import graphic_primitives as gp
-import click
+from gerbonara import graphic_objects as go
 
 
 __version__ = '1.0.0'
@@ -64,7 +67,7 @@ def traces_to_gmsh(traces, mesh_out, bbox, model_name='gerbonara_board', log=Tru
     trace_tags = {}
     trace_ends = set()
     render_cache = {}
-    for i, tr in enumerate(traces):
+    for i, tr in enumerate(traces, start=1):
         layer = tr[1].layer
         z0 = 0 if layer == 'F.Cu' else -(board_thickness+copper_thickness)
 
@@ -141,6 +144,141 @@ def traces_to_gmsh(traces, mesh_out, bbox, model_name='gerbonara_board', log=Tru
     gmsh.write(str(mesh_out))
 
 
+def traces_to_gmsh_mag(traces, mesh_out, bbox, model_name='gerbonara_board', log=True, copper_thickness=0.035, board_thickness=0.8, air_box_margin=5.0):
+    import gmsh
+    occ = gmsh.model.occ
+    eps = 1e-6
+
+    board_thickness -= 2*copper_thickness
+
+    gmsh.initialize()
+    gmsh.model.add('gerbonara_board')
+    if log:
+        gmsh.logger.start()
+
+    trace_tags = []
+    trace_ends = set()
+    render_cache = {}
+    for i, tr in enumerate(traces, start=1):
+        layer = tr[1].layer
+        z0 = 0 if layer == 'F.Cu' else -(board_thickness+copper_thickness)
+
+        objs = [obj
+                 for elem in tr
+                 for obj in elem.render(cache=render_cache)]
+
+        tags = []
+        for ob in objs:
+            if isinstance(ob, go.Line):
+                length = dist((ob.x1, ob.y1), (ob.x2, ob.y2))
+                w = ob.aperture.equivalent_width('mm')
+                box_tag = occ.addBox(0, -w/2, 0, length, w, copper_thickness)
+                angle = atan2(ob.y2 - ob.y1, ob.x2 - ob.x1)
+                occ.rotate([(3, box_tag)], 0, 0, 0, 0, 0, 1, angle)
+                occ.translate([(3, box_tag)], ob.x1, ob.y1, z0)
+                tags.append(box_tag)
+
+                for x, y in ((ob.x1, ob.y1), (ob.x2, ob.y2)):
+                    disc_id = (round(x, 3), round(y, 3), round(z0, 3), round(w, 3))
+                    if disc_id  in trace_ends:
+                        continue
+
+                    trace_ends.add(disc_id)
+                    cylinder_tag = occ.addCylinder(x, y, z0, 0, 0, copper_thickness, w/2)
+                    tags.append(cylinder_tag)
+            
+        for elem in tr:
+            if isinstance(elem, kicad_pcb.Via):
+                cylinder_tag = occ.addCylinder(elem.at.x, elem.at.y, 0, 0, 0, -board_thickness, elem.drill)
+                tags.append(cylinder_tag)
+
+        print('fusing', tags)
+        tags, tag_map = occ.fuse([(3, tags[0])], [(3, tag) for tag in tags[1:]])
+        print(tags)
+        assert len(tags) == 1
+        (_dim, tag), = tags
+        trace_tags.append(tag)
+
+    print('fusing top-level', trace_tags)
+    tags, tag_map = occ.fuse([(3, trace_tags[0])], [(3, tag) for tag in trace_tags[1:]])
+    print(tags)
+    assert len(tags) == 1
+    (_dim, toplevel_tag), = tags
+
+    first_geom = traces[0][0]
+    interface_tag_top = occ.addDisk(first_geom.start.x, first_geom.start.y, 0, first_geom.width, first_geom.width)
+    interface_tag_bottom = occ.addDisk(first_geom.end.x, first_geom.end.y, -board_thickness, first_geom.width, first_geom.width)
+
+    (x1, y1), (x2, y2) = bbox
+    substrate = occ.addBox(x1, y1, -board_thickness, x2-x1, y2-y1, board_thickness)
+
+    x1, y1 = x1-air_box_margin, y1-air_box_margin
+    x2, y2 = x2+air_box_margin, y2+air_box_margin
+    w, d = x2-x1, y2-y1
+    z0 = -board_thickness-air_box_margin
+    ab_h = board_thickness + 2*air_box_margin
+    airbox = occ.addBox(x1, y1, z0, w, d, ab_h)
+
+    print('Cutting airbox')
+    occ.cut([(3, airbox)], [(3, toplevel_tag)], removeObject=True, removeTool=False)
+
+    print('Cutting substrate')
+    occ.cut([(3, substrate)], [(3, toplevel_tag)], removeObject=True, removeTool=False)
+
+    print('Synchronizing')
+    occ.synchronize()
+
+    substrate_physical = gmsh.model.add_physical_group(3, [substrate], name='substrate')
+    airbox_physical = gmsh.model.add_physical_group(3, [airbox], name='airbox')
+    trace_physical_surfaces = [
+            gmsh.model.add_physical_group(2, list(gmsh.model.getAdjacencies(3, tag)[1]), name=f'trace{i}')
+            for i, tag in trace_tags.items()]
+    
+    airbox_adjacent = set(gmsh.model.getAdjacencies(3, airbox)[1])
+    in_bbox = {tag for _dim, tag in gmsh.model.getEntitiesInBoundingBox(x1+eps, y1+eps, z0+eps, x1+w-eps, y1+d-eps, z0+ab_h-eps, dim=3)}
+    airbox_physical_surface = gmsh.model.add_physical_group(2, list(airbox_adjacent - in_bbox), name='airbox_surface')
+    
+    points_airbox_adjacent = set(gmsh.model.getAdjacencies(0, airbox)[1])
+    points_inside = {tag for _dim, tag in gmsh.model.getEntitiesInBoundingBox(x1+eps, y1+eps, z0+eps, x1+w-eps, y1+d-eps, z0+ab_h-eps, dim=0)}
+    gmsh.model.mesh.setSize([(0, tag) for tag in points_airbox_adjacent - points_inside], 10e-3)
+
+    gmsh.option.setNumber('Mesh.MeshSizeFromCurvature', 90)
+    gmsh.option.setNumber('Mesh.Smoothing', 10)
+    gmsh.option.setNumber('Mesh.Algorithm3D', 10)
+    gmsh.option.setNumber('Mesh.MeshSizeMax', 1)
+    gmsh.option.setNumber('General.NumThreads', multiprocessing.cpu_count())
+
+    print('Meshing')
+    gmsh.model.mesh.generate(dim=3)
+    print('Writing')
+    gmsh.write(str(mesh_out))
+
+
+def traces_to_magneticalc(traces, out, pcb_thickness=0.8):
+    coords = []
+    last_x, last_y, last_z = None, None, None
+    def coord(x, y, z):
+        nonlocal coords, last_x, last_y, last_z
+        if (x, y, z) != (last_x, last_y, last_z):
+            coords.append((x, y, z))
+
+    render_cache = {}
+    for tr in traces:
+        z = pcb_thickness if tr[1].layer == 'F.Cu' else 0
+        objs = [obj
+                 for elem in tr
+                 for obj in elem.render(cache=render_cache)
+                 if isinstance(elem, (kicad_pcb.TrackSegment, kicad_pcb.TrackArc))]
+
+        # start / switch layer
+        coord(objs[0].x1, objs[0].y1, z)
+
+        for ob in objs:
+            coord(ob.x2, ob.y2, z)
+
+    np.savetxt(out, np.array(coords) / 10) # magneticalc expects centimeters, not millimeters.
+
+
 class SVGPath:
     def __init__(self, **attrs):
         self.d = ''
@@ -235,16 +373,19 @@ def print_valid_twists(ctx, param, value):
 @click.option('--keepout-zone/--no-keepout-zone', default=True, help='Add a keepout are to the footprint (default: yes)')
 @click.option('--keepout-margin', type=float, default=5, help='Margin between outside of coil and keepout area (mm, default: 5)')
 @click.option('--twists', type=int, default=1, help='Number of twists per revolution. Note that this number must be co-prime to the number of turns. Run with --show-twists to list valid values. (default: 1)')
+@click.option('--circle-segments', type=int, default=64, help='When not using arcs, the number of points to use for arc interpolation per 360 degrees.')
 @click.option('--show-twists', callback=print_valid_twists, expose_value=False, type=int, is_eager=True, help='Calculate and show valid --twists counts for the given number of turns. Takes the number of turns as a value.')
 @click.option('--clearance', type=float, default=None)
 @click.option('--arc-tolerance', type=float, default=0.02)
 @click.option('--mesh-out', type=click.Path(writable=True, dir_okay=False, path_type=Path))
+@click.option('--mag-mesh-out', type=click.Path(writable=True, dir_okay=False, path_type=Path))
+@click.option('--magneticalc-out', type=click.Path(writable=True, dir_okay=False, path_type=Path))
 @click.option('--clipboard/--no-clipboard', help='Use clipboard integration (requires wl-clipboard)')
 @click.option('--counter-clockwise/--clockwise', help='Direction of generated spiral. Default: clockwise when wound from the inside.')
 @click.version_option()
 def generate(outfile, turns, outer_diameter, inner_diameter, via_diameter, via_drill, via_offset, trace_width, clearance,
              footprint_name, layer_pair, twists, clipboard, counter_clockwise, keepout_zone, keepout_margin,
-             arc_tolerance, pcb, mesh_out):
+             arc_tolerance, pcb, mesh_out, magneticalc_out, circle_segments, mag_mesh_out):
     if 'WAYLAND_DISPLAY' in os.environ:
         copy, paste, cliputil = ['wl-copy'], ['wl-paste'], 'xclip'
     else:
@@ -256,6 +397,9 @@ def generate(outfile, turns, outer_diameter, inner_diameter, via_diameter, via_d
     if mesh_out and not pcb:
         raise click.ClickException('--pcb is required when --mesh-out is used.')
 
+    if magneticalc_out and not pcb:
+        raise click.ClickException('--pcb is required when --magneticalc-out is used.')
+
     outer_radius = outer_diameter/2
     inner_radius = inner_diameter/2
     turns_per_layer = turns/2
@@ -466,7 +610,7 @@ def generate(outfile, turns, outer_diameter, inner_diameter, via_diameter, via_d
         end_angle = fold_angle + sweeping_angle
 
         x = inverse[i]*floor(2*sweeping_angle / (2*pi)) * 2*pi
-        (x0, y0), (xn, yn), clen = do_spiral(0, outer_radius, inner_radius, start_angle, fold_angle, (x + start_angle)/total_angle, (x + fold_angle)/total_angle)
+        (x0, y0), (xn, yn), clen = do_spiral(0, outer_radius, inner_radius, start_angle, fold_angle, (x + start_angle)/total_angle, (x + fold_angle)/total_angle, circle_segments)
         do_spiral(1, inner_radius, outer_radius, fold_angle, end_angle, (x + fold_angle)/total_angle, (x + end_angle)/total_angle)
 
         xv, yv = inner_via_ring_radius*cos(fold_angle), inner_via_ring_radius*sin(fold_angle)
@@ -578,6 +722,12 @@ def generate(outfile, turns, outer_diameter, inner_diameter, via_diameter, via_d
         if mesh_out:
             traces_to_gmsh(traces, mesh_out, ((-r, -r), (r, r)))
 
+        if mag_mesh_out:
+            traces_to_gmsh_mag(traces, mesh_out, ((-r, -r), (r, r)))
+
+        if magneticalc_out:
+            traces_to_magneticalc(traces, magneticalc_out)
+
 #        for trace in traces:
 #            print(f'Trace {i}', file=sys.stderr)
 #            print(f'  Length: {len(trace)}', file=sys.stderr)