Трёхмерная схема установки кодом

Преимущества программного описания геометрии

Трёхмерная модель установки обычно создаётся вручную в CAD-редакторе; для машиностроения такой подход подходит, а для ускорителя, перестраиваемого каждый сезон, — нет.

Магнитная структура хранится в текстовом файле со списком элементов, снабжённых длинами, градиентами и углами поворота, и этот файл изменяется постоянно: добавляется квадруполь, переносится диагностика, меняется угол, заданный поворотному магниту. Каждое такое изменение означает, что построенная вручную модель устарела, и сотрудник, отвечающий за чертежи, должен вспомнить о необходимости её перестроить. Как правило, этого не происходит.

В библиотеке CadQuery геометрия описывается кодом на Python. Модель размещается в репозитории рядом с расчётом, различия между версиями видны в git diff, сборка чертежа включается в конвейер тестов, запускаемый на каждом коммите, и тело установки строится непосредственно из описания структуры, заданного один раз, а не параллельно ему.

Структура как данные

Рассмотрим участок тракта, содержащий два квадруполя, дрейфы между ними и два поворотных магнита по 45°.

import numpy as np
import cadquery as cq

# тип элемента, длина в метрах, параметр
LATTICE = [
    ("drift", 0.30, None), ("quad", 0.20, +1), ("drift", 0.30, None),
    ("bend",  0.50, np.pi / 4),
    ("drift", 0.30, None), ("quad", 0.20, -1), ("drift", 0.30, None),
    ("bend",  0.50, np.pi / 4),
]

В реальном проекте такой список не задаётся вручную, а читается из файла магнитной структуры, по которому работает расчётный код. Здесь он задан явно, чтобы пример был самодостаточным, однако источником геометрии и источником физики служит один и тот же файл, прочитанный дважды разными программами.

Один тип элемента, одна функция

На каждый тип элемента заводится своя функция-конструктор, знающая о создаваемом ею теле и ничего не знающая о соседях.

def drift(length, _):
    return cq.Workplane("XY").circle(0.02).extrude(length)

def quad(length, sign):
    body = cq.Workplane("XY").rect(0.24, 0.24).extrude(length)
    return body.edges("|Z").fillet(0.03) if sign > 0 else body

def bend(length, angle):
    r = length / angle                       # радиус по длине дуги и углу
    return (cq.Workplane("XY").rect(0.20, 0.16)
              .revolve(np.degrees(angle), (r, 0, 0), (r, 1, 0)))

BUILDERS = {"drift": drift, "quad": quad, "bend": bend}

Тела здесь строятся выдавливанием и вращением плоского контура, а не составлением из примитивов. Дрейф задан окружностью, вытянутой на длину, а поворотный магнит прямоугольником, проведённым вокруг оси на требуемый угол, и геометрия остаётся параметрической: при изменении угла в описании структуры тело перестраивается автоматически.

Словарь BUILDERS в конце заменяет длинную цепочку if, разрастающуюся с каждым новым типом элемента. Добавление секступоля сводится к написанию функции и одной строки в словаре; код расстановки, написанный один раз, изменять не требуется.

Размещение элементов

Элементы расставляются вдоль траектории пучка, изгибающейся на каждом поворотном магните.

Введём две величины: положение начала текущего элемента и направление, задаваемое пучку. Для каждого элемента необходимо знать, куда сместится его конец относительно начала, унаследованного от предыдущего элемента.

def local_step(kind, length, param):
    """Смещение конца элемента относительно начала, в его собственных осях."""
    if kind != "bend":
        return np.array([0.0, 0.0, length]), 0.0
    r = length / param
    return np.array([r * (1 - np.cos(param)), 0.0, r * np.sin(param)]), param

Прямой элемент сдвигает пучок вперёд на отведённую ему длину и направления не меняет. Поворотный ведёт пучок по дуге радиуса \(R = L/\theta\), с концом, смещённым вперёд на \(R\sin\theta\) и вбок на \(R(1-\cos\theta)\), а направление доворачивается на \(\theta\).

