Actually, I'd say the data structure comes out last from how
you need to use it.
On one hand we are
solving for the position and velocity (linear and angular) for
each object which can be done analytically, but on the other
hand the objects are interacting, which requires solving a
system of equations.
For example, with two
objects connected by a damped spring, you have a "phantom
object" of the spring, connected to the two object. I call
it a phantom object because we may not wish to model it as a
deformable object which is being collided with. For the
collision detector it's not really there, hence "phantom".
The damped spring doesn't constraint the motion of the objects,
as such, but provides a known force (it doesn't need to be
solved for as do constraint forces) based on the relative
distance between two fixed points on the objects (spring) and
how this distance is changing (damper).
This
requires the damped spring to query the system to determine the
location and velocities of the two points on those objects it's
connected to. From that it can add in the two forces on the two
objects for our calculation of the state of the system at time +
delta_t.
If we replace the spring with a
rigid rod (phantom or not), then we also need to determine the
acceleration of the points on the two objects so that they're
not accelerating away from each other to satisfy the constraint
of the rod's length. As the force we're applying with the rod
will affect the points' accelerations we need to calculate how a
unit force along the rod will change those accelerations. This
creates a system of equations we need to solve.
For
just two objects we might be able to hard-code solving a system
of just two equations, but the number of objects in the "contact
chain" may be greater than just two. I think this requires
inserting all the parameters into a matrix to solve the system
of n equations. If we do model the rod as a real object then
there will be at least three equations.
Since
a force applied at one constraint point may acclerate a body
such that the constraints at other points are violated, and
across several different contact points between bodies in the
contact chain, we have to iteratively converge on a solution for
all the forces. At least I think we have to do it iteratively.
With a force between objects A and B affecting the constraints
between objects B and C, and any force applied between B and C
then affecting the constraints between A and B, we have to
numerically converge on a solution rather than solve it
analytically.
We have the set of all objects
in the system, some of which may be "phantom" objects
and not analyzed by the collision detector. Each object has its
inherent physical properties and its state in global space
(other than the phantom objects). The forces in the system form
a set of force pairs between objects. Some of these forces may
directly calculable from the state of the system (springs,
gravity, jet exhaust, gyroscopic forces), while others need to
be solved for subject to constraints. The constrained forces
would be in the class of contact forces needed to prevent
interpenetration or joint forces needed to preserve the
constraints of the joints; I think both can be treated together
in a general way. These force pairs are either introduced by the
objects themselves (one object dynamically interacts with the
other objects in a way not covered by collisions or joints), the
collision detector, or through constant joints.
Oh,
and throw in the collision impulse calculations to prevent
inter-penetration of objects. Ah, I guess the same method by
which we calculate how a force/torque produces acceleration in
an object can be used for how an impulse will affect its
velocity. I think the equations would be identical.
In
general:
1. Start with a system with no
inter-penetrating objects.
2. Examine
collision points and their relative velocities, solving for
impulses necessary to prevent velocities which are "toward"
interpenetration. May be messy if this includes friction. Apply
those impulses.
3. Calculate and apply all
forces which do not depend on any other forces (unconstrained),
such as springs, dampers, gravity, motors, gyroscopic forces,
magnetism, etc.
4. Examine all the
constrained interaction between objects in terms of contact
points and force pairs. Solve numerically for a solution that
satisfies all constraints.
5. Advance the
simulation time. May use the simple Euler for
force/acceleration, assuming they're not changing over time, or
the modified Euler for a simple average of how forces are
changing over time, the Euler Predictor-Corrector to solve for
an average force over the time step more accurately, or a method
such as fourth-order Runge-Kutta for more initial accuracy. Halt
the simulation time before inter-prenetration has occured,
possibly needing to "back up". I'm not sure how this
step should interact with the previous steps. Anything other
than simple Euler evaluates the forces at (t + delta_t) to
estimate (dF/dt) and higher derivatives. Calculating the
unconstrained forces is not so time consuming, but re-solving
for constrained forces may be so. On the other hand, the more
costly calculations allow us to use a larger step size in the
simulation, so the simulation rate may be higher. On the
gripping hand, the previous steps may introduce discontinuities
into the forces over time (impulse-forces certainly do) which
aren't good for the ODE solving. However, impulse forces were
applied before we got here, and any further collisions
necessarily halt this step to re-apply impulses.
6.
Finish with a new system state with no inter-penetration at (t +
delta_t).
I'm working on exactly how to
formulate constraints and constraint solving in equations/code.
If a generic quadratic program-solver were used the geometry of
everything would all have to go in one big matrix. Since we know
conceptually the Big Matrix is formed out of several discrete
objects and forces, perhaps it would be fine to work on those
directly. We may still need a generic simultaneous
equation-solver as the number of force points change dynamically
as the simulation runs.
The calculations
involve querying the position, velocity, and acceleration of a
point on an object to determine if a constraint is satisfied,
and how the object's acceleration responds to force applied at a
particular point on the object to solve for forces to satisfy a
constraint. Once I've really worked out how the calculations are
performed then I'll know the information needed, which will
inform me of how that information should be stored. :)
The
problem is that there's really no simple intermediate case.
Collisions between frictionless rigid bodies is fairly trivial.
If you add static contact forces (which are a subset of general
constraint forces) the problem balloons-up to a quadratic
program where you have to simultaneously solve all the forces to
preserve non-interpenetration. Throw in fiction and it gets even
worse.
One example of some of the
difficulties: imagine a cube sitting on a table. To prevent
interpenetration with the table we must determine forces which
provide both zero velocity relative to the table at the four
vertexes of the cube in contact with the table. There is however
no unique solution. As long as the upward forces sum to the
force from gravity to give zero linear acceleration to the cube,
and the forces from diagonally-opposite points are equal to
provide no torque to the cube, the cube will not penetrate the
table. Each point in contact could have a quarter of the weight
of the cube on it, or two could have half-each and the other two
none. Depending on how we're attempting to solve for the forces
we could either divide by zero somewhere along the way, or
circle about never actually converging on a solution.
Anyway,
I'm working on doing the simple pendulum with more sophisticated
constraints. Because of the angular motion I'm afraid that
constraints on linear velocities and accelerations at the pivot
point will still allow the positional constraint to be violated.
I'm wondering how to fix that. If it were considered a collision
then an impulse would need to be applied just before the
positional constraint were violated. I'm not sure if I should
wait for the "slop" value in the position to be
violated to apply an impulse, or considered applying an impulse
whenever the relative velocity is non-zero. Even then, with high
accelerations the position may integrate to being outside the
allowed "slop" in one timestep.
Oh,
I did also play about with varying the time step based on the
forces/velocities at play. In the gyroscope simulation the
angular acceleration only prompted a change in time step when
the angular velocity was low. This makes sense, but also points
out that if you want to maintain a given accuracy for objects'
positions you may still have to scale the time step relative to
acceleration even though most of the time it's secondary to
velocity. |
