Jakobsen particles+Verlet -> breakable joints?

Please don't post Bullet support questions here, use the above forums instead.
Post Reply
pinballwizard
Posts: 10
Joined: Tue Jun 19, 2007 5:52 am

Jakobsen particles+Verlet -> breakable joints?

Post by pinballwizard »

Hi all,

Recently I dusted off my 3D Jakobsen particle-based simulation code, implemented basically as Jakobsen described (embedded tetrahedrons in 3d collision geometry, barycentric coordinate/projection for collision resolution). I recently started implementing joints. To start off, I implemented a hinge joint by adding 2 new particles to represent the hinge axis, fixing these 2 new particles to the tetrahedron by adding 3 distance constraints per hinge particle (distance constraints between the hinge particle and a tetrahedron particle), and finally making the 2 hinge particles infinite mass so they would never change position (I also excluded the hinge particles from gravity/other forces). This works great and forces the tetrahedron to swing around the axis formed by the 2 hinge particles.

Now, I'm wondering about making breakable joints. To do this I would need to measure the force required to enforce a constraint and delete the constraint if the force exceeds a certain threshold.

Is it possible to measure the constraint force, or otherwise make breakable joints with the Jakobsen particle approach? Since the Jakobsen approach handles constraints by projection onto the constraint manifold, and not by calculation of constraint forces, it's not clear to me if I can somehow extract (or approximate) the constraint force with Jakobsen's particles and the Verlet integrator.

It seems there must be a way (the approach is so elegant, that the constraint force must be in there somewhere), but it eludes me at the moment. Any ideas?

Thanks.
Dirk Gregorius
Posts: 861
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA

Post by Dirk Gregorius »

I think Barenbrug derives forces from position constraints in his thesis. So you maybe want to have a look here:

http://home.iae.nl/users/starcat/dynamo
pinballwizard
Posts: 10
Joined: Tue Jun 19, 2007 5:52 am

Post by pinballwizard »

I think Barenbrug is explicitly solving for constraint forces, so he doesn't have a closed-form equation for the constraint force, but instead has to plug in the derivative/gradient/Jacobian/whatever-you-call-it :oops: of the constraint force into his quasi-Newton solver, and gets the constraint force as output.

With Jakobsen's particles there is no explicit expression for force or solving for force so I don't think Barenbrug's results are useful, as far as I can tell.

However I have another idea.

Jakobsen's suggestion for making a pin joint is to share a particle between two bodies. This shared particle gets tugged back and forth between the two bodies, and conversely tugs on the other tetrahedron particles of each body, until a constraint-satisfying configuration is found.

Instead of sharing the pin joint particle between the bodies, just duplicate it for each body and add another breakable zero-distance constraint between the two pin-joint particles. Then, each pin joint particle, unattached to the other, will be tugged along with its own body's tetrahedron during constraint solving, and at some point the zero-distance constraint between the two will be solved. At this point, check if the zero-distance constraint is exceeded by greater than some threshold, and if so, delete the constraint.

In theory this sounds like it will work (treating unconstrained joint position error as an indicator of constraint force), but I have to implement it to find out. I'm not sure how this will work with the constraint iteration loop - I can imagine that constraint ordering might be important e.g.:

1. Iterate over and solve all non-breakable constraints first (non-penetration constraints, rigid stick constraints between tetrahedron particles, and non-breakable joints).

After step 1, we have the configuration of the system after one timestep assuming all breakable joints are broken.

2. for all breakable joints, determine joint error
2a. if joint error < threshold, fix joint
2b. if joint error > threshold, break joint

3. Iterate over other constraints again (because step 2 might have violated other constraints again).

I wonder if it's possible to solve/break the breakable joints in the same loop as with all the other non-breakable constraints, or if you have to separate them in order as I outlined above. My fear is if you don't separate them, then the breakable joints might never break because the error is always kept small.

Thoughts?
raigan2
Posts: 197
Joined: Sat Aug 19, 2006 11:52 pm

Post by raigan2 »

I've been meaning to implement something similar but haven't gotten around to it yet -- is there something wrong with simply keeping a per-frame sum for each 0-distance-breakable-constraint?

i.e when calculating the error vector (during each solve of the breakable constraint), add the vector to a sum-of-vectors, and the vector's length (or squared length) to a sum-of-lengths; use these values (length of sum-of-vectors and absval of sum-of-lengths) at the end of the frame/after solving all constraints to measure stress/strain at that "joint". Possibly add some smoothing since this measurement is a bit noisy/random.

I haven't tried this yet of course, though something similar was recommended by someone doing impulse-based simulation.
pinballwizard
Posts: 10
Joined: Tue Jun 19, 2007 5:52 am

