Используем силу уравнений Ньютона и численных методов для моделирования динамики простых плоских мешей в реальном времени! В конце вы сможете моделировать падение ножниц ✂️ как на анимации
Зачем? Это наглядный, простой и интересный метод решения задачи многих тел с простым взаимодействием типа точка-палка-точка, что придется как раз для пятничного развлекательного программирования (recreational programming session!).
Среда разработки
Все примеры показаны в WLJS Notebook, так как автору лень изучать ipywidgets для создания простой интерактивщины, а Marimo слишком медленный.
Блокнот-версия доступна по ссылке
Набор точек
Создадим случайный набор точек:
pts0 = RandomReal[{-1,1}, {10, 2}];
pts1 = pts0; pts2 = pts0;
Graphics[{Rectangle[{-10,-10}, {10,10}], Red, Point[pts0 // Offload]}]
Чтобы понять, где будут эти точки в следующий момент времени, нам нужно решать как минимум, задачу кинематики, а так как в системе присутствуют внешние силы (гравитация), то и задачу динамики тоже. Метод Стёрмера–Верле — это лишь одна из численных схем для дискретизации уравнений движения Ньютона (II з.):
Чем она отличается от всех других? Для неё требуется хранить только текущую и две предыдущие позиции: 
, 
и 
. Как видно: скорость не нужна, что сыграет особую роль позже. Попробуйте выполнить следующий код:
With[{dt = 0.5},
pts0 = 2 pts1 - pts2 + Table[{0,-1}, {Length[pts0]}] Power[dt,2];
pts2 = pts1;
pts1 = pts0;
]Добавим случайные связи между точками (связи — это как твёрдые стержни):
bonds = RandomGraph[Length[pts0]{1, 2}, VertexCoordinates->pts0]
Метод Cтёрмера–Верле особенно удобен при наличии ограничений типа жёстких/мягких связей (как выше). Например, если хотим сохранить фиксированные расстояния между точками согласно графу для этого потребовалась бы нетривиальная работа со скоростью и координатой, но метод Верле оставляет нам только 
Жесткость таких связей задаётся не через дифференциальные уравнения, а специальным алгоритмом, который аппроксимирует реальные потенциалы этих связей. На этот счёт есть целая наука, как не решать диффуры, а обойтись манипуляциями координат. Здесь же представлен наиболее наивный и простой вариант.
Пусть начальная длина "стержня" соединяющего две точки — 
, где
В какой-то момент при движении точек это расстояние будет отличным от того, что было в предыдущий момент времени. Тогда мы можем его нахально скорректировать
где k это параметр задающий жёсткость стержня (связи):
бесконечно жесткая палка
мягкая связь (как будто пружина aka закон Гука)
Реализация на чистых функциях:
getDir[bonds_, pts_] := Map[(pts[[#[[1]]]] - pts[[#[[2]]]]) &, bonds]
applyBond[pairs_, k_:1.0, maxDelta_:0.1][p_] := Module[{pts = p},
MapThread[Function[{edge, length}, With[{
diff = pts[[edge[[1]]]] - pts[[edge[[2]]]]
},
With[{
delta = Min[k ( ((length)/(Norm[diff]+0.001)) - 1), maxDelta] diff
},
pts[[edge[[1]]]] += delta/2.0;
pts[[edge[[2]]]] -= delta/2.0;
]
] ], RandomSample[pairs] // Transpose];
pts
] Зададим функция для динамической визуализации наших точек и стержней
showGeometry[points_, bonds_, opts___] := Graphics[{
Gray, Rectangle[5{-1,-1}, 5{1,1}], Blue,
Table[
With[{i = edge[[1]], j = edge[[2]]},
Line[With[{p = points}, {p[[i]], p[[j]]}]] // Offload
]
, {edge, bonds}],
Red, Point[points // Offload]
}, opts]
SetAttributes[showGeometry, HoldFirst];Нарисуем наши точки со стержнями
showGeometry[pts0, EdgeList @ bonds]
Теперь можно собрать основной цикл анимации:
pts0 = RandomReal[{-1,1}, {10, 2}];
pts1 = pts0; pts2 = pts0;
bonds = EdgeList[RandomGraph[Length[pts0]{1, 3}]];
lengths = Norm /@ getDir[bonds, pts0];
pairs = {bonds, lengths} // Transpose;
Do[With[{\[Delta]t = 0.05},
pts0 = Clip[applyBond[pairs, 1.0][
2 pts1 - pts2 + Table[{0,-1}, {i, Length[pts0]}] Power[\[Delta]t,2]
], {-5,5}];
pts2 = pts1;
pts1 = pts0;
Pause[0.01];
], {i, 200}]
Процедурная генерация мешей
Даже без оптимизации (например, векторизации или OpenCLLink) можно попробовать более сложные конструкции и взаимодействие с мышью пользователя.
Обод
Создаём обод из точек:
rim = Join @@
(Table[#{Sin[x], Cos[x]}, {x,-Pi, Pi, Pi/6}] &/@ {1, 0.9} // Transpose);
rbonds = Graph[
Join[
Table[i<->i+1, {i, 1, Length[rim]-1}],
Table[i<->i+2, {i, 1, Length[rim]-2, 2 }],
Table[i<->i+2, {i, 2, Length[rim]-4, 2 }],
{Length[rim]-2 <-> 2},
{Length[rim]-1 <-> 1}
]
, VertexCoordinates->rim]
rbonds = EdgeList[rbonds];
rpairs = {rbonds, Norm /@ getDir[rbonds, 2.0 rim]} // Transpose;
Начальные значения координат (для каждого шага по времени):
rim0 = 2.0 rim; rim1 = rim0; rim2 = rim0;
rimR = rim0;
mousePos = {100.,0.};Выведем их на экран, подпишемся на "анимацию", а также на движение мыши пользователя
EventHandler[
showGeometry[rimR, rbonds, TransitionType->None, "Controls"->False, Epilog->{
AnimationFrameListener[rimR // Offload, "Event"->"anim"],
Circle[mousePos // Offload, 0.5]
}, PlotRange->{{-5,5}, {-5,5}}]
, {"mousemove" -> Function[xy,
mousePos = xy;
]}]Здесь мы подписываемся на цикл отрисовки браузера. Каждый раз когда браузер заканчивает рисовать кадр (обычно 1/60 сек), он отправляет событие "anim". В это время обработчик на стороне WL рассчитывает новое положение точек и отправляет его обратно браузеру по WebSocket. Важно отметить, что новое событие не будет инициировано до тех пор, пока не отработает старое. Таким образом обеспечивается максимально возможная частота кадров синхронизированная с обновлением экрана в текущей конфигурации системы (с учетом задержки сети, где находится сервер, клиент, производительности CPU, лагов JS и т.п.)
Теперь добавим также три кнопки управления анимацией, включая их обработчиков
Button["Start",
EventHandler["anim", With[{}, Do[With[{\[Delta]t = 0.006},
rim0 = applyBond[rpairs, 1.0, 0.8][
Clip[(2 rim1 - rim2) + Table[
If[Max[Abs[i - mousePos]] < 2.0, {0,-1} - 7.0 (i - mousePos), {0,-1}]
, {i, rim0}] Power[\[Delta]t,2], {-5,5}]
];
rim2 = rim1;
rim1 = rim0;
], {i, 10}];
rimR = rim0;
]&];
rimR = rim0;
]
Button["Stop", EventRemove["anim"]]
Button["Reset", With[{}, rim0 = 2.0 rim; rim1 = rim0; rim2 = rim0;
rimR = rim0;]]Результат можете увидить ниже
Меш из изображения
Было бы здорово нарисовать что-то мышкой, и потом преобразовать в сетку вершин и связей автоматически
EventHandler[InputRaster[mask], (mask = #)&]
Нарисовали ножницы ✂️ хорошо. Теперь триангулируем и преобразуем в граф:
mesh = TriangulateMesh[
ImageMesh[ImageResize[mask , 100] // ColorNegate]
, MaxCellMeasure->{"Area"->100}];
graph = MeshConnectivityGraph[mesh, 0] // IndexGraph
vertices = graph // GraphEmbedding;
edges = (List @@ #) &/@ (graph // EdgeList) ;
vertices = Map[Function[x, x - {0.4,-1.2}], 0.04 vertices];
Сделаем цикл анимации, как и раньше, и используем другие имена переменных
vertices0 = vertices; vertices1 = vertices0; vertices2 = vertices0;
verticesR = vertices0;
vpairs = {edges, Norm /@ getDir[edges, vertices0]} // Transpose;
mousePos = {100.,0.};EventHandler[
showGeometry[verticesR, edges, TransitionType->None, "Controls"->False, Epilog->{
AnimationFrameListener[verticesR // Offload, "Event"->"anim2"],
Circle[mousePos // Offload, 0.5]
}, PlotRange->{{-5,5}, {-5,5}}]
, {"mousemove" -> Function[xy,
mousePos = xy;
]}]Button["Start",
EventHandler["anim2", With[{}, Do[With[{\[Delta]t = 0.003},
vertices0 = applyBond[vpairs, 1.0, 0.8][
Clip[(2 vertices1 - vertices2) + Table[
If[Max[Abs[i - mousePos]] < 2.0, {0,-1} - 7.0 (i - mousePos), {0,-1}]
, {i, vertices0}] Power[\[Delta]t,2], {-5,5}]
];
vertices2 = vertices1;
vertices1 = vertices0;
], {i, 3}];
verticesR = vertices0;
]&];
verticesR = vertices0;
]
Button["Stop", EventRemove["anim2"]]
Button["Reset", With[{}, vertices0 = vertices; vertices1 = vertices0; vertices2 = vertices0;
verticesR = vertices0;]]Ура! У нас получилось ✨
Очевидно скриптовый язык плох в таких численных расчётах, и стоит поместить бОльшую часть логики в выражение Compile. Тогда вся тяжесть будет автоматически перенесена на сгенерированный и динамически подгруженный C-код. Но об этом в другой раз 🧙🏼♂️
