Add cuttable solid layer support by approx-infinity · Pull Request #378 · fusion-energy/paramak

approx-infinity · 2026-05-26T21:47:39Z

Summary

This PR introduces several improvements to the plasma-layer handling:

Adds LayerType.SOLID(cut=True) support to split layers across the plasma boundary with stable naming (layer_N.1, layer_N.2)
Adds validation for missing inner SOLID layers before the plasma.
Adds comprehensive tests for cut validation behavior.
Ensure cut only happens when layers exist on both sides of the plasma.

Motivation

Most of the time the Inner Board and Outer Board layers of a tokamak have different material compositions. Till now we need to create a continuous layer which can contain a single material composition. So, this PR improves SOLID layer handling by enabling deterministic splitting of layers while validating invalid configurations early.

1. Add `LayerType.SOLID(cut=True)` support

model = paramak.tokamak_from_plasma(
    radial_build = [
        (LayerType.GAP, 10),
        (LayerType.SOLID, 20),     # layer_1
        (LayerType.SOLID, 20),     # layer_2
        (LayerType.GAP, 5),
        (LayerType.SOLID(cut=True), 20),     # layer_4.1
        (LayerType.GAP, 5),
        (LayerType.SOLID, 10),     # layer_3
        (LayerType.GAP, 10),
        (LayerType.PLASMA, 100),     # plasma
        (LayerType.GAP, 10),
        (LayerType.SOLID, 20),     # layer_3
        (LayerType.GAP, 5),
        (LayerType.SOLID, 40),     # layer_4.2
    ],
    triangularity=0.55,
    rotation_angle=180,
)

print(model.names())

output:
['layer_1', 'layer_2', 'layer_3', 'layer_4.1', 'layer_4.2', 'plasma']

FreeCAD

New behavior
SOLID layers can now be split across the plasma boundary.

Example naming:

layer_4.1
layer_4.2

Validation behavior
Splitting only occurs when geometry exists on both sides of the plasma.
LayerType.SOLID(cut=True) can be called once or twice(both side of plasma) but result will be same.

…single side of plasma

Removed the _copy_metadata method and split_solids method from the Assembly class. I have found a better solution and will add it in next PR.

Removed the test for splitting solids.

shimwell · 2026-05-27T16:48:07Z

Hi @approx-infinity, thanks for putting this together — the PR description and screenshots are really clear, and I can see the use case (different materials inboard vs. outboard).

While I look into the best way to add this upstream, I wanted to mention that splitting a matched layer across the plasma boundary can also be done as a post-process on the returned cq.Assembly, without any changes to paramak itself. It might give a bit more flexibility than the radial-build approach — you pick which layers to split, where to split them, and the cut shape isn't tied to the layer API.

Something like:

import cadquery as cq
import paramak
from paramak.utils import sum_up_to_plasma

radial_build = [
    (paramak.LayerType.GAP, 10),
    (paramak.LayerType.SOLID, 30),
    # ...
    (paramak.LayerType.PLASMA, 300),
    (paramak.LayerType.GAP, 60),
    (paramak.LayerType.SOLID, 20),
    (paramak.LayerType.SOLID, 120),   # the layer you want split
    (paramak.LayerType.SOLID, 10),
]

model = paramak.tokamak_from_plasma(radial_build=radial_build, ...)

# Cylinder at the plasma's inboard edge — natural inboard/outboard boundary
inboard_edge_r = sum_up_to_plasma(radial_build)
boundary = cq.Workplane("XY").circle(inboard_edge_r).extrude(5000, both=True)

target = "layer_7"   # whichever .names() shows
layer = next(p[0] for p in model if p[1].split('/')[-1] == target)

inboard  = layer.intersect(boundary)
outboard = layer.cut(boundary)

model = model.remove(target)
model.add(inboard,  name=f"{target}_inboard")
model.add(outboard, name=f"{target}_outboard")

Would something like this cover your use case in the meantime? Happy to keep discussing the right shape of the upstream API.

approx-infinity · 2026-05-27T18:24:08Z

Thanks for the clarification.

Yes, doing it as a post-process on the returned cq.Assembly should work for my current use case as well, and it actually gives more flexibility as you mentioned.

I will not continue modifying the current PR, since the post-process already covers the use case.

shimwell · 2026-05-27T20:29:16Z

ok super i shall close this pr and make an example of this post processing for the docs. I will just keep it open for now till i have that example ready

approx-infinity · 2026-05-27T21:05:44Z