Post by pinballwizard »

Raigan, right, your suggestion also occurred to me last night and sounds reasonable. Basically every time you are "fixing up" a constraint, you are "effectively" applying a force/impulse, so by keeping a vector sum of these position correction vectors, you are effectively keeping a vector sum of the total "constraint correction force" or "constraint impulse" you applied in the current frame.

Also I think my previous suggestion to separate normal and breakable joints is not needed. However, this implies you need at least two iterations over all constraints (so that the effect of all constraints is felt by all other constraints). I mention this because John Schultz mentioned in another post that he had a one-iteration 100Hz stable-stacking Jakobsen particle simulator working. However I don't think he was worrying about joints, or breakable joints.

I'm thinking of coding up a simple test case of a double-pendulum with the weight initially located on the ceiling, then starting the simulation and letting gravity pull down the weight. The stress on the middle joint should be greatest when the weight is at its lowest point. Also by increasing the gravity you should be able to see the stress increase on the joint.

It's going to be a while before I get around to coding this, so feel free to post your results/thoughts here! :)
ewjordan
Posts: 26
Joined: Sat Jun 30, 2007 4:34 am

Post by ewjordan »

Having tried something like this before, I can tell you that this method works, kind of sort of, depending on what you want to use it for. Under a Verlet scheme, you're using the position and last position as stand-ins for the velocity, so there is a linear relation between the displacement of the current position and the impulse delivered; with a fixed time step, this is proportional to the force delivered to the particle.

The problem you'll run into is if you imagine an N-pendulum, with N particles all connected by linear distance constraints. Let's suppose for simplicity that it's hanging motionless at the stable lowest position. Physically the tension at each joint should increase the closer we get to the top, since the higher joints have to support all the weight of the lower particles. Using Verlet integration, however, it is easy to see that after one time step under the influence of gravity, each particle will have fallen the exact same amount. Hence the error in the pre-corrected values will be identical, and the error doesn't function as a reasonable proxy for joint tension.

The real problem is that forces are transmitted through chains of connected particles in this type of simulation over time - if you were to put a spring as the top connector, it would indeed get pulled more than if you put a spring as the bottom one, so the system of constraints does actually transmit the appropriate effective weights and forces. But even if you did use a spring instead of a rigid connection, you would only be able to read off the force after everything settled into equilibrium; with rigid constraints we cut that process out and jump directly to the equilibrium state, so I'm not sure that there's actually enough information coming into a single joint to get the force information out.

All that said, though, it does work for a lot of situations, and can look pretty good. Just bear in mind that you're not reading the true force/impulse, you're just reading the force/impulse contributions from the immediately connected particles. For a tetrahedron representation, this will usually mean that you're connected to all the particles in the two relevant bodies, which is pretty good, but don't expect to get everything absolutely correct this way.

Though if you have any ideas on this, I'd love to hear them - I kind of gave up fighting this one a while ago, so there may yet be a reasonable solution. Keep the N-pendulum in mind as a test case...
pinballwizard
Posts: 10
Joined: Tue Jun 19, 2007 5:52 am

Post by pinballwizard »

Very interesting post, ewjordan...
ewjordan wrote:Hence the error in the pre-corrected values will be identical, and the error doesn't function as a reasonable proxy for joint tension.

The real problem is that forces are transmitted through chains of connected particles in this type of simulation over time
Right. So, I wonder if it is possible to keep a connectivity graph of particles and somehow transmit the force (which, of course, would need to be computed/approximated somehow) by traversing the graph.
- if you were to put a spring as the top connector, it would indeed get pulled more than if you put a spring as the bottom one, so the system of constraints does actually transmit the appropriate effective weights and forces. But even if you did use a spring instead of a rigid connection, you would only be able to read off the force after everything settled into equilibrium;
Unless you used very stiff springs, which transforms this into a penalty problem with all the associated stiffness problems.

