Finite Elements Method: are there proble,s with stability?

[bone]

So you end up with this formula:

(J*W*JT + CFM) * lambda = -ERP * C / dt - J * v

Based on how CFM is defined in the links you provided, I'm pretty sure you need to divide it by the timestep in the above equation, so it should read:

(J*W*JT + CFM / dt) * lambda = -ERP * C / dt - J * v

[Dirk Gregorius]

No, it depends on how you derive it. Numerically they are identical as you can imagine since you can choose CFM independently. If you want to use the formulas from the ODE to relate CFM/ERP and ks/kd you need to devide by dt. (http://www.ode.org/ode-0.5-userguide.html#sec_3_8_2)

Here Erin derives everything using the same formulation as I did:
http://www.bulletphysics.org/Bullet/php ... f=4&t=1354
http://www.bulletphysics.org/Bullet/php ... +cfm#p9888



HTH,
-Dirk

[bone]

I said "based on how CFM is defined in the links provided," and I stand by that statement. You gave the formula for CFM based on Erin's definition, but gave links to ODE, which defines CFM differently. I am just trying to clarify it for people, so if they calculate their CFM based on real-life coefficients, the equations come out correctly.

EDIT: It's a reasonable point that CFM is just a number to some people, tweaked until it looks right. I guess I'm looking at this from a different point of view, in that I want it to be related to known physical quantities.

[Dirk Gregorius]

I agree and it is good that you pointed this out.


Anyway, CFM has nothing to do with spring modeling. It just happens that you can use it for this since you weaken the constraints. I lately found a good explanation of CFM in some papers of Claude Larcousiere. Maybe you find it useful. Basically a constrained mechanical system looks like this:

dx/dt = v
dv/dt = M^-1 * ( f_ext + f_c ) // and f_c = JT * lambda
C(x) = 0

You have three unknowns (dx/dt, dv/dt and lambda) and three equations, but you cannot relate them. This is why we solve on the velocity or acceleration level. Sadly this is not equivalent anymore with the original definition above since we *don't* consider the original constraint manifold anymore but its tangent space. This leads to the drift problems and is why we need stabilization.

It is a well know fact that constraint forces have (theoretically) infinite stiff potentials U = 0.5 * k * C(x)^2. If k = inf this can potentially only have a finite solution if C(x) = 0 (0*inf is undefined but can be finite). I now define k = 1 / eps and let eps -> 0.

f_c = -dU/dx = -1/eps * JT * C(x) = JT * lambda -> eps * lambda + C(x) = 0

Formally nothing has changes, but lambda can now be eliminated from the system. If e -> 0 nothing changes to the original formulation. I am no expert in this, but I think this goes into the realm of index reduction of DAE. If you are interested in this I recommend looking into the papers below. Claude is really the lord commander of this stuff:

http://www.cs.umu.se/kurser/5DV058/VT09 ... ure12a.pdf
http://www.cs.umu.se/kurser/5DV058/VT09 ... ure12b.pdf
http://www.cs.umu.se/kurser/5DV058/VT09 ... knotes.pdf

Starting from this formulation you can easily derive many steppers easily. (E.g. ODE, Box2D and Spook)


HTH,
-Dirk

[bone]

Funny you should mention Claude, as I just sent him an e-mail on his "Interactive Simulation of Elastic Deformable Materials" paper (which is older than the above) ...

Well I'll chew on the above papers for awhile. Working on anything interesting, Dirk?

[Dirk Gregorius]

I left Havok beginning of the year and work now as free consultant. Ironically my first job was again for Havok, but currently I am helping one of the major publishers to develop their inhouse engine. This is quite interesting. Besides this I nearly finished my engine which I use as reference implementation. I finished the convex-conves SAT with the Gauss-Map optimization and results and quality are great, but I need to implement some more optimizations to make it competetive to GJK/EPA. In any case I am now convinced that SAT is at least the better penetration solver than EPA. I also have a tendency looking into engineering again and spent some time looking into other formulations (e.g. Claude's work) and I do research direct sparse sub-solvers.


Cheers,
-Diirk

[fishboy82]

Hi Dirk:
what is " SAT with the Gauss-Map optimization ", could you give me a short explanation or some reference ? Thanks a lot

[Dirk Gregorius]

http://bulletphysics.org/Bullet/phpBB3/ ... 918&p=7437

[fishboy82]

Thanks Dirk, I will take a close look

[dust]

http://www.caelinux.com

there are a bunch of programs which are open source and might help about engineering and real world physics.

[bone]

M ( v2 - v1 ) = JT * dP
J * v2 + CFM * (P + dP) + ERP * C / dt = 0

J*M^-1*JT * dP = -ERP * C / dt - J * v1 - CFM * P

Sorry, I've been slowly working through these equations while doing other stuff. I think there's a missing term in the final equation above. It should be:

(J*M^-1*JT + CFM) * dP = -ERP * C / dt - J * v1 - CFM * P

If this is incorrect, I don't know where the CFM * dP term went from the equation above that one.

[Dirk Gregorius]

Yes, this is a typo. I will fix this. Thanks for pointing out.

[chevydevil-weimar]

Hello. I'm doing a very similar simulation. But in 3D and on the gpu and OpenCL. Forces are calculated for every point of tetrahedron, stored and finally summed up(Tedrahedron share Masspoints). Then I integrate the law of motion for every point. This works fine for the classic verlet approach. Now I want to try the velocity-verlet approach. But it seems to be very unstable and works only for very small time steps and/or spring stiffness. I figure I'm doing something wrong. Here is what I'm doing in pseudo code:

Code: Select all

force += gravitation;

position += velocity * t + acceleration * t * t *0.5f;
velocity += acceleration * t * 0.5f;

acceleration = force/mass;
acceleration += gravitation;

velocity += acceleration * t * 0.5f;

Am I overseeing something here?
Thx