inboard_edge_r = sum_up_to_plasma(radial_build)

boundary = cq.Workplane("XY").circle(inboard_edge_r).extrude(13000, both=True)
boundary = boundary.val()

targets = ["layer_2", "layer_3", "layer_4", "layer_5", "layer_6"]

new_model = []

for p in eu_demo:
    name = p[1].split("/")[-1]

    if name in targets:
        layer = p[0]

        inboard = layer.intersect(boundary)
        outboard = layer.cut(boundary)

        new_model.append((inboard, f"{name}_inboard"))
        new_model.append((outboard, f"{name}_outboard"))
    else:
        new_model.append(p)

eu_demo = new_model

assembly = cq.Assembly()
for item in eu_demo:
    shape = item[0]
    name = item[1] if isinstance(item[1], str) else str(item[1])
    # strip path separators if present
    name = name.split("/")[-1]
    assembly.add(shape, name=name)

assembly.export("EU_DEMO_Symmetry_180.step")
print("Saved as EU_DEMO_Symmetry_180.step")
print(assembly.name)

this outputs:

(openmc) hridoy@approx:~/Programs/openmc/Fusion/EU_DEMO$ python EU_DEMO_D_symmetry.py
/home/hridoy/Programs/openmc/Fusion/paramak/src/paramak/workplanes/toroidal_field_coil_princeton_d.py:30: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  asol = integrate.odeint(solvr, [z_0, 0], a_R)
Saved as EU_DEMO_Symmetry_180.step
c6142c74-5a0a-11f1-83f7-78465c886f03

I also tried print(assembly.name()) and print(assembly.names()) but they don't work.

FreeCAD

So, here you can see that the same problem like #325, #327, #377 occurs. I applied split_solids() method from #379 to get separate solids with names
assembly = cq.Assembly()

for item in eu_demo:
    shape = item[0]
    name = item[1] if isinstance(item[1], str) else str(item[1])
    # strip path separators if present
    name = name.split("/")[-1]
    assembly.add(shape, name=name)

assembly = assembly.split_solids()

but it shows a error:

(openmc) hridoy@approx:~/Programs/openmc/Fusion/EU_DEMO$ python EU_DEMO_D_symmetry.py
/home/hridoy/Programs/openmc/Fusion/paramak/src/paramak/workplanes/toroidal_field_coil_princeton_d.py:30: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  asol = integrate.odeint(solvr, [z_0, 0], a_R)
Traceback (most recent call last):
  File "/home/hridoy/Programs/openmc/Fusion/EU_DEMO/EU_DEMO_D_symmetry.py", line 128, in <module>
    assembly = assembly.split_solids()
               ^^^^^^^^^^^^^^^^^^^^^
  File "/home/hridoy/miniconda3/envs/openmc/lib/python3.13/site-packages/cadquery/assembly.py", line 820, in __getattr__
    raise AttributeError(f"{name} is not an attribute of {self}")
AttributeError: split_solids is not an attribute of <cadquery.assembly.Assembly object at 0x76defab3e7b0>

So, if we cut solids at the post-processing step then getting the parts name is again confusing.
With LayerType.SOLID(cut=True) it is possible to get stable names through print(model.names()), with or without applying the split_solids(). So, I am thinking to add a position parameter beside cut, like LayerType.SOLID(cut=True, cut_at=N) sothat user gets more flexibility.

Should I apply this?

approx-infinity · 2026-05-27T21:13:04Z

Ok now i got it!

from paramak.assemblies.assembly import Assembly
assembly = Assembly()
for item in eu_demo:
    shape = item[0]
    name = item[1] if isinstance(item[1], str) else str(item[1])
    # strip path separators if present
    name = name.split("/")[-1]
    assembly.add(shape, name=name)

assembly = assembly.split_solids()

assembly.export("EU_DEMO_Symmetry_180.step")
print("Saved as EU_DEMO_Symmetry_180.step")
print(assembly.names())

