Skip to main content
gamedev.net
🔒 Locked

stiffness matrix

Started by Viper173 16 replies 2k views
Original Post
Viper173
Viper173
I finally found a paper that describes how to build a nice mass spring system !!!!! YES! But the very first point to do is way to much for me again. Of course! Could anyone please tell me a site where I can understand the idea of stiffness matrices and how I can precompute those? A Stiffness matrix, as I understand it, is a matrix that includes all necessary information for my entire mass spring system, right? That would explain all those 0's in stiffness matrices, because a single mass just has a few mass points that it's connected to (its neighbors). The other values represent the spring constant k for that mass point. .....??? Still, I dont have a clue how to use them!
Eelco
Eelco
i wont be of much help just yet, but could you please give me a linky to the paper you mentioned? id be very interested, since ive never found any even half decent explanation of implict solvers.

maybe after ive read through it i can help you out some.
Viper173
Viper173
there you go !!!

www.multires.caltech.edu/pubs/GI99.pdf

The method described in this paper is called an
implicit predictor/corrector scheme. So they say...

They claim stability between k's of 0 to 10^6 at 50 frames per second, with much less computanional efford than the actual implicit method.
Not bad, ehy? Do you see any disadvantages? I don't.

I just have extremly big problems to understand it!
If you are interested we might try to get it running. At least I would feel way more comfortable to be sure that some of my assumptions are right...

For a long time now I'm trying to implement a spring-system with high stiffness.
And pushing the update frequency to 300hz or so isn't quite the nicest way, right?
Viper173
Viper173
wait a minute, Eelco.

We had the pleasure about half a year ago.
At that time I started trying to find a solution.
The forum's moderator then suggested his paper which (I admitt it) I haven't understood quite yet up to now.

I was using an explicit method, which (now I know it) isn't usfull for small time steps and high stiffness.
Eelco
Eelco
hmm thats not exactly leisurly reading material... it has a nice explanation of the advantages of implicit methods, but after that it quickly gets more abstract.

i also cant say i completely follow them. i think a more easy-going explanation paper would be nice, as it wouldnt only fill in the gap of the stiffness matrix, but maybe all those other parts where they skip exactly the parts they consider trivial but id like explained aswell.

ill go and have a look. it something that ive wanted to understand fully for a long time.
Viper173
Viper173
well, me too.

Take a look at the provided pseudo code at the end:

First thing it says is: precompute W = ( I - [(dt*dt)/m] * H ) ^ (-1)

Well, I'm stuck there already :)

--> I stands for the identity matrix (that's for sure...I asked the author)
--> [(dt*dt)]/m is clear too. (turns F into a)

--> and now...the stiffness matrix!

let's say we have 4 mass points. we then get a 4 by 4 matrix.
Row 1 column 4 stands for the stiffness k between masses 4 and 1. ...and so on
But what about Row 3 column 3 ??? In the paper's example there's a -2 and a -1.
That's what I don't get...
Eelco
Eelco
yeah i didnt really get that either.

note though that this isnt the stiffness matrix, but its referred to with some other name. maybe it lends those properties from there, and my guess is, that for interpreting it as a stiffness matrix, it doesnt matter a thing what the values on the diagonal are.
Viper173
Viper173
Yes, H is referred to as the negated hessian matrix. Great!

For Calculating W however it does matter what the values on the diagonal are.

Think about that 4 mass point system again:
starting in the top left corner with mass 1 go clockwise to 2, 3 and 4.
(a rectangle)
Each mass is connected with the others with stiffness k except 2 and 3.

soooo.... my guess the matrix H would be :

-3 1 1 1
1 -2 0 1
1 0 -2 1 times k
1 1 1 -3


what's your guess?
Viper173
Viper173
anyway, the rows of my W matrix don't sum to 1. (as it's said to be the case in the paper)
Eelco
Eelco
first i looked at your matrix and wondered: where did he get that from. then i remembered something from the paper:

'meaning that each line sums to 1'

seems you are right, and you knew the answer all along (without knowing? :P)
Viper173
Viper173
geee, I don't know what's going on here...I'm right.

Yes, the matrix H seems to be correct. The negative numbers on the diagonal represent the number of connections that this mass point has.

