Basic Verlet Integration implementation

*Using the power of Newton's equations and numerics to solve dynamics of arbitary planar meshes in real-time*

Bunch of points

Define a set of random points

pts0 = RandomReal[{-1,1}, {10, 2}];
pts1 = pts0; pts2 = pts0;

Graphics[{Rectangle[{-10,-10}, {10,10}], Red, Point[pts0 // Offload]}]

The **Størmer–Verlet method** is a variant of the discretization of Newton's equations of motion:

x_{i+1} = 2x_i - x_{i-1} + a_i\, \delta t^2

This method only requires storing the current position $x_{i+1}$ and the two previous positions $x_i$ and $x_{i-1}$; velocity is not needed. Try to run this cell below

With[{\[Delta]t = 0.5},
  pts0  = 2 pts1 - pts2 + Table[{0,-1}, {Length[pts0]}] (*SpB[*)Power[\[Delta]t(*|*),(*|*)2](*]SpB*);
  pts2 = pts1;
  pts1 = pts0;
]

We can add random connections to our points (forming bonds) using graphs

bonds = RandomGraph[Length[pts0]{1, 2}, VertexCoordinates->pts0]

The **Størmer–Verlet method** really shines when it comes to solving equations with constraints. In our case, we want to keep the distances between points fixed according to a connection graph. Normally, this would involve approximating motion using a combination of acceleration, velocity, and position.

However, things can become complicated quickly. With the Størmer–Verlet method, we only work with positions $ x_i $. By directly adjusting the positions of two points connected by a bond to maintain the bond length, we effectively account for all forces in a single step.

The stiffness of the bond is not directly in the differential equation It’s introduced through how you enforce the constraint — either strictly (hard) or partially (soft). Let's define the initial length $L=|\\mathbf{d}|$

\\mathbf{d} := \\mathbf{x}_a - \\mathbf{x}_b

To enforce the bond we compute the correction vector

\\Delta \\mathbf{d} := \\mathbf{d} \\cdot k \\Big( \\frac{L}{|\\mathbf{d}| } - 1 \\Big)

and then apply it to each point

\\begin{align*}

\\mathbf{x}_a &= \\mathbf{x}_a + \\frac{1}{2} \\Delta \\mathbf{d} \\\\

\\mathbf{x}_b &= \\mathbf{x}_b - \\frac{1}{2} \\Delta \\mathbf{d} \\\\

\\end{align*}

where $k$ is a stiffness parameter. There are two cases:

- $k=1$ hard bond

- $k < 1$ soft bond (kinda bouncy like a spring)

Let's try to implement this in pure functions

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 ( (*FB[*)((length)(*,*)/(*,*)(Norm[diff]+0.001))(*]FB*) - 1), maxDelta] diff
    },
      pts[[edge[[1]]]] += delta/2.0;
      pts[[edge[[2]]]] -= delta/2.0;
    ]
  ] ], RandomSample[pairs]//Transpose];

  pts
]

Here is a helper function to draw bonds, keeping them coupled to a dynamic symbol so that we can see how they change in real time.

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];

Show it using our randomly generated points

showGeometry[pts0, EdgeList @ bonds]

We are almost set for our animation! Let's regenerate random points and connections each time we evaluate the next cell and run the main update loop, which consists of Verlet methods combined with hard constraints

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]}] (*SpB[*)Power[\[Delta]t(*|*),(*|*)2](*]SpB*)
  ], {-5,5}];
  
  pts2 = pts1;
  pts1 = pts0;
  
  Pause[0.01];
  
], {i, 200}]

Animation widget

To avoid extra steps in evaluation, we can make a widget out of that with scoped variables

Module[{pts0, pts1, pts2, trigger, button, bonds, lengths, pairs, cache},
 trigger = EventObject[];
 button = InputButton["Restart"];

 EventHandler[button, Function[Null,
    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;
 ]];

 button // EventFire;

 {
  button,
  "","",
  Animate[showGeometry[pts0, EdgeList @ bonds], {n,1,100,1}, 

    "UpdateFunction"->Function[n,
      
      With[{\[Delta]t = 0.05},
        pts0 = Clip[applyBond[pairs, 1.0][
          2 pts1 - pts2 + Table[{0,-1}, {i, Length[pts0]}] (*SpB[*)Power[\[Delta]t(*|*),(*|*)2](*]SpB*)
        ], {-5,5}];
  
        pts2 = pts1;
        pts1 = pts0;
      ];
      (* prevent default action of Animate*)
      False
    ], 
  
    "TriggerEvent"->EventClone[button], AnimationRate->30
  ]
 } // Column // Panel 
]

Procedurally generated structures

Even without any optimizations of present code (normally you would need to use vectorization tricks or `Compile` or switch to `OpenCLLink`) we can still try to model some more complex than random graph of 10 points and add interaction with a user's mouse.

Rim

Here is a quick way of making this using just tables

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;

Prepare variables for storing positions

rim0 = 2.0 rim; rim1 = rim0; rim2 = rim0;
rimR = rim0;

mousePos   = {100.,0.};

Here we will automate our update cycle by attaching an animation frame listener to the `rimR` symbol. It will trigger an event to prepare the next frame of animation within the update cycle of your screen (actually - the app window). Each time, the timer is reset by an update of the `rimR` symbol, ensuring that the next frame will only come after the previous one has been finished.

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;
]}]

Animation controls

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}] (*SpB[*)Power[\[Delta]t(*|*),(*|*)2](*]SpB*), {-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;]]

Generate a Mesh from an Image

The more challenging, the more fun! It would be great to sketch a bitmap image using a brush and transform it into a mesh that can be modeled.

mask = NotebookStore["reburying-4d2"];
EventHandler[InputRaster[mask], (mask = #)&]

Now convert this image to a mesh, triangulate it and render as an index graph. Using `Graph` object is a bit easier, than working with triangulated mesh, since all bonds are sorted and can easily be extracted

mesh = TriangulateMesh[
  ImageMesh[ImageResize[mask , 100]// ColorNegate]
, MaxCellMeasure->{"Area"->100}];

graph = MeshConnectivityGraph[mesh, 0] // IndexGraph
vertices = graph // GraphEmbedding;
edges = (List @@ #) &/@ (graph // EdgeList) (*BB[*)(*convert to list of lists*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*);

vertices = Map[Function[x, x - {0.4,-1.2}], 0.04 vertices] (*BB[*)(*scale and translate*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*);

Animation cycle (*as before*)

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}] (*SpB[*)Power[\[Delta]t(*|*),(*|*)2](*]SpB*), {-5,5}]
  ];

  
  vertices2 = vertices1;
  vertices1 = vertices0;
  
 ], {i, 3}];
  verticesR = vertices0;

 ]&]; 
 verticesR = vertices0;
]

Button["Stop", EventRemove["anim"]]

Button["Reset", With[{}, vertices0 =  vertices; vertices1 = vertices0; vertices2 = vertices0;
verticesR = vertices0;]]