I am trying to implement the "Position Based Elastic Rods" paper:
http://www.nobuyuki-umetani.com/Positio ... icRod.html
All the constraints are nicely derived in the paper's appendix except the bending and twisting constraint, for which the authors do not provide the final position update.
I know that the derivatives of the Darboux vector given in the appendix need to be plugged into equation (17) in the paper, but I am not sure how to do this.
This is what I got so far using low constraint stiffness: https://www.youtube.com/watch?v=mQdwS2k4AZ4
I provide the full source code and Unity3D project below. I triple checked all the equations code as well as matrix ops should be fine (although I can't 100% guarantee).
Big thanks in advance for help
korzen303
EDIT: I have just reveived a hint: "Equation (17) is an generalization of the typical PBD projection into vector valued constraint. In the bending-twisting constraint, the constraint is vector value (C is a 3-vector). Hence \nabla C is a 3x(3*5) matrix {I think he meant (3x3)*5, korzen303}. Please compute 3x3 normal matrix (\nabla C)^T(nabla C). And compute inverse.
The appendix shows how to compute the nabla C:
http://www.nobuyuki-umetani.com/Positio ... pendix.pdf"
I have edited the following source code accordingly, but still getting unstable behaviour.
EDIT: Get the working code here https://github.com/korzen/PositionBased ... ElasticRod
The constraint method for a quick glance:
Code: Select all
/// <summary>
/// Projects bending And twisting constraint
/// </summary>
/// <param name="n">Constraint index</param>
/// <param name="pi">pA, pB and pC are positions of subsequent centreline points. pD and pE are ghost points</param>
/// <param name="wi">Inverse masses of the corresponding points</param>
/// <param name="dA">Material frame A consisting of pA, pB and pD. Must be recalculated at each constraint update</param>
/// <param name="dB">Material frame B consisting of pB, pC and pE. Must be recalculated at each constraint update</param>
/// <param name="darbouxVector">Darboux vector of this element. Must be recalculated at each constraint update</param>
void ProjectBendingAndTwistingConstraint(int n, ref Vector3 pA, ref Vector3 pB, ref Vector3 pC, ref Vector3 pD, ref Vector3 pE,
float wA, float wB, float wC, float wD, float wE, ref Matrix3x3 dA, ref Matrix3x3 dB, ref Vector3 darbouxVector)
{
ComputeMaterialFrame(n, pA, pB, pD, ref dA);
ComputeMaterialFrame(n + 1, pB, pC, pE, ref dB);
ComputeDarbouxVector(n, dA, dB, ref darbouxVector);
ComputeMaterialFrameDerivative(n, pA, pB, pD, dA,
ref m_dajpi[0, 0], ref m_dajpi[0, 1], ref m_dajpi[0, 2],
ref m_dajpi[1, 0], ref m_dajpi[1, 1], ref m_dajpi[1, 2],
ref m_dajpi[2, 0], ref m_dajpi[2, 1], ref m_dajpi[2, 2]);
ComputeMaterialFrameDerivative(n + 1, pB, pC, pE, dB,
ref m_dbjpi[0, 0], ref m_dbjpi[0, 1], ref m_dbjpi[0, 2],
ref m_dbjpi[1, 0], ref m_dbjpi[1, 1], ref m_dbjpi[1, 2],
ref m_dbjpi[2, 0], ref m_dbjpi[2, 1], ref m_dbjpi[2, 2]);
ComputeDarbouxGradient(n, darbouxVector, m_midEdgeLength, m_dajpi, m_dbjpi,
dA, dB, m_bendTwistKs,
ref m_deltaOmegas[0],
ref m_deltaOmegas[1],
ref m_deltaOmegas[2],
ref m_deltaOmegas[3],
ref m_deltaOmegas[4]);
Vector3 C = new Vector3( m_bendTwistKs[0] * (darbouxVector[0] - m_restDarbouxVector[0]),
m_bendTwistKs[1] * (darbouxVector[1] - m_restDarbouxVector[1]),
m_bendTwistKs[2] * (darbouxVector[2] - m_restDarbouxVector[2]));
Matrix3x3 factor_matrix = new Matrix3x3();
float[] invMasses = new float[5] { wA, wB, wC, wD, wE };
for (int i = 0; i < 5; i++)
{
//cannot invert the matrix if multiplied by w = 0 (i.e. locked point)
// factor_matrix += (m_deltaOmegas[i] * m_deltaOmegas[i].Transpose()) * invMasses[i];
factor_matrix += (m_deltaOmegas[i] * m_deltaOmegas[i].Transpose());
}
factor_matrix = factor_matrix.Inverse();
Vector3[] deltaPos = new Vector3[5];
for (int i = 0; i < 5; ++i)
{
deltaPos[i] = m_deltaOmegas[i] * (factor_matrix * C);
}
pA -= (m_bendTwistMult * deltaPos[0]);
pB -= (m_bendTwistMult * deltaPos[1]);
pC -= (m_bendTwistMult * deltaPos[2]);
pD -= (m_bendTwistMult * deltaPos[3]);
pE -= (m_bendTwistMult * deltaPos[4]);
}