Making systems of distance constraints come to rest

[asbjornlystrup]

I'm developing a 2D game with physics based on position-based dynamics and it's been working great, but I have some issues with "intersections" of distance constraints, and was wondering if there are any methods to mitigate this (apart from reducing the time step). This video shows the issue at hand, in the case of ropes: https://giant.gfycat.com/DizzyJaggedDikkops.mp4

In the video example, it eventually resolves itself, but often it stays intersected. The biggest problem is that it keeps on adding energy to the system, as I in general want the system to come to rest. I plan on using these distance constraints in a more general way to create some soft-body materials and cloth, so I have a feeling this could be even more problematic there.

I accept that the intersections/tunneling can happen. I'm more concerned with the system coming to rest. One way to mitigate it might be to check for intersections of distance constraints and temporarily increase their rest distance if they're not resolving, so it satifies both the collision constraints and the distance constraints at that point, but this sounds potentially slow and complicated?

This doesn't happen in the same way with ropes in 3D, but I guess an analogue there would be intersecting planes, so I'm curious if there exists any methods for mitigating the energy problems in this analogue.

This is done in a Jacobian manner on the GPU. Maybe Gauss-Seidel deals with this better?

Thanks!

[janbender]

Typically there should be no problems if distance constraints intersect. However, maybe you have additional constraints which might effect this?

I also don't think that the Jacobi solver should be a problem.

Can you give me more info about that problem? Are there other constraints involved? Do you have high mass ratios? How is the distance constraint implemented?

[asbjornlystrup]

All particles have a collision constraint; they can all collide with each other, which I think is the reason for the "knots" adding energy. All particles' mass is the same (apart from mass 0 on the static particles, but that can very likely be ignored), and there's no special "rope entity" code, there are just general particles and distance constraints. Here's an illustration: The green lines are (practically) satisfied distance constraints, and the red unsatisfied. I think the problem is that the distance constraints marked in red here are incapable of getting to rest, because to satisfy the 4 particles' collision constraints, the distance constraint can't be satisfied, so I guess there's oscillation there that adds energy to the system. Maybe there's some way to smooth out this oscillation? The rest distance on all distance constraints are the same, which is set to be slightly above the diameter of the particles.

At the beginning of each simulation step, I first sort all particles for neighbor search. Then I iterate this pipeline a few times:

• For each particle, all neighbors are iterated to sum up a delta required to separate potential circle intersections. This delta is divided by the number of intersections that affected it and then added to the predicted position.
• Then the distance constraints are solved in a constraint-centric manner. For each constraint, I atomically gather the delta required to satisfy the rest distance between the two particles in the constraint. I also atomically sum up the constraint count. When all constraints have been iterated, all particles are iterated and their new delta is divided by the distance constraint count and then added to the predicted position.

Everything is Jacobian in the way that the kernel reads from one buffer and writes to another. After each kernel dispatch, I swap the read and write buffers.

[asbjornlystrup]

I was reading more into "A Survey on Position Based Dynamics, 2017" and noticed that in the linear systems presented in 4.2.1, each different constraint type works on the same input position.

I assume in a Gauss-Seidel implementation, you update the predicted positions as immediately as possible? And the most "Jacobian way" would probably be to not update the predicted positions until all constraint types are processed, making the constraint types' execution order-independent.

The "Unified Particle Physics for Real-Time Applications" paper seems to have a concept of "constraint groups" wherein each group, the constraint types' execution is order-independent. I combined the collision and distance constraints into such a group, so that the predicted position is not updated until both these constraint types are processed. This seems to have fixed most instabilities, but it converges very slow.

I'm curious about the practical methodology for determining the grouping and order of constraint types within each solver iteration, and ways to solve instabilities related to this. I have some ideas for what to try next, such as the successive over-relaxation mentioned in the latter paper. And I also wonder if it could help to divide the position deltas by the constraint counts related to that delta only, as opposed to dividing the two constraint types' collective delta by their collective constraint count.

Any information is greatly appreciated!

[janbender]

I think you have to avoid interpenetrations to solve the problem. When we simulated elastic rods in 3D we used capsules to detect the collisions and not spheres. Here you can find information how to compute the distance between two capsules:

https://iquilezles.org/www/articles/distfunctions/distfunctions.htm

In the PositionBasedDynamics we also use constraint groups to parallelize the Gauss Seidel method. Maybe you could take a look how we implemente that. Note that in most cases the Gauss Seidel method converges faster than the Jacobi method.

[asbjornlystrup]

The capsule idea sounds interesting. Does this mean you have the ends of the capsules overlapping, and filter out the 2 closest capsules from collision per capsule? I have other types of materials in the game, so I might stick to particles to make unification easier, but implementing a collision solve filter to allow some overlap in the rest state of the ropes might be nice.

I couldn't find any relaxation code related to the constraint groups in the repo, but maybe I didn't look hard enough. Have you implemented any constraint relaxation schemes in this repo?

[janbender]

Yes, the capsules are overlapping and we don't test for collisions between two adjacent elements.

What exactly do you mean by "constraint relaxation schemes"?

[asbjornlystrup]

PhysX Flex has two relaxation modes for example, that they call global and local (source). Global means you simply multiply position deltas by a constant factor, and local is the more stable "constraint averaging with successive over-relaxation", as in the delta is divided by the number of constraints that affected it and multiplied by a constant factor (usually around 1 to 2). I think Flex has collision, density, and spring constraints in the same constraint group; as in they add them to the same delta (so local mode is very jacobian), but not completely sure.

Curious if there are other approaches to relaxation that fits into this system. E.g. if you could divide the delta of each of the different constraint types by only its own constraint count, but you still wait with updating the predicted position until all are solved or something. I have yet to try this. Just looking for ways to avoid instabilities similar to this "constraint averaging" method, but with faster convergence. Updating positions between each constraint type seems to be more unstable. Maybe implementing the second-order integration described in the PBD survey will make it converge faster. Are there any stability problems or other issues with second-order integration that you know of?