Как реализовать решатель ограничений для двумерной геометрии?

Я бы попробовал полевой подход, аналогичный этому.

1. Каждый ползунок отодвигает все ползунки

   с силой, масштабируемой расстоянием ^ 2, как будто все они будут иметь электрический заряд одинаковой полярности или пружины, прикрепленные друг к другу.

2. Кроме того, добавьте трение в зависимости от скорости

   на самом деле не имеет значения, воздух v^2 или жидкость v^3

3. внедрить кинематические ограничения

   только для горизонтального и вертикального скольжения это должно быть очень легко.

4. Выполните физическое моделирование и подождите, пока оно не сойдется в стабильное состояние v=~0

   если вы попали в локальный минимум / максимум, немного встряхните все это или расположите все это случайным образом и повторите попытку. Вы также можете сделать это, чтобы получить другое решение.

[Edit4] Пример решателя C++

1. структуры/классы для представления системы ползунков

   Чтобы облегчить последующий код, я не буду поддерживать замкнутые циклы или двойную привязку. Вот почему ползунок i1 (самый правый) ни к чему не привязан (просто обеспечивает силовое поле). Я закончил с этим определением слайдера:

   

   посмотрите источник class _slider для получения дополнительной информации.

2. рендеринг

   Тире-тире означает фиксированный ползунок. Серебряные — горизонтальные, цвета морской волны — вертикальные, а желтые выбираются мышью. Возможно, позже красный цвет будет означать какую-то ошибку/зависание или что-то в целях отладки. Для решателей силового поля я иногда добавляю напряженность поля в виде красно-синей шкалы, но не уверен, буду ли я реализовывать это здесь или нет.

   Для простоты я не буду реализовывать функции масштабирования/панорамирования, поскольку ваши размеры удобны для прямого рендеринга без преобразований.

   

3. выполнить первоначальную настройку

   sliders sys;
   int i0,i1,a0,a1,a2,a3,a4,b1,b2,b3,b4,b5;
   sys.slider_beg();//ia,ib,   x,    y,    a0,    a1,    b0,    b1,_horizontal
   i0=sys.slider_add(-1,-1, 25.0, 25.0,  -5.0, 405.0,   0.0,   0.0, 0);
   a0=sys.slider_add(i0,-1,  0.0,  0.0,   0.0, 400.0,   0.0,   0.0, 1);
   a1=sys.slider_add(i0,-1,  0.0,100.0,   0.0, 400.0,   0.0,   0.0, 1);
   a2=sys.slider_add(i0,-1,  0.0,200.0,   0.0, 400.0,   0.0,   0.0, 1);
   a3=sys.slider_add(i0,-1,  0.0,300.0,   0.0, 400.0,   0.0,   0.0, 1);
   a4=sys.slider_add(i0,-1,  0.0,400.0,   0.0, 400.0,   0.0,   0.0, 1);
   b1=sys.slider_add(a0,a2, 20.0,  0.0,   0.0, 125.0, 125.0, 250.0, 0);
   b2=sys.slider_add(a3,-1, 40.0,  0.0, -70.0,  30.0,   0.0,   0.0, 0);
   b3=sys.slider_add(a1,-1, 60.0,  0.0, -70.0,  30.0,   0.0,   0.0, 0);
   b4=sys.slider_add(a2,-1, 80.0,  0.0, -30.0,  70.0,   0.0,   0.0, 0);
   b5=sys.slider_add(a3,a1,100.0,  0.0,-125.0,   0.0,-125.0,-250.0, 0);
   i1=sys.slider_add(-1,-1,425.0, 25.0,  -5.0, 405.0,   0.0,   0.0, 0);
   sys.slider_end();
   

   Где ia — родительский индекс, а ib — дочерний индекс (сам класс ползунка содержит ib в качестве родителя, но это может привести к путанице при инициализации, поскольку вам нужно будет связать элемент, который еще не существует, поэтому преобразование ib обрабатывается в функции sys.add) . sys — это класс, содержащий все это, а sys.add просто добавьте к нему новый ползунок и верните его индекс, считая с нуля. x,y - это относительное положение по отношению к родителю.

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

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

