Constraint forces and implementation

[Daniel Krueger]

The following two questions on Simbodys kinematic constraints were originally discussed between Sherm an me. However since others might be interested as well, I start this new topic:

1) In Simbody a constraint provides the method Constraint:getConstraintForcesAsVectors() to obtain the corresponding forces and moments on the constrained bodies expressed in ground frame. My question is: Is it possible to get the point of application of these forces or are they expressed as if they were acting on the bodys reference frame?

2) We found, that the SimTK::Constraint::Weld can be satisfied with multiple orientations of the constrained bodies. This ambiguity is also mentioned in the documentation. (e.g. at doxygen) I am thinking about implementing an umambiguous version of a weld constraint. However, to avoid running into problems you probably already considered while implementing the weld constraint I am interested in your opinion. Why has the current weld constraint been implemented this way?

[Daniel Krueger]

This is Sherm's reply:
1) The Constraint base class methods can only deal with features common to all constraints, but a "point of application" only makes sense for certain constraints. Some constraints act on joint axes rather than bodies (couplers), some apply pure torques to bodies (orientation constraints), and some may not have an obvious force location (say a cylinder constrained to roll on a plane). So the base class treats the constraint body forces as though they apply a pure moment to the body and a force at the body origin. For any specific Constraint type, however, like weld or point constraint, there might be a well-defined application point. In addition to inheriting the base class methods, each specific Constraint type has its own API that returns quantities that make sense for that type of constraint. So if you know the type of Constraint you are looking at, you can downcast the general Constraint to a specific one like Constraint::Weld and then have access to the Weld API.

It might be a nice base class feature to have some kind of "effective application point" that could be used for visualization, or could be the real application point in cases where that makes sense for the specific constraint type. Defining what that would mean precisely would have to be left to the specific constraint, though. If you think that would be useful, please file a GitHub issue describing what problem you would like it to solve and ideas for the API if you have thoughts on that.

2) The Weld's orientation constraints are implemented as they are for numerical stability during simulation. There are three dot product constraints enforcing perpendicularity between the axes of frames 1 and 2, with the constraints looking like x1.y2=0, y1.z2=0, z1.x2=0. As you noticed, these have the unfortunate property of being satisfied by several undesirable orientations! The obvious alternative that would be unambiguous (and probably what you have in mind) is the "parallel axes" constraint x1.x2=1, y1.y2=1, z1.z2=1. Those equations assemble beautifully but perform very poorly during a simulation, because they are extremely difficult to maintain. You can see why if you think of the dot product as computing a cosine. Cosine is very flat near one, and has maximum slope near zero. Thus when near one, small changes in angle have almost no effect on the constraint value. But near zero, small changes in angle have a large effect on the constraint. That means that numerically it is much easier to hold a dot product constraint near zero than to hold it near one. The runtime consequence is that the constrained angles would either drift from their desired values or an excessive amount of compute time will go into trying to maintain these poorly-conditioned equations.

So the end result is that the perpendicular-axes constraint performs poorly when unassembled but beautifully once assembled. The parallel-axes constraint performs beautifully when unassembled and poorly when assembled. Since a simulation must keep the constraint satisfied at all times, the perpendicular-axes constraint is the one to use.

My plan has long been to switch constraints during assembly, but I haven't gotten around to doing that. Any easy way to fix the problem in a specific model (if this is workable for you) would be to use three point-to-point constraints instead. That provides 9 constraints that are all translational and thus unambiguous in assembly. The 3 redundant constraints are what allow this set to function well in both assembly and simulation. Simbody will spread the forces over the 9 constraints; there is a small runtime cost to the additional constraints but it would probably be unnoticeable.

[Daniel Krueger]

Ok, and here what I did with it so far:
We are currently working on a project dealing with FEA simulations of human bones. Therefore we want to use force boundary conditions obtained from OpenSim models. This includes all muscle forces, joint reactions and constraint forces directly applied to a body. In case of the constraint forces it is necessary to know a physiologically reasonable point of application since in FEA the body (bone) is considered compliant. I think this can be achieved using the APIs offered by concrete classes derived from SimTK::Constraint. For example Constraint::Ball offers a method that will report its actual stations on the contrained bodies. Since in praxis I do not know the type of the constraints I need to check for some types before downcasting. Here is another problem. How can I check the type of constraint given only a generic SimTK::Constraint reference.
My first idea was using the isInstanceOf method, e.g.:
SimTK::Constraint::Ball::isInstanceOf(aGenericConstraint)-However this will never return true:-(
The second approach would be to actually do the downcast like:
SimTK::Constraint::Ball& aBall=SimTK::Constraint::Ball::downcast(aGenericConstraint)
- and then check if the returned reference is valid. Again, this seems always to be the case.
Any ideas?
Best regards,
Daniel

[Michael Sherman]

In the absence of a base class method that does what you want, isInstanceOf() and downcast() seem like reasonable choices. I'm glad you were able to find them -- they are very obscure!

I am surprised that didn't work for you -- maybe there is a usage problem, likely due to lack of good documentation. I built a model with 6 constraints, two of them ball constraints, and then wrote this code:

Code: Select all

    for (int i=0; i < matter.getNumConstraints(); ++i) {
        const Constraint& cons = matter.getConstraint(ConstraintIndex(i));
        bool isBall = Constraint::Ball::isInstanceOf(cons);
        printf("Constraint %d: %s a Ball\n", i, isBall ? "IS" : "is NOT");
        if (isBall) {
            const Constraint::Ball& ball = Constraint::Ball::downcast(cons);
            std::cout << "P1=" << ball.getDefaultPointOnBody1()
                << " P2=" << ball.getDefaultPointOnBody2() << std::endl;
        }
    }

Here is the output:

Code: Select all

Constraint 0: IS a Ball
P1=~[2,-1,0] P2=~[0,1,0]
Constraint 1: is NOT a Ball
Constraint 2: is NOT a Ball
Constraint 3: IS a Ball
P1=~[0,-1,0] P2=~[0,1,0]
Constraint 4: is NOT a Ball
Constraint 5: is NOT a Ball

Please check to make sure that you have your code in a try/catch block that reports error messages, and (if you aren't already) try running it in a Debug build where more extensive error checking will be done. My thought is that you may be getting an exception thrown but aren't seeing the message. That would occur if you call downcast() before checking isInstanceOf() for example.

Regards,
Sherm