Extract individual contact forces and positions from Moco solution

[Nicholas Bianco]

Hi Pasha,

I think then this also defines all the ground contact deformation to take place in the contact sphere, with an undeformable contact halfspace.

This is not quite true. The model still accounts for deformation from the half space and assumes that the half space and contact sphere have the same material properties. Here is where the contact points are calculated in Simbody by taking into account half the indentation distance (i.e., since the material properties are assumed equal).

However, when the contact forces are applied to the half space body, they are transformed into body forces via Simbody's applyForceToBodyPoint. This computes the set of spatial forces and torques expressed in ground equivalent to the applied contact forces. The torques are computed by taking the cross product of the contact point in the half body (re-expressed in ground) with the applied contact force vector (Simbody implementation here).

All that to say: when you compute the spatial body forces on the half space from the contact sphere, the body torques will be about the origin, already accounting for true contact point location due to the Hertz force deformation. That's maybe more detail than you wanted, but sometime helpful to know where the devil is in these types of calculations.

As a check that I'm fully understanding this, essentially you mean Tz = Mz - x * Fy + y * Fx [using Z is up] or Ty = My + x *Fz - z *Fx [using Y is up]. All this assumes we're talking about the halfspaceforces and torques, all expressed in the ground frame with respect to the ground origin.

That is correct. I think this is what most users would expect from this utility: a set of external loads that mimic the loads usually applied to OpenSim models (e.g., when running the inverse dyanmics tool).

Best,
Nick

[Pasha van Bijlert]

Hi Nick,

Thanks for point me towards where everything is calculated, I think diving into the source code is quite useful for my own understanding of what's happening under the hood.

The model still accounts for deformation from the half space and assumes that the half space and contact sphere have the same material properties.

I was imprecise in what I'd written down, but I meant that we set COP_y = 0, so we're projecting it onto the ground without deformation.

For all my future simulations, I think it's most straightforward to simply compute the HalfSpaceTorques and HalfSpaceForces, as outlined in your example code, and then use this to compute COP. I still have the issue that I have a large batch of simulations where I didn't compute it in that manner, and ideally I would ensure I do everything in the same way. I did save all the contact sphere centroids, individual force components, and even the sphere torques as outputs from those simulations. But - I've only saved the SphereTorques, not the HalfSpaceTorques. The X and Y components are negatives of each other, but if I plot the Z components, the HalfSpace & Sphere torques aren't mirror images of each other (see below) I thought that force.getRecordValues(state).get(5) = - force.getRecordValues(state).get(11), but this appears not to be true. What does force.getRecordValues(state).get(5) represent?

sphere_and_halfspace_torques.jpeg (120.1 KiB) Viewed 813 times

As an alternative, I tried reverse engineering the contact point location (as per your links), so that I could compute the GRMs myself using the cross product between the contact points & the contact forces. I went through the steps (if it's helpful I can attach the code). I was having trouble understanding how the actual normal direction (with respect to the contact sphere origin) is computed. I think the relevant bit in the source code is "normalContactHalfSpace =
(bodyHalfSpace.getBodyRotation(state)*contactHalfSpaceFrame.x())", which I think in global (y is up) coordinates means the -y direction, unless you're not simulating flat walking. What I tried instead was subtracting the radius of the sphere from the y-position, and dividing this by 2 if it's negative (to represent the distributed contact deformation), but this again did not give me an exact match when I used it for the cross-product, and used the moments acquired that way.

I'm pretty much almost there, I either need to transform SphereTorque_z to HalfSpaceTorque_z, or more accurately reverse engineer how the contact point is being computed.

Any suggestions are welcome! As a last resort, I could reload all the .sto files and the respective models and compute the HalfSpaceTorques directly, but I'm hesitant to do this because it's error-prone, since I sometimes add programmatic modifications to the model before simulating it, and if I overlook one then the simulation won't be dynamically consistent.

Cheers,
Pasha

[Nicholas Bianco]

Hi Pasha,

While the the half space and contact sphere forces will be equal and opposite, the torques will not be since that are computed about different points (i.e., each frame's origin), even though they're both expressed in ground.

As a last resort, I could reload all the .sto files and the respective models and compute the HalfSpaceTorques directly

While possibly error prone, I think this is probably the most reliable approach for recomputing COPs. Re-engineering the applied spatial forces based on estimating the contact points, indentation, etc is likely to be similarly error prone (and similarly dependent on the model and solution files you use to compute them).

That's probably not quite the answer you're hoping for, but I do think it's the best approach.

Best,
Nick

[Pasha van Bijlert]

Hi Nick<

Fair enough, that's what I ended up doing. Thanks for working through this with me!

Cheers,
Pasha