4. взаимодействие с мышью

   простое движение ползунка для отладки и настройки начальных значений настройки. И или обработка застрявших дел. Вам нужно обрабатывать события мыши, выберите ближайший ползунок, если он еще не редактируется. И если кнопка мыши нажата, переместите выбранный ползунок в положение мыши...

5. физическое ограничение/взаимодействие

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

   Использование - это немного больше кода. Сначала сохраните фактическую позицию для обновленного ползунка. Затем обновите ползунок до новой позиции/состояния. После этого, если ограничения не выполняются, остановите фактическую скорость ползунка и восстановите его исходное положение.

   Это будет немного медленнее, но я слишком ленив, чтобы кодировать полное обновление ограничений (этот код может стать очень сложным...).

   Я различаю 2 взаимодействия параллельно и перпендикулярно. Параллель идет прямо. Но перпендикуляр - это взаимодействие между краем ползунка и перпендикулярными ползунками рядом с ним, не включая уже пересекающиеся ползунки (привязанные a, b или просто пересекающиеся) в исходном состоянии. Поэтому я создал список пересекающихся ползунков (ic) при запуске, которые будут игнорироваться для этого взаимодействия.

6. физическое моделирование

   Подойдет простая физика Ньютона-Д'Аламбера для нерелятивистских скоростей. Просто на каждой итерации устанавливайте ускорения ax,ay на напряженность поля и трения.

7. решатель полей

   Это набор правил/уравнений для установки ускорений симуляции для каждого ползунка, чтобы сходиться к решению. Я остановился на электростатической втягивающей силе F = -Q/r^2 и линейном демпфировании скорости. Также реализованы ограничители абсолютной скорости и ускорения, чтобы избежать числовых проблем.

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

Вот полный код класса C++/VCL для этого:

//---------------------------------------------------------------------------
//--- Sliders solver ver: 1.01 ----------------------------------------------
//---------------------------------------------------------------------------
#ifndef _sliders_h
#define _sliders_h
//---------------------------------------------------------------------------
#include <math.h>
#include "list.h"   // linear dynamic array template List<T> similar to std::vector
//---------------------------------------------------------------------------
const double _slider_w   =   3.00;  // [px] slider half width (for rendering)
const double _slider_gap =   4.00;  // [px] min gap between sliders (for colisions)
const double _acc_limit=   100.00;  // [px/s^2]
const double _vel_limit=   100.00;  // [px/s]
const double _friction =     0.90;  // [-]
const double _charge   =250000.00;  // [px^3/s^2]
//---------------------------------------------------------------------------
class _slider   // one slider (helper class)
    {
public:
    // properties
    double x,y;             // actual relative pos
    bool _horizontal;       // orientation
    double a0,a1;           // slider vertexes 0 is anchor point
    double b0,b1;           // anchor zone for another slider
    int ia;                 // -1 for fixed or index of parrent slider
    int ib;                 // -1 or index of parrent slider
    // computed
    List<int> ic;           // list of slider indexes to ignore for perpendicular constraints
    double a,b;             // force field affected part
    double X,Y;             // actual absolute position
    double vx,vy,ax,ay;     // actual relative vel,acc
    // temp
    int flag;               // temp flag for simulation
    double x0,x1;           // temp variables for solver
    // constructors (can ignore this)
    _slider()           {}
    _slider(_slider& a) { *this=a; }
    ~_slider()          {}
    _slider* operator = (const _slider *a) { *this=*a; return this; }
    //_slider* operator = (const _slider &a) { ...copy... return this; }
    };