outputs:
['toroidal_field_coil_1_1', 'toroidal_field_coil_1_2', 'toroidal_field_coil_1_3', 'toroidal_field_coil_1_4', 'toroidal_field_coil_1_5', 'poloidal_field_coil_2', 'poloidal_field_coil_case_3', 'poloidal_field_coil_4', 'poloidal_field_coil_case_5', 'poloidal_field_coil_6', 'poloidal_field_coil_case_7', 'poloidal_field_coil_8', 'poloidal_field_coil_case_9', 'poloidal_field_coil_10', 'poloidal_field_coil_case_11', 'poloidal_field_coil_12', 'poloidal_field_coil_case_13', 'divertor_lower_1_1', 'divertor_lower_1_2', 'layer_1', 'layer_2_inboard', 'layer_2_outboard_1', 'layer_2_outboard_2', 'layer_3_inboard', 'layer_3_outboard_1', 'layer_3_outboard_2', 'layer_4_inboard', 'layer_4_outboard_1', 'layer_4_outboard_2', 'layer_5_inboard', 'layer_5_outboard_1', 'layer_5_outboard_2', 'layer_6_inboard', 'layer_6_outboard_1', 'layer_6_outboard_2', 'layer_7', 'layer_8', 'layer_9', 'layer_10']

FreeCAD

So, we need to use Assembly from paramak.

approx-infinity · 2026-05-27T21:44:57Z

inboard_edge_r = sum_up_to_plasma(radial_build) + 2000

boundary = cq.Workplane("XY").circle(inboard_edge_r).extrude(15000, both=True)
boundary = boundary.val()

targets = ["layer_2", "layer_3", "layer_4", "layer_5", "layer_6"]

new_model = []

for p in eu_demo:
    name = p[1].split("/")[-1]

    if name in targets:
        layer = p[0]

        inboard = layer.intersect(boundary)
        outboard = layer.cut(boundary)

        new_model.append((inboard, f"{name}_inboard"))
        new_model.append((outboard, f"{name}_outboard"))
    else:
        new_model.append(p)

eu_demo = new_model

assembly = Assembly()
for item in eu_demo:
    shape = item[0]
    name = item[1] if isinstance(item[1], str) else str(item[1])
    # strip path separators if present
    name = name.split("/")[-1]
    assembly.add(shape, name=name)

assembly = assembly.split_solids()

assembly.export("EU_DEMO_Symmetry_180.step")
print("Saved as EU_DEMO_Symmetry_180.step")
print(assembly.names())

material_tags = ["tf_coil1", "tf_coil2", "tf_coil3", "tf_coil4", "tf_coil5",
                 "pf_coil1", "pf_coil_case1", "pf_coil2", "pf_coil_case2", "pf_coil3", "pf_coil_case3", "pf_coil4", "pf_coil_case4", "pf_coil5", "pf_coil_case5", "pf_coil6", "pf_coil_case6",
                 "divertor_lower", "divertor_upper",
                 "cs_coil",
                 "w_armor",
                 "ib_first_wall", "ob_first_wall",
                 "ib_breeding_zone", "ob_breeding_zone",
                 "ib_pbli_mainfolds", "ob_pbli_mainfolds",
                 "ib_water_mainfolds", "ob_water_mainfolds",
                 "ib_bss", "ob_bss",
                 "vv_shield_1",
                 "vv",
                 "vv_shield_2"]

eu_demo_model = CadToDagmc()
eu_demo_model.add_cadquery_object(cadquery_object=eu_demo, material_tags=material_tags)

If i follow this work flow, it gives a error:

(openmc) hridoy@approx:~/Programs/openmc/Fusion/EU_DEMO$ python EU_DEMO_D_symmetry.py
/home/hridoy/Programs/openmc/Fusion/paramak/src/paramak/workplanes/toroidal_field_coil_princeton_d.py:30: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  asol = integrate.odeint(solvr, [z_0, 0], a_R)
Saved as EU_DEMO_Symmetry_180.step
['toroidal_field_coil_1_1', 'toroidal_field_coil_1_2', 'toroidal_field_coil_1_3', 'toroidal_field_coil_1_4', 'toroidal_field_coil_1_5', 'poloidal_field_coil_2', 'poloidal_field_coil_case_3', 'poloidal_field_coil_4', 'poloidal_field_coil_case_5', 'poloidal_field_coil_6', 'poloidal_field_coil_case_7', 'poloidal_field_coil_8', 'poloidal_field_coil_case_9', 'poloidal_field_coil_10', 'poloidal_field_coil_case_11', 'poloidal_field_coil_12', 'poloidal_field_coil_case_13', 'divertor_lower_1_1', 'divertor_lower_1_2', 'layer_1', 'layer_2_inboard', 'layer_2_outboard', 'layer_3_inboard', 'layer_3_outboard', 'layer_4_inboard', 'layer_4_outboard', 'layer_5_inboard', 'layer_5_outboard', 'layer_6_inboard', 'layer_6_outboard', 'layer_7', 'layer_8', 'layer_9', 'layer_10']
Traceback (most recent call last):
  File "/home/hridoy/Programs/openmc/Fusion/EU_DEMO/EU_DEMO_D_symmetry.py", line 150, in <module>
    eu_demo_model.add_cadquery_object(cadquery_object=eu_demo, material_tags=material_tags)
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/hridoy/miniconda3/envs/openmc/lib/python3.13/site-packages/cad_to_dagmc/core.py", line 1340, in add_cadquery_object
    iterable_solids = cadquery_compound.val().Solids()
                      ^^^^^^^^^^^^^^^^^^^^^
