[Christiane]

I’m might be using totally the wrong English term, the slope is the first derivative of the spline. So if I want the local value of the slope at one point, I just calculate the vector between the two adjacent points. Because I used the Matlab spline function to create a series of closely set points, I can have a fairly localized estimation of the orientation. That worked well for me.

On the other hand, your current beam is set along the x-axis, so the orientation that you gave is correct.

Regards,
Christiane

[Christiane]

Hello again,

The beam points: I had 20 points that described my curved beams. Those are the points I used to mesh the beam, and for which I needed an orientation.

The Matlab spline function: I had the x, y and z coordinates of each point, so what I had was really 3 parametric splines, one for x, for y and for z, with t varying from 1 to 20. I simply input each spline into the spline function with t = 0:0.01:21. That gave me xx, yy and zz with a lot more points than my 20 initial points, and with extrapolation at both ends so I could determine the slope of those points as well.

With those points series, it was then very easy to calculate the 3D slope at my beam points. I then just transformed that into quaternions, and that is was I gave to SOFA as orientations. And yes, those orientations represent the x-axis of your beam.

Regards,
Christiane

[Christiane]

Hello Xiaojuan,

I used Matlab to define my spline and beam orientations. You can use the spline function with your beam points; it will help you define a continuous point series that represents your curved beam. Then you just use these points to determine the instantaneous slope at each point of your beam: that is the orientation you give to SOFA for your beam points.

Voilà!
Christiane

[Christiane]

Hello again,

I believe I found the problem, however I’m unable to correct it myself…

On line 288 of AffineMovementConstraint.inl, the rotation is calculated using quaternions. The final rotation is obtained by adding the initial orientation quaternion and the rotation quaternion. Normally, combining rotations is done by multiplying quaternions, not by adding them; in fact, a quaternion addition will not yield a unity quaternion, hence the result is not a rotation. There’s probably a check on the Rigid3Types which prevents keeping this non-unity quaternion as the orientation, and the end result is that I see my properly oriented beam segments go back to not being oriented until the end of the displacement, when they are reset to their original orientation.

Not being a progammer, this is beyond my meager c++ skills to correct.

Regards,
Christiane

[Christiane]

Hello Wong,

I had the same issue; I managed to stabilize it by setting the orientations of the DOFs in the mechanical object. Basically, you need to orient each node so that the x-axis is the slope of the curved beam.

I did this by creating a spline of my curved beams and then calculating the immediate slope at each node, aligning the x-axis with the spline. I didn’t bother with the y and z-axis, just oriented them randomly, because I was using beams with round sections anyway. I compared the resulting SOFA beams with an equivalent Ansys beam, and I had good agreement between the two.

Regards,
Christiane

[Christiane]

Hello Hugo,

Thanks, solving with the same solver for all rigids seems to solve the problem.

The next problem I’m having is with the JointSpringForceField. I’m able to specify the initial length using REST_T in its definition, which gets rid of any initial force in the springs; however, I’m unable to specify its initial orientation. I use REST_R to specify it, but the interface doesn’t seem to read it properly and it’s not applied. It just initializes its rotation by itself, and then when I solve I have an initial rotational load that changes my vertebrae orientations.

Any thoughts on that? Is it a bug, or do I have a typo someplace?

Regards,
Christiane

[Author]

Posts