//---------------------------------------------------------------------------
class sliders   // whole slider system main class
    {
public:
    List<_slider> slider;           // list of sliders

    double vel_max;                 // max abs velocity of sliders for solver precision control
    double charge;                  // actual charge of sliders for solve()
    int    mode;                    // actual solution precision control mode

    // constructors (can ignore this)
    sliders();
    sliders(sliders& a) { *this=a; }
    ~sliders()          {}
    sliders* operator = (const sliders *a) { *this=*a; return this; }
    //sliders* operator = (const sliders &a) { ...copy... return this; }

    // VCL window API variables (can ignore this)
    double mx0,my0,mx1,my1; // last and actual mouse position
    TShiftState sh0,sh1;    // last and actual mouse buttons and control keys state
    int sel;

    // API (this is important stuff)
    void slider_beg(){ slider.num=0; }  // clear slider list
    int  slider_add(int ia,int ib,double x,double y,double a0,double a1,double b0,double b1,bool _h); // add slider to list
    void slider_end();              // compute slider parameters
    bool constraints(int ix);       // return true if constraints hit
    void positions();               // recompute absolute positions
    void update(double dt);         // update physics simulation with time step dt [sec]
    void solve(bool _init=false);   // set sliders accelerations to solve this
    void stop();                    // stop all movements
    // VCL window API for interaction with GUI (can ignore this)
    void mouse(int x,int y,TShiftState sh);
    void draw(TCanvas *scr);
    };
//---------------------------------------------------------------------------
sliders::sliders()
    {
    mx0=0.0; my0=0.0;
    mx1=0.0; my1=0.0;
    sel=-1;
    }
//---------------------------------------------------------------------------
int sliders::slider_add(int ia,int ib,double x,double y,double a0,double a1,double b0,double b1,bool _h)
    {
    _slider s; double q;
    if (a0>a1) { q=a0; a0=a1; a1=q; }
    if (b0>b1) { q=b0; b0=b1; b1=q; }
    s.x=x; s.vx=0.0; s.ax=0.0;
    s.y=y; s.vy=0.0; s.ay=0.0;
    s.ia=ia; s.a0=a0; s.a1=a1;
    s.ib=-1; s.b0=b0; s.b1=b1;
    s.ic.num=0;
    if ((ib>=0)&&(ib<slider.num)) slider[ib].ib=slider.num;
    s._horizontal=_h;
    s.a=a0; // min
    if (s.a>a1) s.a=a1;
    if (s.a>b0) s.a=b0;
    if (s.a>b1) s.a=b1;
    s.b=a0; // max
    if (s.b<a1) s.b=a1;
    if (s.b<b0) s.b=b0;
    if (s.b<b1) s.b=b1;
    slider.add(s);
    return slider.num-1;
    }
//---------------------------------------------------------------------------
void sliders::slider_end()
    {
    int i,j;
    double a0,a1,b0,b1,x0,x1,w=_slider_gap;
    _slider *si,*sj;
    positions();
    // detect intersecting sliders and add them to propriet ic ignore list
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
     for (sj=si+1   ,j=i+1;j<slider.num;j++,sj++)
      if (si->_horizontal!=sj->_horizontal)
        {
        if (si->_horizontal)
            {
            a0=si->X+si->a; a1=sj->X-w;
            b0=si->X+si->b; b1=sj->X+w;
            x0=si->Y;       x1=sj->Y;
            }
        else{
            a0=si->Y+si->a; a1=sj->Y-w;
            b0=si->Y+si->b; b1=sj->Y+w;
            x0=si->X;       x1=sj->X;
            }
        if (((a0<=b1)&&(b0>=a1))||((a1<=b0)&&(b1>=a0)))
         if ((x0>x1+sj->a-w)&&(x0<x1+sj->b+w))
            {
            si->ic.add(j);
            sj->ic.add(i);
            }
        }
    }