AND, my matrix W comes out correct too. (rows sum up to ONE)
(I just used an identity matrix with ONE's everywhere instead just the diagonal...)

Well, that's W. I guess the next steps in the pseudo code aren't that hard.
The only two tough things to do are

1. Post correction of angular momentum
2. inverse dynamics
Viper173
Viper173
ok, that's it for me for now....it's about 1:18 am here and I had enough for today.


Have fun :)
Eelco
Eelco
Quote:
Original post by Viper173
geee, I don't know what's going on here...I'm right.

Yes, the matrix H seems to be correct. The negative numbers on the diagonal represent the number of connections that this mass point has.

AND, my matrix W comes out correct too. (rows sum up to ONE)
(I just used an identity matrix with ONE's everywhere instead just the diagonal...)

oh... entries of W need to sum to one?

anyway, an identity matrix as i know it only has 1's on the diagonal... well just wait and see what the result is.
Viper173
Viper173
Yes. First I got an incorrect result because I used a different I-Matrix.

But now I managed to get W for the described 4 mass point system, with k = 1000 N/m, dt = 0.01s, and mass = 1kg. And now the rows sum to 1. Then I guess it's right.

I'm moving on the computation of the internal forces now.
Viper173
Viper173
Do you have any idea what the big ^ sign in section 3 stands for?

You need that operator to get your global torque d(in old letter)T .
I have no idea. Is delta T a vector or a scalar???

in the code it says: delta T += F(filtered) ^ x(i)

I think "to the power of" doesn't make sense, because x(i) is a vector.
Eelco
Eelco
its quite an ambigious operator, but im rather sure its a cross product in this case.
grhodes_at_work
grhodes_at_work
I'll try to shed some light on the original question, though this might not really help use them.

If you're familiar with linear springs, then you know the formula below:

f = -kx

Here, F is the force generated by the spring, x is the displacement of the spring from its rest length, and k is often called the "spring constant."

In engineering terms, k is also called the spring stiffness, or just stiffness. It is a more formal term for the spring constant.

Suppose you have two springs, 1 and 2. Then, you can write their force equations as:

f1 = -k1*x1
f2 = -k2*x2

And, you can actually write this in matrix/vector form:

[f1]    [-k1   0 ][x1][f2] =  [ 0   -k2][x2]


Now, lets say that F is the vector [f1 f2], X is the displacement vector [x1 x2], and K is the 2 x 2 matrix above. Then, you can write the force equation as:

F = KX

Here, K is the stiffness matrix!

The thing to note is that the stiffness matrix is independent (in the linear spring model) of how far the springs are stretched, so it can be precomputed.

That's your introduction to the concept of a stiffness matrix and how they can be precomputed!

Note that the f's also are still equal to mass times acceleration. So you can write equations of motion for the spring endpoints as follows:

[m1  0][a1]    [-k1   0 ][x1][0  m2][a2] =  [ 0   -k2][x2]


Here, the 2x2 matrix on the left is called the mass matrix.

When the system is written that way, you can just integrate the rows in sequence using your integrator.

I will CAUTION you that I have deliberately oversimplified a bit. That example is fine for two springs that are not connected to each other. For connected springs, x1 and x2 are not independent, and so the equations as written as above (with the stiffness matrix being a diagonal matrix) can't be dealt with as is. The following PDF document has what I consider to be a reasonable example of a problem with coupled springs. In dealing with the fact that the x's are not independent, the stiffness matrix becomes a bit more complex, with the k's intermixed and the stiffness matrix no longer a diagonal matrix. This is the reality of coupled spring systems.

A Simple One-Dimensional Example of Discretization

I will also note that this whole idea is extensible to deal with structures more complex than mass-spring systems. For the more general case, the rows of the equation are not independent and you have to use an implicit or iterative solver at each time step.
Graham Rhodes Moderator, Math & Physics forum @ gamedev.net
Viper173
Viper173
Thanks, Mr Rhodes.

I'll take a look at the paper tomorrow.
I just wanna note somthing. Yes, the stiffness matrix doesn't stay constant.
But that's basically what the paper is about that we are currently trying to "translate". They used an approximation, wich as a matter of fact CAN be precomputed.

OK, Eelco. I take the ^ - sign as a cross product now. I will try it out tomorrow.
Gee, if I really get this algorithm working, then I will make an extra line in my life's TO-DO-LIST and a big cross right next to it. Right between "Sex with Twins" and "Walking on the moon"

I just can't believe, that these people made up all this stuff and even tested it in a virtual reality environment...

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.

Sign in to reply to this topic.