Остаётся пройти по описанной структуре, поворачивая локальное смещение на накопленный угол.

def build(lattice):
    assembly, pos, heading = cq.Assembly(), np.zeros(3), 0.0
    for i, (kind, length, param) in enumerate(lattice):
        assembly.add(BUILDERS[kind](length, param), name=f"{kind}{i}",
                     loc=cq.Location(cq.Vector(*pos),
                                     cq.Vector(0, 1, 0), np.degrees(heading)))
        step, turn = local_step(kind, length, param)
        c, s = np.cos(heading), np.sin(heading)
        pos = pos + np.array([c * step[0] + s * step[2], 0.0,
                              -s * step[0] + c * step[2]])
        heading += turn
    return assembly, pos, heading

Порядок трёх действий, выполняемых внутри цикла, существен. Сначала элемент ставится в текущее положение с текущим направлением, затем сдвигается положение, и только после этого доворачивается направление. При перестановке двух последних строк каждый элемент, поставленный после поворотного магнита, окажется смещённым.

Проверка

Геометрию, построенную кодом, легко проверить, не открывая просмотрщик. Два поворота по 45° должны развернуть пучок в точности на 90°, а сумма длин элементов, выписанных в структуре, даёт длину траектории.

acc, pos, heading = build(LATTICE)
print("элементов:", len(acc.children))
bb = acc.toCompound().BoundingBox()
print(f"габарит: X {bb.xlen:.2f} м, Z {bb.zlen:.2f} м")
print(f"конец траектории: X {pos[0]:.3f} м, Z {pos[2]:.3f} м")
print(f"поворот: {np.degrees(heading):.1f}°")
assert bb.zmax >= pos[2] and bb.xmax >= pos[0]   # тела обязаны накрывать траекторию

элементов: 8
габарит: X 1.32 м, Z 2.10 м
конец траектории: X 1.202 м, Z 2.002 м
поворот: 90.0°

Два поворота по 45° дают в точности 90°, полученная длина составляет 2.60 м, а конец траектории лежит внутри габарита. Таким образом, чертёж покрывается тестом.

assert в последней строке требует, чтобы тела элементов накрывали траекторию, по которой они расставлены. В процессе написания этой главы магнит строился вращением вокруг неверной оси, и обнаружено это было не визуально (изображение выглядело правдоподобно), а по тому, что Z-габарит, измеренный по сборке, составил 1.83 м при конце траектории на 2.002 м. Тело оказалось смещённым относительно пучка.

Готовая сборка выгружается в обменный формат, читаемый любым CAD.

cq.exporters.export(acc.toCompound(), "layout.step")

Выводы

Геометрия остаётся таким же продуктом кода, как график. К ней применимы контроль версий, сборка в конвейере и тесты. Модель, собранная из описания структуры, от него не отстаёт. Построенная вручную отстаёт всегда.

Разбиение на конструкторы повторяет декомпозицию, описанную в главе «От скрипта к приложению». Одна функция, отвечающая за одно тело, ничего не знает о соседях; расстановка не знает, из чего собран квадруполь. Граница проведена там, где ожидаются изменения: новые типы элементов появляются часто, а правило расстановки, записанное однажды, практически никогда не меняется.

Проверяемость важнее наглядности изображения. Ошибка, скрытая в трёхмерной модели, визуально почти не обнаруживается. Установка, выведенная на экран, выглядит правдоподобно, даже если один магнит смещён на полметра. Арифметика — сумма углов, длина траектории и габарит — обнаруживает её сразу.

Полезные ссылки

Документация CadQuery, описание операций, применимых к телам.
CQ-editor, редактор с трёхмерным просмотром для отладки.
build123d, развитие той же идеи, переписанное с другим синтаксисом.

Идея строить трёхмерную схему ускорителя непосредственно из файла магнитной структуры принадлежит А. В. Петренко (ИЯФ СО РАН). Пример, приведённый в этой главе, написан заново и устроен иначе.

Разработка и применение программного обеспечения в физических исследованиях