//---------------------------------------------------------------------------
bool sliders::constraints(int ix)
    {
    int i,j;
    double a0,a1,b0,b1,x0,x1,x,w=_slider_gap;
    _slider *si,*sj,*sa,*sb,*s;
    s=slider.dat+ix;
    // check parallel neighbors overlapp
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
     if ((i!=ix)&&(si->_horizontal==s->_horizontal))
        {
        if (s->_horizontal)
            {
            a0=s->X+s->a; a1=si->X+si->a;
            b0=s->X+s->b; b1=si->X+si->b;
            x0=s->Y;      x1=si->Y;
            }
        else{
            a0=s->Y+s->a; a1=si->Y+si->a;
            b0=s->Y+s->b; b1=si->Y+si->b;
            x0=s->X;      x1=si->X;
            }
        if (((a0<=b1)&&(b0>=a1))||((a1<=b0)&&(b1>=a0)))
            {
            if ((i<ix)&&(x0<x1+w)) return true;
            if ((i>ix)&&(x0>x1-w)) return true;
            }
        }
    // check perpendicular neighbors overlapp
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
     if ((i!=ix)&&(si->_horizontal!=s->_horizontal))
        {
        // skip ignored sliders for this
        for (j=0;j<s->ic.num;j++)
         if (s->ic[j]==i) { j=-1; break; }
          if (j<0) continue;
        if (s->_horizontal)
            {
            a0=s->X+s->a; a1=si->X-w;
            b0=s->X+s->b; b1=si->X+w;
            x0=s->Y;      x1=si->Y;
            }
        else{
            a0=s->Y+s->a; a1=si->Y-w;
            b0=s->Y+s->b; b1=si->Y+w;
            x0=s->X;      x1=si->X;
            }
        if (((a0<=b1)&&(b0>=a1))||((a1<=b0)&&(b1>=a0)))
         if ((x0>x1+si->a-w)&&(x0<x1+si->b+w))
          return true;
        }
    // conflict a anchor area of parent?
    if (s->ia>=0)
        {
        si=slider.dat+s->ia;
        if (s->_horizontal)
            {
            x0=si->Y+si->a0;
            x1=si->Y+si->a1;
            x=s->Y;
            }
        else{
            x0=si->X+si->a0;
            x1=si->X+si->a1;
            x=s->X;
            }
        if (x<x0+w) return true;
        if (x>x1-w) return true;
        }
    // conflict b anchor area of parent?
    if (s->ib>=0)
        {
        si=slider.dat+s->ib;
        if (si->_horizontal)
            {
            x0=si->X+si->b0;
            x1=si->X+si->b1;
            x=s->X;
            }
        else{
            x0=si->Y+si->b0;
            x1=si->Y+si->b1;
            x=s->Y;
            }
        if (x<x0+w) return true;
        if (x>x1-w) return true;
        }
    // conflict b anchor area with childs?
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
     if ((i!=ix)&&(si->ib==ix))
        {
        if (s->_horizontal)
            {
            x0=s->X+s->b0;
            x1=s->X+s->b1;
            x=si->X;
            }
        else{
            x0=s->Y+s->b0;
            x1=s->Y+s->b1;
            x=si->Y;
            }
        if (x<x0+w) return true;
        if (x>x1-w) return true;
        }

    // check childs too
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
     if ((i!=ix)&&(si->ia==ix))
      if (constraints(i)) return true;
    return false;
    }
//---------------------------------------------------------------------------
void sliders::positions()
    {
    int i,e;
    _slider *si,*sa;
    // set flag = uncomputed
    for (si=slider.dat,i=0;i<slider.num;i++,si++) si->flag=0;
    // iterate until all sliders are computed
    for (e=1;e;)
     for (e=0,si=slider.dat,i=0;i<slider.num;i++,si++)
      if (!si->flag)
        {
        // fixed
        if (si->ia<0)
            {
            si->X=si->x;
            si->Y=si->y;
            si->flag=1;
            continue;
            }
        // a anchored
        sa=slider.dat+si->ia;
        if (sa->flag)
            {
            si->X=sa->X+si->x;
            si->Y=sa->Y+si->y;
            si->flag=1;
            continue;
            }
        e=1; // not finished yet
        }
    }
