NewtonUserJointAddAngularRow : what does it do ?

[Julio Jerez]

since the bodies moves on an inertial frame and the value of J is fixed in the frame, the J derivative is zero.

I see, I mean the pin is fixed to the frame of the body.
there is the physical interpretation, but it also comes from the geometrical and calculus derivation.
If you write the angular constraint equation, this is.

J . w = b

. Dot product operator
J is the pin
W the angular velocity
b some desired angular velocity function.

Remember we say we solve on the force and torque domain
So, we need to get the angular acceleration from the constraint equation. And that is achieved by taking the time derivative. So, you get

D(j.w) = D(b)

D is the time derivative operator.

We get

D(j) . W + j . D (w) = c

Let us go again over the easy terms.
C is the acceleration that the joint specifies ex (motor, free, fix, etc)

J. D(w) = J . a

That is the angular velocity of the body projected over the pin j.

The last term is where I confused you, It is very simple.
We just need to expand the expression and see what it leads.
To calculate that term, we must know what the value of D(j) is.

And that's where the meaning of the frame of reference is important.
We assume that j is of fix length, so it is invariance in any frame.
But j can change orientation with time, therefore it's time derivative is the
cross product of the pin and the pins angular velocity.
Let us call the pin angular velocity Wp
Now we can write.

D(j). W = (Wp x J) . W

Now we need to know how to get Wp.
And there are at least three common classic cases.

1- inertial frame, one that is not rotating, and not acceleratind
On this frame, the pin angular velocity is the body angular velocity,
simply because the pin is fix to the body.

Wp = w

So, we get

Wp. X J . W = W x J . W

But by the triple product identity

a x b . c= c x a . b

We then get

W x J . W = W x W . J = 0

Because a x a = 0
x is the cross product operator.

2-case. moving on a rotation frame.
Here the pin omega is the angular velocity of the body plus the angular velocity of the frame
wp = w + wf

Wf is the instantaneous angular velocity of the frame and may not be zero
plugging that into the equation and the frame inflict an acceleration of the body.

3-case. the pin is fixed to the frame of reference, example for this are contact joints.
Here D(j) is zero since is constant.

There can be other scenarios like pin that change in length, and things like that but, I stop covering those.

[Julio Jerez]

I did not write the pascal wrapper, a community use did, but I think he move on.
on this

Haven't you tested http://www.fmxexpress.com/create-wrapper-interfaces-for-c-and-c-libraries-using-swig-with-delphi-support/?

that's awesome, I will update swig this weekend and add the start wrapper for Delphi.
after the core base is added it is just a patter for adding the header of what you need.

technically, I should be supported already, the site says Delphi support was added on version 3.xx and I am using version 4.00 so I will try adding the script see what happens.

[blackbird_dream]

WOWWWW
AWESOME !!!!
Many thks.
I know 2 other people from this forum (at least) who will enjoy.

[JernejL]

If you can compile a dll with exposed C api i can make the pascal (and maybe other languages too) header translations easily.

[blackbird_dream]

if NDv4 can be parsed to free pascal, forget this reply.
I did a test with a body fixed to a double hinge(x,y)=angularRow along z and moving due to gravity. Its final position should be (2,2,-2srqt(2)) in theory. I observe a few percents error in position with my demo (1.97,2.06,-2.8 ).

Code: Select all

  NewtonBodySetLinearDamping(body,0);
  v:=vectormake(0,0,0);
  NewtonBodySetAngularDamping(body,@v.V[0]);     

[Julio Jerez]

yes, now you are seeing the difference between the theory and the implementation.

the first point is that the linear system the result from the derivation isn't well behave
they then to have singular values and they are also high condition number. you can check on those terms in wikipedia.
For the implementation, it means that the solution will have a large number or pivoting iterations.
and pivoting is a linear algebra operation that is uses to combine tow row or tow Colum.
the expression A = B + s * C
where the upper case means vectors and lower-case scalar.
the problem with the pivoting is that is it a high numerical error operation, so the more it is use
so it is in the best interest of the app to use algorism the reduce the among or pivot operation

the secund issues is that the integration is numerical, right there the solution will be an approximation by design.

so what a solver ha to do is to find a way to minimize the errors, while at the same time trying to be reasonably efficient so that can be uses for Realtime simulation.

Let us take the secund part. the numerical integration.
In the engine, a body goes to rest when the velocity and net force is zero.
that the theory, but the numeral negation in very few cases can achieve that goal.
so a practical heuristic approach is to be use instead.

These are different form version to version, but in general in 3.xx the heuristic is this.