Still, hmm... a hidden penalty method with stiff springs for "approximating" the joint force (we don't care if the system blows up since its computations are not used for the final particle positions), and a standard Jakobsen projection scheme for the "real" movement of particles?
with rigid constraints we cut that process out and jump directly to the equilibrium state
Good observation. More generally, I would say that we cut the process out and jump to "a" constraint-satisfying state on the constraint manifold, but not necessarily "the" state which would have been found by a force computation method (e.g. Baraff).
so I'm not sure that there's actually enough information coming into a single joint to get the force information out.
Wow, this is a much more subtle and interesting problem than I thought. I still somehow think particle connectivity graphs and/or a hidden penalty scheme may be of benefit. Need to give this some more thought...
raigan2
Posts: 197
Joined: Sat Aug 19, 2006 11:52 pm

Post by raigan2 »

ewjordan wrote: Just bear in mind that you're not reading the true force/impulse, you're just reading the force/impulse contributions from the immediately connected particles.
This is only true if you iterate a single time; if you use several constraint iterations, the "forces" should propagate through the chain in a single frame.. the amount that they actually do propagate is dependent on the number of iterations and the order in which constraints are applied though.
pinballwizard
Posts: 10
Joined: Tue Jun 19, 2007 5:52 am

Post by pinballwizard »

raigan2 wrote:
This is only true if you iterate a single time; if you use several constraint iterations, the "forces" should propagate through the chain in a single frame.. the amount that they actually do propagate is dependent on the number of iterations and the order in which constraints are applied though.
OK, this is starting to get interesting.

I coded up a little 1D simulator (Jakobsen-style constraints on a line) in Python. I can post it here in its entirety if you want; it might serve as a useful baseline for discussion. In my simulator, I had a fixed infinite mass particle "X" at x=0, a particle "A" of mass 1.0 at x=1, and a particle "B" of mass 1.0 at x=2. Two constraints, a 1-distance constraint between X-A, and another 1-distance constraint between A-B. Gravity was set to +1. Verlet integrator with timestep 1.0 was used (artificially large to exaggerate gravity effects and force large constraint correction).

To compute the "constraint force", for each constraint between two particles p1 and p2, both p1 and p2 are moved to correct the error, scaled based on their masses. The movement of p1 is, we interpret, due to a "pulling correction" coming from p2; similarly, the movement of p2 is due to a "pulling correction" of p1. Everytime the constraint was corrected, I added the "pulling correction" to an accumulator for each particle. So we have constraint.totalp1correction and constraint.totalp2correction.

I ran the simulation with different mass values to see what would happen.

I can't interpret the results yet, but I can say two interesting things.

1. The total sum of the "pulling correction" values for a given constraint and particle (e.g. constraint1.totalp1correction) asymptotically approach some value dependent on the number of iterations needed to reach a constraint-satisfying state.
2. Large mass ratio differences can really slow down constraint correction. E.g. when I tried X with mass infinite, A with mass 0.01, B with mass 1.99, the system took very long to reach a stable state. This is because the X-to-A constraint yanks A back towards X without moving X (infinite mass). Then, the next constraint A-to-B yanks A back in the opposite direction towards B, moving B only a tiny amount (since B is much more massive than A).

So because of (2), the mass ratios of the particles can greatly affect how many iterations are needed to satisfy the constraints, and how much error correction is required at each step. This in turn affects the asymptotic values approached by (1), which are total sums of all error correction values.

I'm still not sure what the asymptotic values mean. Also there seems to be some interesting relationship between the total error correction at X and the total error correction at B (total error correction at X minus total error correction at B is 1.0). Also, this is a simple example - cyclic constraints, multiple constraints at a point, etc would be even more complex. So understanding this simple case is essential to understanding the more complicated cases.
ewjordan
Posts: 26
Joined: Sat Jun 30, 2007 4:34 am

Post by ewjordan »

Sorry, lost track of this thread for a while, just checked in today.

Upon review I have to take back a lot of what I said earlier about this not giving a good measure when interconnections are present. I opened up a debugging level for a Verlet-based game I've been working on and altered it so that the drawn width of a connection is proportional to the summed error over the time step. I figured it would be interesting to see what comes out of a realistic use of this idea (previously I'd used this quantity to trigger breaking, but I'd never actually checked to see what it looks like). Honestly, it seems to look just about exactly as it should, with the tension decreasing apparently linearly as you go further down the rope, so I think I was far too pessimistic before about this situation.

Having looked over my earlier post, I think my misgivings were based on a misunderstanding of the suggested method - I mistakenly thought that you meant to merely look at the total error per particle, not the cumulative absolute (or signed, perhaps, but still scalar) error per constraint. Taking the cumulative error is actually very reasonable, and has some solid physical motivation, I suspect, in the sense that in the limit of a perfectly consistent set of constraints iterated continuously until they are satisfied, we would probably see that quantity represent the correct tension (once factors of mass, etc. are properly handled). By continuous iteration I mean that instead of exactly applying each constraint at each iteration, we move towards it by an infinitesimal amount. Even if this is the case in theory, though (and I think it probably is, though I'm not about to go about proving it at the moment), we will lose something in practice because a) we can't do everything totally continuously, and b) a lot of configurations we specify will not be exactly consistent due to rounding errors, so we will introduce errors of about that magnitude during each iteration (this is where the signed sum of errors might help, because inconsistent configurations will tend to cause oscillations that would add to the unsigned error but average out to zero in a signed error sum). The issues with ordering of constraint evaluation are very relevant to a) - infinitesimal moves towards various constraint surfaces should commute, but finite ones may not.