//---------------------------------------------------------------------------
void sliders::update(double dt)
    {
    int i;
    _slider *si,*sa;
    double x,X;
    // D'Lamnbert integration
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
     if (si->_horizontal)
        {
        x=si->y; si->vy+=si->ay*dt;     // vel = Integral(acc*dt)
                 si->vy*=_friction;     // friction k*vel
        X=si->Y; si->y +=si->vy*dt;     // pos = Integral(vel*dt)
        positions();                    // recompute childs
        if ((si->ia<0)||(constraints(i))) // if fixed or constraint hit (stop and restore original position)
            {
            si->vy=0.0;
            si->y =x;
            si->Y =X;
            positions();                // recompute childs
            }
        }
    else{
        x=si->x; si->vx+=si->ax*dt;     // vel = Integral(acc*dt)
                 si->vx*=_friction;     // friction k*vel
        X=si->X; si->x +=si->vx*dt;     // pos = Integral(vel*dt)
        positions();                    // recompute childs
        if ((si->ia<0)||(constraints(i))) // if fixed or constraint hit (stop and restore original position)
            {
            si->vx=0.0;
            si->x =x;
            si->X =X;
            positions();                // recompute childs
            }
        }
    }
//---------------------------------------------------------------------------
void sliders::solve(bool _init)
    {
    int i,j,k;
    double a0,a1,b0,b1,x0,x1;
    _slider *si,*sj,*sa;
    // init solution
    if (_init)
        {
        mode=0;
        charge=_charge;
        }
    // clear accelerations and compute actual max velocity
    vel_max=0.0;
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
        {
        si->ax=0.0;
        si->ay=0.0;
        x0=fabs(si->vx); if (vel_max<x0) vel_max=x0;
        x0=fabs(si->vy); if (vel_max<x0) vel_max=x0;
        }
    // precision control of solver
    if ((mode==0)&&(vel_max>25.0)) { mode++; }                  // wait until speed raises
    if ((mode==1)&&(vel_max<10.0)) { mode++; charge*=0.10; }    // scale down forces to lower jitter
    if ((mode==2)&&(vel_max< 1.0)) { mode++; charge*=0.10; }    // scale down forces to lower jitter
    if ((mode==3)&&(vel_max< 0.1)) { mode++; charge =0.00; stop(); } // solution found
    // set x0 as 1D vector to closest parallel neighbor before and x1 after
    for (si=slider.dat,i=0;i<slider.num;i++,si++) { si->x0=0.0; si->x1=0.0; }
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
     for (sj=si+1   ,j=i+1;j<slider.num;j++,sj++)
      if (si->_horizontal==sj->_horizontal)
        {
        // longer side interaction
        if (si->_horizontal)
            {
            a0=si->X+si->a; a1=sj->X+sj->a;
            b0=si->X+si->b; b1=sj->X+sj->b;
            x0=si->Y;       x1=sj->Y;
            }
        else{
            a0=si->Y+si->a; a1=sj->Y+sj->a;
            b0=si->Y+si->b; b1=sj->Y+sj->b;
            x0=si->X;       x1=sj->X;
            }
        if (((a0<=b1)&&(b0>=a1))||((a1<=b0)&&(b1>=a0)))
            {
            x0=x1-x0;
            if ((si->ia>=0)&&(x0<0.0)&&((fabs(si->x0)<_slider_gap)||(fabs(si->x0)>fabs(x0)))) si->x0=-x0;
            if ((si->ia>=0)&&(x0>0.0)&&((fabs(si->x1)<_slider_gap)||(fabs(si->x1)>fabs(x0)))) si->x1=-x0;
            if ((sj->ia>=0)&&(x0<0.0)&&((fabs(sj->x0)<_slider_gap)||(fabs(sj->x0)>fabs(x0)))) sj->x0=+x0;
            if ((sj->ia>=0)&&(x0>0.0)&&((fabs(sj->x1)<_slider_gap)||(fabs(sj->x1)>fabs(x0)))) sj->x1=+x0;
            }
        // shorter side interaction
        if (si->_horizontal)
            {
            a0=si->Y-_slider_gap; a1=sj->Y+_slider_gap;
            b0=si->Y+_slider_gap; b1=sj->Y+_slider_gap;
            x0=si->X;             x1=sj->X;
            }
        else{
            a0=si->X-_slider_gap; a1=sj->X+_slider_gap;
            b0=si->X+_slider_gap; b1=sj->X+_slider_gap;
            x0=si->Y;             x1=sj->Y;
            }
        if (((a0<=b1)&&(b0>=a1))||((a1<=b0)&&(b1>=a0)))
            {
            if (x0<x1) { x0+=si->b; x1+=sj->a; }
            else       { x0+=si->a; x1+=sj->b; }
            x0=x1-x0;
            if (si->ia>=0)
                {
                sa=slider.dat+si->ia;
                if ((sa->ia>=0)&&(x0<0.0)&&((fabs(sa->x0)<_slider_gap)||(fabs(sa->x0)>fabs(x0)))) sa->x0=-x0;
                if ((sa->ia>=0)&&(x0>0.0)&&((fabs(sa->x1)<_slider_gap)||(fabs(sa->x1)>fabs(x0)))) sa->x1=-x0;
                }
            if (sj->ia>=0)
                {
                sa=slider.dat+sj->ia;
                if ((sa->ia>=0)&&(x0<0.0)&&((fabs(sa->x0)<_slider_gap)||(fabs(sa->x0)>fabs(x0)))) sa->x0=+x0;
                if ((sa->ia>=0)&&(x0>0.0)&&((fabs(sa->x1)<_slider_gap)||(fabs(sa->x1)>fabs(x0)))) sa->x1=+x0;
                }
            }
        }
    // set x0 as 1D vector to closest perpendicular neighbor before and x1 after
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
     for (sj=si+1   ,j=i+1;j<slider.num;j++,sj++)
      if (si->_horizontal!=sj->_horizontal)
        {
        // skip ignored sliders for this
        for (k=0;k<si->ic.num;k++)
         if (si->ic[k]==j) { k=-1; break; }
          if (k<0) continue;
        if (si->_horizontal)
            {
            a0=si->X+si->a; a1=sj->X-_slider_w;
            b0=si->X+si->b; b1=sj->X+_slider_w;
            x0=si->Y;
            }
        else{
            a0=si->Y+si->a; a1=sj->Y-_slider_w;
            b0=si->Y+si->b; b1=sj->Y+_slider_w;
            x0=si->X;
            }
        if (((a0<=b1)&&(b0>=a1))||((a1<=b0)&&(b1>=a0)))
            {
            if (si->_horizontal)
                {
                a1=sj->Y+sj->a;
                b1=sj->Y+sj->b;
                }
            else{
                a1=sj->X+sj->a;
                b1=sj->X+sj->b;
                }
            a1-=x0; b1-=x0;
            if (fabs(a1)<fabs(b1)) x0=-a1; else x0=-b1;
            if ((si->ia>=0)&&(x0<0.0)&&((fabs(si->x0)<_slider_gap)||(fabs(si->x0)>fabs(x0)))) si->x0=+x0;
            if ((si->ia>=0)&&(x0>0.0)&&((fabs(si->x1)<_slider_gap)||(fabs(si->x1)>fabs(x0)))) si->x1=+x0;
            if (sj->ia<0) continue;
            sa=slider.dat+sj->ia;
            if ((sa->ia>=0)&&(x0<0.0)&&((fabs(sa->x0)<_slider_gap)||(fabs(sa->x0)>fabs(x0)))) sa->x0=-x0;
            if ((sa->ia>=0)&&(x0>0.0)&&((fabs(sa->x1)<_slider_gap)||(fabs(sa->x1)>fabs(x0)))) sa->x1=-x0;
            }
        }
    // convert x0,x1 distances to acceleration
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
        {
        // driving force F = ~ Q / r^2
        if (fabs(si->x0)>1e-10)  x0=charge/(si->x0*si->x0); else x0=0.0; if (si->x0<0.0) x0=-x0;
        if (fabs(si->x1)>1e-10)  x1=charge/(si->x1*si->x1); else x1=0.0; if (si->x1<0.0) x1=-x1;
        a0=x0+x1;
        // limit acc
        if (a0<-_acc_limit) a0=-_acc_limit;
        if (a0>+_acc_limit) a0=+_acc_limit;
        // store parallel acc to correct axis
        if (si->_horizontal) si->ay=a0;
         else                si->ax=a0;
        // limit vel (+/- one iteration overlap)
        if (si->_horizontal) x0=si->vy;
         else                x0=si->vx;
        if (x0<-_vel_limit)  x0=-_vel_limit;
        if (x0>+_vel_limit)  x0=+_vel_limit;
        if (si->_horizontal) si->vy=x0;
         else                si->vx=x0;
        }
    }
