[SOLVED] Position Based Elastic Rods implementation

[korzen303]

Good to hear that. I have made a fix in CMake files to turn the OpenMP off by default.

[mobeen]

Korzen, I think there is at least one place in your code where u r not testing for division by zero which might be the reason the solution diverges. Let me finish my implementation and then I will let u know.

[mobeen]

Hi Korzen,
As expected this was the reason. You are not checking for division by zero in line 88 to 92 in file PositionBasedElasticRod.cpp. These lines

Code: Select all

float p2pm_mag = p2pm.norm();
p2pm *= 1.0f / p2pm_mag;
lambda = (p2pm_mag - ghostEdgeRestLength) / wSum * edgeKs;

should be

Code: Select all

float p2pm_mag = p2pm.norm();
if (p2pm_mag > EPSILON)
	p2pm *= 1.0f / p2pm_mag;
lambda = (p2pm_mag - ghostEdgeRestLength) / wSum * edgeKs;

Doing so solves the issue and your code works for openmp as well.
Also you should also revert the permutation array to match those in the paper lines 29 to 35 in PositionBasedElasticRod.cpp file.

Code: Select all

const int permutation[3][3] = {
	0, 1, 2,
	1, 2, 0,
	2, 1, 0
	//0, 2, 1,
	//1, 0, 2,
	//2, 1, 0
};

Just a tip you could also replace the need for this array altogether by doing something like this in ComputeDarbouxVector function.

Code: Select all

for (int c = 0; c < 3; ++c)
{
//const int i = permutation[c][0];
//const int j = permutation[c][1];
//const int k = permutation[c][2];
const int i = c;
const int j = (c+1) % 3;
const int k = (c+2) % 3;
//...rest of code

Thanks for sharing.

[mobeen]

A quick question Korzen, in your code in function ComputeMaterialFrameDerivative, you call a function CrossProductMatrix (lines 306, 309, 312). I dont understand why you call it because in the given appendix in eq 43, 44, 45, there is no such thing. In eq. 43 and 44 there is a tensor product for d3 and d2. If I am not wrong, eq. 44 can only be solved if the second argument is a column vector. In that case, it should give me a 3x1 vector not a 3x3 matrix because the identity matrix subtracted from d2 tensor d2 should give me a 3x3 matrix. Then this matrix is multiplied to a 3x1 column vector which should give me a column vector? Why do you use the cross product matrix?

[bone]

Just a tip you could also replace the need for this array altogether by doing something like this in ComputeDarbouxVector function.

Code: Select all

<snip>
const int i = c;
const int j = (c+1) % 3;
const int k = (c+2) % 3;
</snip>

Are you trying to make it cleaner or faster? If the performance of that code is critical, the array is probably faster than taking the modulus, depending on compiler optimizations and platform. Faster still (but almost unreadable) would be adjusting your version with some sort of stupid rollover trick to simulate the modulus, like:

Code: Select all

// this code assumes 32-bit int's and that c is always in the range 0-2:
const int i = c;
const int j = static_cast< unsigned int >((c+1) * 1431655766) >> 30;
const int k = static_cast< unsigned int >((c+2) * 1431655766) >> 30;

[mobeen]

Hi bone,
Thanks for that I did not say anything about being faster but the permuted indices appeared straight forward mod with 3 thats why I wrote it but as you say, array indices should be fast so I trust you. Your constant based cryptic mod should be faster than my mod 3 version but that was just my 2 cents

[bone]

Ahh, no problem, I see your point.

[korzen303]

I did not have time to optimize the code. From a quick comparison it seems that nearly direct port to Eigen is 5-10 times (!) slower than the original implementation by Serphertoh. This is most likely because the lack of vectorization on Eigen's Vector3f and Matrix3f structures (I don't know why PBD library is not using SIMD optimized structures). Maybe changing layout of some matrices from col-major to row-major will also speed-up the calculations. I will get back to this next month but feel free to make any changes and I will pull them into the repo.

[mobeen]

Korzen, I have forked your repo here https://github.com/mmmovania/PositionBa ... ElasticRod
I have just removed all possible division by zeros from your fork.

[korzen303]

Thanks Mobeen, will merge it together with the OpenMP fix

[mobeen]

Hi Korzen,
I have finally got my own implementation of position based elastic rods. Attached is the snapshot from my demo.

snap.png

[mobeen]

I am now onto frame attachment constraint for orienting the tip of rod. I provide the frame attachment's initial frame such that the local z axis is pointing in the direction of world X axis. This way I would expect the rod's tip to be oriented in the direction of the world x axis. I get the desired result as shown below.

snap.png

This works fine but as soon as I move the tip around such that it orients in the opposite direction of initial frame, the rod starts to wiggle around unstably and does not converge. The behavior is like someone has added life to the rod and it wiggles like a snake hard to show in a snapshot but I will try below

wiggle.png

Anyone got any clues on what might be causing this. Is there a corner case I need to care about? Note that I store the initial frame which is always oriented such that the local X of initial frame is oriented to the world X axis.

[bone]

I had tried out Korzen's original demo and also experienced that wiggle exactly once. I couldn't figure out how to reproduce it so I never mentioned it. IIRC, the paper detailed a specific order of solving the constraints ... did you follow that?

[mobeen]

Hi Bone,
Thanks for your reply. No i did not follow the interleaving order for solving the constraints. let me try that and revert. Thanks for the pointers.

EDIT: Thanks bone that was the problem I added the interleaving order and the elastic rod is very stable no more wiggles

Now I will be on triangle attachments and then rigid body attachments to finish this up.
After all these are done, I will be adding this to my github repository soon for others to try.

[bone]

Glad I was of some minor help!

I look forward to checking out the implementation (someday, at least - it's not pertinent to my current project but I'd like to fiddle around with it in the near future).