In any case, you can take a look at http://www.gamefight.org/tensionTest/ and see if you have any thoughts - you'll have to accept an OpenGL security certificate first and maybe click the applet a few times to give it keyboard focus. Controls are r, g, b, left, and right, and you can grab the red rope if you're holding down r while the red ball touches it. I set the player location so you should have no trouble finding the rope. FWIW, the simulation is taking 5 steps per frame and running through the constraints twice per step, and the rope width doesn't go past 15 pixels because on my computer, at least, OpenGL starts to get mad at me and make the whole thing run really slow if I make it draw lines that are too fat, so you're not quite seeing the entire range of tensions, but it's close enough to get a feel for this.
pinballwizard
Posts: 10
Joined: Tue Jun 19, 2007 5:52 am

Post by pinballwizard »

Unfortunately I couldn't get the demo at http://www.gamefight.org/tensionTest/ to run on my Linux computer ("Error: Start failed: class not found: rgb2").

One question, though about your demo: did you try making the particles different masses? As I said, particles with different masses and large mass ratios can cause slow convergence, which affects the cumulative error values. After some experimentation I found out how to properly scale the cumulative error values to handle differing particle masses, but I can't recall right now exactly what I had to do.

In particular, the top-most link in the rope should feel a tension/force equal to the total mass*acceleration of all particles below it. This force on the top-most link should be independent of the mass distribution beneath it; it should only depend on the mass quantity beneath it. I.e. particles of mass 1+1+1+1+1 should yield the same top-link tension as do particles of mass 2.5 + 0.1 + 1.0 + 0.2 + 1.2.

I successfully got my 1D example to give these results, although a 1D example is somewhat simplistic. If you can get the same results in 2D or 3D, then it seems that indeed the cumulative error values can be used to reasonably estimate the constraint force.

Another example that would be interesting would be a stack of boxes (pushing forces instead of pulling forces in a rope). Assuming you are using embedded tetrahedrons with rigid stick constraints to represent your rigid bodies, then a stack of boxes should increase the tension/force felt by the rigidity constraints in the bottom box(es). Can you verify this?
pinballwizard
Posts: 10
Joined: Tue Jun 19, 2007 5:52 am

Post by pinballwizard »

pinballwizard wrote: Another example that would be interesting would be a stack of boxes (pushing forces instead of pulling forces in a rope).
Upon further reflection I realized that Jakobsen's paper suggests, due to the expense of collision detection, that collision detection and resolution only be run once per frame in the constraint-satisfying loop, while other constraints (joints, sticks) can be looped over multiple times per frame since they are cheap to compute.

This implies that in a box stack, the "pushing" force from a top box might not propagate downwards to the bottom box (because collision detection is only being run once per frame). Or, alternatively, the pushing force might take more than one frame to propagate. This might be OK if we sum up cumulative constraint correction over several frames to measure the constraint force.

Unfortunately my 3D Jakobsen/Verlet simulator isn't in a runnable state right at this moment so I can't run any tests right now, but this does seem like an interesting question (if box stacking results in appropriate pushing force on the bottom boxes).
ewjordan
Posts: 26
Joined: Sat Jun 30, 2007 4:34 am

Post by ewjordan »

Seeing as you can't load that applet (it's a general JOGL applet problem on some computers, I have no idea what causes it since I've never come across it myself, but it's not exclusively a Linux thing, it's been known to happen on OS X as well), I figured I'd run a test and post the screenshot:
Image

Ignore the three pronged thing, that's just the player, and the white things are just background crap, it's the two red strings we're interested in. Each of the particles on the left string has mass 1, and each on the right has mass 1 except the bottom one, which has mass 5 to make the total weight the same as the left one. As you can see, above the point where the weight differs, the tension (the fatness of each line, which is capped at 25 for unimportant reasons) seems identical in each rope.

The sum that I'm using to represent tension is the sum of each particle's displacement times its mass, and I'm computing it per linear constraint once a frame. This doesn't cover sign issues (a pull from a constraint should be offset by a push, so we need to be careful about this), but apart from that it seems to be what we want. Currently I have enough constraint iterations that the propagation happens in one frame, at some point I'll see what happens if I turn that down - for statics, it shouldn't make much of a difference if I'm correct, but I'll actually try it just to be sure.