//---------------------------------------------------------------------------
void sliders::stop()
    {
    int i;
    _slider *si;
    for (si=slider.dat,i=0;i<slider.num;i++,si++)
        {
        si->vx=0.0;
        si->vy=0.0;
        si->ax=0.0;
        si->ay=0.0;
        }
    }
//---------------------------------------------------------------------------
void sliders::mouse(int x,int y,TShiftState sh)
    {
    int i,q0,q1;
    double d,dd;
    _slider *si;
    // update mouse state
    mx0=mx1; my0=my1; sh0=sh1;
    mx1=x;   my1=y;   sh1=sh;
    // slider movement with left mouse button
    q0=sh0.Contains(ssLeft);
    q1=sh1.Contains(ssLeft);
    if ((sel>=0)&&(q1))
        {
        si=slider.dat+sel;
        // stop simulation for selected slider
        si->vx=0.0;
        si->vy=0.0;
        si->ax=0.0;
        si->ay=0.0;
        // use mouse position instead
        if (si->ia>=0)
            {
            if (si->_horizontal){ d=si->y; dd=si->Y; si->y+=my1-si->Y; si->Y=my1; si->vy=0.0; si->ay=0.0; positions(); if (constraints(sel)) { si->y=d; si->Y=dd; positions(); }}
             else               { d=si->x; dd=si->X; si->x+=mx1-si->X; si->X=mx1; si->vx=0.0; si->ax=0.0; positions(); if (constraints(sel)) { si->x=d; si->X=dd; positions(); }}
            }
        }
    // select slider (if not left mouse button used)
    if (!q1)
     for (sel=-1,d=_slider_w+1.0,si=slider.dat,i=0;i<slider.num;i++,si++)
        {
        dd=_slider_w+1.0;
        if (si->_horizontal){ if ((mx1>=si->X+si->a)&&(mx1<=si->X+si->b)) dd=fabs(my1-si->Y); }
         else               { if ((my1>=si->Y+si->a)&&(my1<=si->Y+si->b)) dd=fabs(mx1-si->X); }
        if ((dd<d)&&(dd<=_slider_w)) { sel=i; d=dd; }
        }
    }