AttributeError: 'list' object has no attribute 'val'

So, I need to convert the list back into a CadQuery compound before passing it to CadToDagmc.

inboard_edge_r = sum_up_to_plasma(radial_build) + 2000

boundary = cq.Workplane("XY").circle(inboard_edge_r).extrude(15000, both=True)
boundary = boundary.val()

targets = ["layer_2", "layer_3", "layer_4", "layer_5", "layer_6"]

new_model = []

for p in eu_demo:
    name = p[1].split("/")[-1]

    if name in targets:
        layer = p[0]

        inboard = layer.intersect(boundary)
        outboard = layer.cut(boundary)

        new_model.append((inboard, f"{name}_inboard"))
        new_model.append((outboard, f"{name}_outboard"))
    else:
        new_model.append(p)

eu_demo = new_model

assembly = Assembly()
for item in eu_demo:
    shape = item[0]
    name = item[1] if isinstance(item[1], str) else str(item[1])
    # strip path separators if present
    name = name.split("/")[-1]
    assembly.add(shape, name=name)

assembly = assembly.split_solids()

assembly.export("EU_DEMO_Symmetry_180.step")
print("Saved as EU_DEMO_Symmetry_180.step")
print(assembly.names())

material_tags = ["tf_coil1", "tf_coil2", "tf_coil3", "tf_coil4", "tf_coil5",
                 "pf_coil1", "pf_coil_case1", "pf_coil2", "pf_coil_case2", "pf_coil3", "pf_coil_case3", "pf_coil4", "pf_coil_case4", "pf_coil5", "pf_coil_case5", "pf_coil6", "pf_coil_case6",
                 "divertor_lower", "divertor_upper",
                 "cs_coil",
                 "w_armor",
                 "ib_first_wall", "ob_first_wall",
                 "ib_breeding_zone", "ob_breeding_zone",
                 "ib_pbli_mainfolds", "ob_pbli_mainfolds",
                 "ib_water_mainfolds", "ob_water_mainfolds",
                 "ib_bss", "ob_bss",
                 "vv_shield_1",
                 "vv",
                 "vv_shield_2"]

shapes_only = [item[0] for item in eu_demo]

# Build a compound from all solids
compound = cq.Compound.makeCompound(
    [s.val() if hasattr(s, 'val') else s for s in shapes_only]
)

eu_demo_model = CadToDagmc()
eu_demo_model.add_cadquery_object(cadquery_object=compound, material_tags=material_tags)

And now this works as expected,