I have some thoughts on the box stacking stuff (shouldn't need to consider multiple frames since equilibrium implies that the relevant quantity is identical in each frame), too, but let me try it first and see what's going on before I think too hard on it.
ewjordan
Posts: 26
Joined: Sat Jun 30, 2007 4:34 am

Post by ewjordan »

Okay, I looked into this a little bit more, and from what I can tell, this actually does seem to work at least to the level of physical plausibility. Here is my new linear constraint function, trimmed of a lot of irrelevant custom bookkeeping/optimization/early exit stuff: [this is Java code, BTW]

Code: Select all

//pA and pB are the particles we are applying the constraint to,
//which I store as members of the constraint class
//breakError sets the breaking threshold, either globally or per-constraint
//sumError is a per-constraint class member
//this.distance is the desired distance between the particles after application of the constraint

	public void applyConstraint() {
		float trueDistSqr = pA.position.distanceSqr(pB.position);
		float trueDist = (float)Math.sqrt(trueDistSqr);
		float alpha = this.distance / trueDist; //WARNING: need to check for div by 0 and handle somehow if you're actually using this code
		float oneminusalpha = 1.0f - alpha;
		float oneoversum = 1.0f / (pA.mass + pB.mass);
		float cmx = (pA.mass * pA.position.x + pB.mass * pB.position.x) * oneoversum;
		float cmy = (pA.mass * pA.position.y + pB.mass * pB.position.y) * oneoversum;
		pA.position.x = (cmx * oneminusalpha + alpha * pA.position.x);
		pA.position.y = (cmy * oneminusalpha + alpha * pA.position.y);
		pB.position.x = (cmx * oneminusalpha + alpha * pB.position.x);
		pB.position.y = (cmy * oneminusalpha + alpha * pB.position.y);

	   this.sumError += 2.0f*pA.mass*pB.mass*oneoversum*(trueDist-this.distance);
		if (this.sumError>breakError||this.sumError<-breakError) breakMe();
	}
The important bits for our purposes are the two lines at the end. Also, once per time step the sumError variable should be zeroed out - in my simulation I do this in the draw function, which I have at the top of the frame (actually, since I do several time steps per frame and I don't currently have a timeStepStart() or End() call to my constraints I am summing the error over several steps, though I'm not entirely sure if that's required - I still stand by my previous comment regarding the steady state assumption of statics leading me to believe this is not necessary), but YMMV. If you were to actually use this code and needed it to perform well you'd probably want to cache some of the computed factors that would be constant from call to call, but I left them written in for clarity here.

In both stacks and hanging chains, the per constraint error appears to appropriately change as you move up/down the chain/stack, so that in a high enough stack, the bottom boxes will collapse under the weight of the top ones, or in a long enough chain, the top pieces will snap off first. However, my simulation is not a pure Jakobsen one, so I'm not sure how exactly this stuff will carry over - in particular, I'm not attaching collision boxes to tetrahedrons, I'm instead treating the distance constraints as sticks and colliding directly against them (so a 2d box is made of four particles and six connecting sticks between them). BTW, I don't necessarily recommend this method, as it's not very efficient; that said, it does the trick for my purposes. But in a more standard tetrahedron representation, you may obtain some artifacts from the fact that the collision envelope does not match the particle-based geometry, I'm not sure. Some of it depends on the particulars of your body-body collision code, which Jakobsen did not explicitly cover (pure projection there can lead to some weird results). Honestly, you'll probably just have to try it and see if it suits your needs.

The other thing to keep in mind is that while these will currently break under enough compression or extension, a single stick will not snap under the normal force of a particle resting on it. For this you would need a different breaking model (I handle edge collisions more just like a standard 2d rigid line vs. point collision, so you could probably just read off the normal component of the impulse and check against a threshold, weighted by the observation that a stick is more likely to snap if hit in the middle than near the ends). But that stuff doesn't really carry over at all to the tetrahedron model, so I'll leave it up to you to decide whether it's worthwhile...

At some point when I get a chance I'll throw up a demo with a bunch of stuff that you can destroy (I amused myself for an embarrassingly long time by setting up a heavy sled that the player can ride and having it slide down a hill and decimate a pyramid of boxes, as well as making a rolling hamster ball of doom - a simple particle system effect when the things break is key here!), and I'll swap out the OpenGL renderer so you should be able to more easily view the applet.

One last thing - coming up with good breaking thresholds is pretty important, and I haven't yet figured out an automatic way to do it so things look realistic in general. Usually some things break too easily, and some don't break at all the way you'd expect them to, and it's all very dependent on the masses of the particles involved. I guess that's step 2, unless you're okay with doing it by hand every time...
Post Reply