there is a table look up table with few entries base of acceleration, and velocities base of ticks.
entries zero has the smallest value, and lower entices has large values until some value
entry zer0 could be
[1 0.001, 0.01] [8 0.01, 0.1] [16 0.1, 0.2] .....

that is if after the integration the body acceleration and velocity is below the acceleration and velocity on entry zero, the body set the sleep flag.
is it does no pass that test, the index moves to the secund empty and if is the acceleration is below those to values, set the sleep flag, and so on until the last entry.

next, it makes the disjoint set of all the connect bodies, static bodies are the break point since there are. these set are called islands.

them for each set in the union of all sleep flags is true, then the entire island is at rest and remove from the simulation until some evert occurs.

this is not perfect, but it is the prize for a practical implementation.

the application has the option to disable this behavior and by setting AutoSleepFlag,
this is a flag in the body that wit set to false, prevents the sleep flag for even become true

so, your first check is to set AutoSleep off on the bodies, because it seems the integrator is reducing the acceleration and velocity, and at some point, it is falling in one of the slots on the Sleep table and decaling the island of one body to be sleeping.


the island algorism for Newton 4 is very different that in 3.xx

in 3.xx an island is the set of all connected bodies boy joints the do no connect to a static body.
as you can see this make for a very spiky behavior, for example you cna have a large pile of object are rest and just touching one body, will cause the entire island to be simulated for not reason gain.

in newton 4 the island is defined as the set of all sleeping bodies connected by a joint, that are surrounded by a fence of sleeping,
a body is can only go to sleep when he is sleeping, the neighbor and the neighbor's neighbor are also sleeping. that's just a sophisticate graph algorithm that paralyzed very well.
and allows for removing the look up table.
so in newton 4 the look up table only has two entries, entry zero when velocity and acceleration fall below so value, in any frame, and secund entry is velocity fall below some larger value for about a secund worth of frames.

Ok try setting auto sleep off and see how that works.

[Julio Jerez]

for what I can see, you say you are not applying any linear constraint, so the body is just on a free fall.

Remember the integrator in newton for force is sympletic euler, you can look at up in wikipedia.
https://en.wikipedia.org/wiki/Semi-impl ... ler_method

this method simply changes the order of the tow equation and for a free fall reduces to this
1- calculate the velocity using the derivative at time t,
2- the use that velocity at time t + dt as if it was the velocity at time t
this is why this method is call semi-implicit,

in code that will be:
v1 = v0 + G * dt
x1 = x0 + v * dt.

yes, it seems it introduces a larger error, but believe or not the is far more stable that RK2,
if you run the equation above, you will see it generate the same result.

and RK2 for a free-falling body. will be

v1/2 = v0 + G * dt / 2
x1 = x0 + v1/2 * dt.
v1 = v1/2 + G * dt / 2

this will produce a perfect free fall integration,
there are also other methods, like predictor corrector, which I think is still in 3.xx and was the method all previous version of newton use, But I are not using any more, because predictor corrector methods are not very good for spiky secund derivative.

after all my attends I settled for semi-implicit.
and the reason is that implicit method tends to extract energy for the system, that is
the integration of a function is just the area under the curve.
explicit methods: RK1, RK2, RK4, predictor corrector, Euler, Verlet, etc all use the derivation of the function at time T, and form a trapezoid, then the integration is the area below the line the derivative
The image below e shows this.

Untitled.png (2.61 KiB) Viewed 12410 times

basically, is the area of the gray region.
For all the methods, what they do is that they bend the derivative by calculating few
more derivatives and averaging them,
but is all cases you can see that the gray area is larger than the actual area under the curve.

[Julio Jerez]

implicit and semi-implicit method on the other hand, calculate the derivative at time T + dt
then translate it to time T, and the integral is the area of the new trapezoid.
as show in the h figure below.

Screenshot_20220628_162731.png (4.28 KiB) Viewed 12408 times

as you can see, the area of the gray region is now smaller than the actual area.
in physics, this has the effect to introduce a natural drag, that in most cases is not necessary to add extra drag.

this is why I chose the the simple semi-implicit method, the stability that is offer far make up for the by design inaccuracy of the method.
here is the Good knew in your case.
the error is quadratic with the time step. It you set the sub sampling to say 4, them that error will be a lot smaller. so you have a way to minimized it, but it will never be zero.

[blackbird_dream]

for what I can see, you say you are not applying any linear constraint, so the body is just on a free fall.

No no:

with autosleep off it's the same. I understand the pb of converging to the steady state but it's not a big deal.