[SOLVED] Position Based Elastic Rods implementation

Please don't post Bullet support questions here, use the above forums instead.
korzen303
Posts: 30
Joined: Thu Dec 19, 2013 12:13 pm

[SOLVED] Position Based Elastic Rods implementation

Post by korzen303 »

Hello All,

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]);

    }
Last edited by korzen303 on Fri Jan 08, 2016 2:35 pm, edited 1 time in total.
korzen303
Posts: 30
Joined: Thu Dec 19, 2013 12:13 pm

Re: Position Based Elastic Rods implementation

Post by korzen303 »

Thanks Mobeen,

I based the above constraint code on Serpheroth's implementation. Remotion's , whose demo from the post works very well, did not reveal the code.

I will keep digging
korzen303
Posts: 30
Joined: Thu Dec 19, 2013 12:13 pm

Re: Position Based Elastic Rods implementation

Post by korzen303 »

The zip files contain only the readMe file, similarly as the GitHub repo. Anyways, I have managed to solve the problem. Will update the source files soon.
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Position Based Elastic Rods implementation

Post by mobeen »

The zip files contain only the readMe file, similarly as the GitHub repo. Anyways, I have managed to solve the problem. Will update the source files soon.
.
Really sorry I had not seen the zip contents. I really don't understand why people create git repos but don't add sources? Whats the point of keeping an emplty open source repo. I m glad that u were able to solve ur problem
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Position Based Elastic Rods implementation

Post by mobeen »

https://bitbucket.org/yukikoyama/view-d ... on-for-3d/

have found an old repo from the original author. Hope this will help.
bone
Posts: 231
Joined: Tue Feb 20, 2007 4:56 pm

Re: Position Based Elastic Rods implementation

Post by bone »

I haven't looked, but it's possible that Jan Bender's PBD code has an implementation, too. But no matter since you've apparently found a solution.
korzen303
Posts: 30
Joined: Thu Dec 19, 2013 12:13 pm

Re: Position Based Elastic Rods implementation

Post by korzen303 »

https://bitbucket.org/yukikoyama/view-d ... on-for-3d/
have found an old repo from the original author. Hope this will help.
This is a different paper (except that the first author is also Japanese). But thanks Mobeen anyway, they also use some sort of elastic rod using PBD and oriented particles [Müller and Chentanez, 2011] which is interesting plus they provide the reference source code.
I haven't looked, but it's possible that Jan Bender's PBD code has an implementation, too. But no matter since you've apparently found a solution.
It is not there, but I am willing to contribute this to Jan's excellent library.
Last edited by korzen303 on Fri Jan 08, 2016 8:35 pm, edited 1 time in total.
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: Position Based Elastic Rods implementation

Post by mobeen »

Korzen,
Just for completeness since you have already found the problem, could you let us know what was wrong? or perhaps highlight how PB elastic rods are implemented for the rest of us? It would be a good idea placing your source code on github (if possible) as a reference for others who might fall in trouble when they try their implementation. (I might give it a try over the weekend :))

Thanks,
Mobeen
korzen303
Posts: 30
Joined: Thu Dec 19, 2013 12:13 pm

Re: Position Based Elastic Rods implementation

Post by korzen303 »

There was a couple of minor bugs. The main thing was to change the indices when computing the Darboux Vector (not entirely sure why). Again, credits to Serpheroth:
const int permutation[3][3] =
{
//0, 1, 2,
//1, 2, 0,
//2, 1, 0
0, 2, 1,
1, 0, 2,
2, 1, 0
};
korzen303
Posts: 30
Joined: Thu Dec 19, 2013 12:13 pm

Re: Position Based Elastic Rods implementation

Post by korzen303 »

I have made an implementation for PBD library. It is basically a rough and dirty port of Serpherot's code to Eigen. I will need to spend some more time on matrix ops optimizations. Very likely that his implementation is much faster.

The repo is here. I have asked Jan to pull it into PBD main branch.
https://github.com/korzen/PositionBased ... ElasticRod

And the binary for a quick look:
https://github.com/korzen/PositionBased ... odDemo.zip
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: [SOLVED] Position Based Elastic Rods implementation

Post by mobeen »

Hi korzen,
Sorry to revive this thread again but I think there are still some problems with the current implementation. The solution diverges in the release build and apart from the first two nodes, the rest of nodes disappear to infinity. The debug build works fine but is awefully slow. I am running it on an Intel Core i7 CPU @ 4GHz with 16GB RAM and an NVIDIA GeForce Titan Black GPU on Windows 8.1 64-bit machine. The demo binaries that you shared work fine though.
korzen303
Posts: 30
Joined: Thu Dec 19, 2013 12:13 pm

Re: [SOLVED] Position Based Elastic Rods implementation

Post by korzen303 »

Check the readMe file in Demos/ElasticRodDemo:

"IMPORTANT: The OpenMP is not supported in this release!
By default, the CMake generated ElasticRodDemo project has OpenMP turned on.
Deactivate it to run the demo, otherwise to rod will explode."

Sorry for the inconvenience
mobeen
Posts: 122
Joined: Thu May 05, 2011 11:47 am

Re: [SOLVED] Position Based Elastic Rods implementation

Post by mobeen »

Check the readMe file in Demos/ElasticRodDemo:

"IMPORTANT: The OpenMP is not supported in this release!
By default, the CMake generated ElasticRodDemo project has OpenMP turned on.
Deactivate it to run the demo, otherwise to rod will explode."

Sorry for the inconvenience
Yep sorry I didn't bother reading the readme :P works fine after turning openmp off :)
Post Reply