(openmc) hridoy@approx:~/Programs/openmc/Fusion/EU_DEMO$ python EU_DEMO_D_symmetry.py
/home/hridoy/Programs/openmc/Fusion/paramak/src/paramak/workplanes/toroidal_field_coil_princeton_d.py:30: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  asol = integrate.odeint(solvr, [z_0, 0], a_R)
Saved as EU_DEMO_Symmetry_180.step
['toroidal_field_coil_1_1', 'toroidal_field_coil_1_2', 'toroidal_field_coil_1_3', 'toroidal_field_coil_1_4', 'toroidal_field_coil_1_5', 'poloidal_field_coil_2', 'poloidal_field_coil_case_3', 'poloidal_field_coil_4', 'poloidal_field_coil_case_5', 'poloidal_field_coil_6', 'poloidal_field_coil_case_7', 'poloidal_field_coil_8', 'poloidal_field_coil_case_9', 'poloidal_field_coil_10', 'poloidal_field_coil_case_11', 'poloidal_field_coil_12', 'poloidal_field_coil_case_13', 'divertor_lower_1_1', 'divertor_lower_1_2', 'layer_1', 'layer_2_inboard', 'layer_2_outboard', 'layer_3_inboard', 'layer_3_outboard', 'layer_4_inboard', 'layer_4_outboard', 'layer_5_inboard', 'layer_5_outboard', 'layer_6_inboard', 'layer_6_outboard', 'layer_7', 'layer_8', 'layer_9', 'layer_10']
Using meshing backend: gmsh
Imprinting assembly for mesh generation
volumes [(3, 1), (3, 2), (3, 3), (3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9), (3, 10), (3, 11), (3, 12), (3, 13), (3, 14), (3, 15), (3, 16), (3, 17), (3, 18), (3, 19), (3, 20), (3, 21), (3, 22), (3, 23), (3, 24), (3, 25), (3, 26), (3, 27), (3, 28), (3, 29), (3, 30), (3, 31), (3, 32), (3, 33), (3, 34)]
/home/hridoy/miniconda3/envs/openmc/lib/python3.13/site-packages/cad_to_dagmc/core.py:806: UserWarning: set_size for w_armor is 50.0 which is below min_mesh_size of 70.0. The mesh size will be clamped to 70.0. Try reducing min_mesh_size to encompass the set_size value.
  warnings.warn(
Boundaries for volume 21: [(0, 179), (0, 180), (0, 181), (0, 182), (0, 209), (0, 210), (0, 211), (0, 212), (0, 213), (0, 214), (0, 215), (0, 216)]
Boundaries for volume 22: [(0, 185), (0, 186), (0, 195), (0, 196), (0, 210), (0, 211), (0, 214), (0, 215)]
Set mesh size 50.0 for boundary (0, 179)
Set mesh size 50.0 for boundary (0, 180)
Set mesh size 50.0 for boundary (0, 181)
Set mesh size 50.0 for boundary (0, 182)
Set mesh size 50.0 for boundary (0, 209)
Set mesh size 50.0 for boundary (0, 210)
Set mesh size 50.0 for boundary (0, 211)
Set mesh size 50.0 for boundary (0, 212)
Set mesh size 50.0 for boundary (0, 213)
Set mesh size 50.0 for boundary (0, 214)
Set mesh size 50.0 for boundary (0, 215)
Set mesh size 50.0 for boundary (0, 216)
Set mesh size 70.0 for boundary (0, 185)
Set mesh size 70.0 for boundary (0, 186)

Now the main concern is, we need do lots of conversions in-between. Though it didn't take much time(3-4 minutes) on my PC to start generating mesh.

approx-infinity added 7 commits May 26, 2026 22:45

Add optional cut=True support for solid layers in tokamak builds

876d9be

Fix the extra sublayer creation problem for cutting layers esists on …

32dbc4a

…single side of plasma

Add split_solids feature

b71397a

Add test for the split_solid() feature

ae301cb

Remove split_solids(), _copy_metadata() method from Assembly class

5137d93

Removed the _copy_metadata method and split_solids method from the Assembly class. I have found a better solution and will add it in next PR.

Correcting mistakes

8ea717f

Remove split_solids test

f07f7b5

Removed the test for splitting solids.

approx-infinity changed the title ~~Add split_solids() to Assembly, add cuttable solid layer support~~ Add cuttable solid layer support May 26, 2026

shimwell force-pushed the main branch from 2c16964 to 5764e76 Compare May 29, 2026 14:26

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add cuttable solid layer support#378

Add cuttable solid layer support#378
approx-infinity wants to merge 7 commits into
fusion-energy:mainfrom
approx-infinity:cut_support

approx-infinity commented May 26, 2026 •

edited

Loading

Uh oh!

shimwell commented May 27, 2026

Uh oh!

approx-infinity commented May 27, 2026 •

edited

Loading

Uh oh!

shimwell commented May 27, 2026

Uh oh!

approx-infinity commented May 27, 2026

Uh oh!

approx-infinity commented May 27, 2026

Uh oh!

approx-infinity commented May 27, 2026 •

edited

Loading

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

approx-infinity commented May 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Motivation

1. Add LayerType.SOLID(cut=True) support

Uh oh!

shimwell commented May 27, 2026

Uh oh!

approx-infinity commented May 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

shimwell commented May 27, 2026

Uh oh!

approx-infinity commented May 27, 2026

Uh oh!

approx-infinity commented May 27, 2026

Uh oh!

approx-infinity commented May 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

approx-infinity commented May 26, 2026 •

edited

Loading

1. Add `LayerType.SOLID(cut=True)` support

approx-infinity commented May 27, 2026 •

edited

Loading

approx-infinity commented May 27, 2026 •

edited

Loading