//---------------------------------------------------------------------------
void sliders::draw(TCanvas *scr)
    {
    int i,j,n;
    double w=_slider_w,r,x,y,a0,a1;
    AnsiString txt;
    _slider *s;
    scr->Brush->Style=bsClear;
    #define _line(aa,bb)           \
    if (s->_horizontal)            \
        {                          \
        scr->MoveTo(s->X+aa,s->Y); \
        scr->LineTo(s->X+bb,s->Y); \
        }                          \
    else{                          \
        scr->MoveTo(s->X,s->Y+aa); \
        scr->LineTo(s->X,s->Y+bb); \
        }
    scr->Pen->Color=clSilver;
    scr->Font->Color=clWhite;
    scr->TextOutA(40,40,AnsiString().sprintf("mode %i",mode));
    scr->TextOutA(40,60,AnsiString().sprintf("vel: %.3lf [px/s]",vel_max));
    scr->TextOutA(40,80,AnsiString().sprintf("  Q: %.3lf [px^3/s^2]",charge));
    scr->Font->Color=clYellow;
    for (s=slider.dat,i=0;i<slider.num;i++,s++)
        {
        if (s->_horizontal) scr->Pen->Color=clSilver;
         else               scr->Pen->Color=clAqua;
        if (i==sel)
            {
            scr->Pen->Color=clYellow;
            txt=AnsiString().sprintf(" ix:%i ia:%i ib:%i ic:",sel,s->ia,s->ib);
            for (j=0;j<=s->ic.num;j++) txt+=AnsiString().sprintf(" %i",s->ic[j]);
            scr->TextOutA(40,100,txt);
            scr->TextOutA(40,120,AnsiString().sprintf("pos: %.1lf %.1lf [px]",s->X,s->Y));
            scr->TextOutA(40,140,AnsiString().sprintf("vel: %.3lf %.3lf [px/s]",s->vx,s->vy));
            scr->TextOutA(40,160,AnsiString().sprintf("acc: %.3lf %.3lf [px/s^2]",s->ax,s->ay));
            scr->Pen->Color=clYellow;
            }
        if (s->ia<0) scr->Pen->Style=psDash;
         else        scr->Pen->Style=psSolid;
        // a anchor loop
        x=s->X;
        y=s->Y;
        if (s->ia>=0) scr->Ellipse(x-w,y-w,x+w,y+w);
        // b anchor loop
        r=0.5*fabs(s->b1-s->b0);
        if (s->_horizontal)
            {
            x=s->X+0.5*(s->b0+s->b1);
            y=s->Y;
            scr->RoundRect(x-r,y-w,x+r,y+w,w,w);
            }
        else{
            x=s->X;
            y=s->Y+0.5*(s->b0+s->b1);
            scr->RoundRect(x-w,y-r,x+w,y+r,w,w);
            }
        // a line cutted by a anchor loop
        a0=s->a0; a1=s->a1;
        if ((s->ia>=0)&&(a0<=+w)&&(a1>=-w))
            {
            if (a0<-w) _line(s->a0,-w);
            if (a1>+w) _line( w,s->a1);
            }
        else _line(s->a0,s->a1);
        }
    scr->Font->Color=clDkGray;
    scr->Pen->Style=psSolid;
    scr->Brush->Style=bsSolid;
    #undef _line
    }
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

Вы можете игнорировать материал VCL, это просто API для взаимодействия с окном моего приложения и рендеринга. Самому солверу от него ничего не нужно. Я использовал свой шаблон динамического линейного массива List<T>, так что вот несколько пояснений:

• List<double> xxx; совпадает с double xxx[];
• xxx.add(5); добавляет 5 в конец списка
• xxx[7] элемент массива доступа (безопасный)
• xxx.dat[7] элемент массива доступа (небезопасный, но быстрый прямой доступ)
• xxx.num - фактический используемый размер массива
• xxx.reset() очищает массив и устанавливает xxx.num=0
• xxx.allocate(100) предварительно выделить место для 100 элементов

Использование просто после правильной инициализации из пули #3 следующим образом:

sys.solve(true);
for (;;)
 {
 sys.solve();
 sys.update(0.040); // just time step
 if (sys.mode==4) break; // stop if solution found or stuck
 }

Вместо цикла for я вызываю это по таймеру и перерисовываю окно, чтобы увидеть анимацию:

• Демо нажмите на медленную загрузку под большой пурпурной кнопкой, введите 4-буквенный буквенно-цифровой код, чтобы начать загрузку без регистрации.

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

• Принудительно направленный макет для ограниченных прямоугольных форм