I'm studying Jan Bender's position-based dynamics library (which is the closest to this paper impl i've found) and have lots of questions.

1) Why don't they have position constraints for contacts? why do they resolve interpenetrations only with velocities?

2) Why do they perform collision detection after velocity update, not after moving bodies to interpolated positions?

Their main loop is as follows:

TimeStepController::step:

1. Apply external forces, advance velocities, integrate positions (first rigid bodies, then particles - is the order important?)

2. Solve ("Project"?) position constraints (for X iterations: solve position constraints for typed constraints, but not for contact constraints (why?))

3. Update velocities (first rigid bodies, then particles)

4. Perform collision detection, create contacts (shouldn't this be done after the first step, after updating the positions?)

5. Solve ("Project"?) velocity constraints (for Y iterations: solve typed constraints, then rigid body contacts, then particle vs rigid body contacts)

3) Can you describe in a few words (or point me to beginner-friendly resources) how to solve position constraints for 3D contacts?

(or how to turn an abstract "Non-Linear Gauss-Seidel Solver" for solving position constraints into working code)

I'm not sure why that implementation would choose to solve collisions in velocity space. The trade off between the two approaches is that restitution > 0 is easier to implement consistently when you resolve collisions with impulses. However, you get much better stability for stacks and objects that are coming to rest if you solve the collisions positionally, and generally that is more important than having things bounce with perfect restitution. In systems I have worked on involving positional collision resolution for rigid bodies, I do a first pass to identify contact constraints, and for each of these I store the points on each body that collide, the normal of the contact, and for that frame the amount of work the constraint has done (depending on the body representation, I may have multiple constraints at different locations for the same body pair). Then when I iterate over the constraints, I remove whatever work that constraint contributed on the previous iteration and if the contact isn't satisfied then I solve it again. This allows collisions that don't contribute to the final solution to remove their influences without sorting by time of collision, and is essentially a positional form of the speculative contacts idea proposed here: http://www.wildbunny.co.uk/blog/2011/03/25/speculative-contacts-an-continuous-collision-engine-approach-part-1/

I can't answer your questions about the library but I can share my knowledge of the Position based dynamics paper.

While the particles that represent the rigid body dynamics have a velocity attribute, the gauss-sidel iterations work directly with positions.

For example a distance constraint is defined as C(p1,p2) = |p1 - p2| - rest_distance = 0. So you have to find a corrective value that will enforce the constraint. Like so C(P + delta_p) = 0. This is a short hand for saying find a value for value for p1 and p2 (p1 and 2 because the cardinality of the distance constraint is two particles), that when added to the predicted position (more on that in a minute) will move the particles into the correct place. The paper mostly talks about the framework used for all PBD, and it gives examples. To implement look at the final form of the constraint where you find delta p (the cardinality of the constraint determines how many deltas or corrective values you will need to solve that constraint. The distance constraint has two values for example).

I don't completely understand how the final constraints equations are derived but I have a rough idea why. In order to solve the constraint numerically using iterative methods (gauss-sidel and jacobi if you want to go parallel) first the Taylor seris is applied to the constraint equation to get its first order form (first derivative), then optimisation (the lambda) is performed to get a scaling value. The scaling value is important because it pushes the corrections into the right direction. Something about manifolds which I don't completely get.

So the bigger picture is this. Iterative methods work better if your initial value is closer to actual solution. So if you look at the algorithm in the paper at the beginning at the main loop you will see velocities and a predicted position calculated using forward euler (or was it verlet). The velocity allows external forces to influence your constraint e.g gravity applied at each step.

After correcting the position in the solver you update the velocities using the new corrected position.

You can use different collision detection methods, but you need to find the point of intersection, and then construct a collision constraint using that value, that will be solved like any other constraints.

I know this is not the best answer, but I hope it helps.

Turns out PBD is not very suitable for simulating absolutely rigid bodies (?).

I quickly had spheres jumping around on a plane, but couldn't get rid of jitter (becomes visible if increase gravity or decrease the timestep).

Maybe, that's why they solve contacts on velocity level (like in most impulse-based engines).

As for constraints equations, these lectures are extremely easy-to-understand:

[2014, 'so this is my *ss', lol, sorry]:

[2015]:

Rigid bodies are possible, I got it working. A few other steps can be taken to fix your problems.