A matrix and df_dx

[GoogleCodeExporter]

Hi, 

First of all, I would like to thank you for your open source. 

While I was looking at your code, I got some questions regarding A matrix and 
df_dx in implicit solver. 

In you code, number of non-zero elements in A matrix is always equal to the 
total number vertices(total_points). I believe it should have more non-zero 
elements. 

In terms of Jacobian matrix df_dx, let's say Kij is an element of Jacobian 
matrix and Kij = df_i / dx_j (i, j are vertex indices). The element of Jacobian 
matrix becomes non-zero when vertex i and j are linked by an edge(spring). In 
case i == j, Kii is always non-zero and Kii = -sum of Kij where j is linked to 
i and j != i. 

Implicit solver is supposed to handle very large stiffness such as 10^4 with 
natural gravity value.

Thank you.

Original issue reported on code.google.com by saggitas...@gmail.com on 21 Mar 2012 at 5:32

[GoogleCodeExporter]

In your code, df_dx is allocated as 'total_springs'. But only 'total_points' is 
used to compute A with df_dx. 

Original comment by saggitas...@gmail.com on 21 Mar 2012 at 5:42

[GoogleCodeExporter]

[deleted comment]

[GoogleCodeExporter]

Below is my code snippet in case it might be helpful..

    for ( int i = 0; i < (int)m_VertexArray.size(); i++ )
    {
        CVertexCloth3D& v_i = m_VertexArray[i];

        double h_over_mass  = h / v_i.m_Mass;
        double h_h_over_mass = h * h_over_mass; 

        I.SetElement(v_i.GetIndex(), v_i.GetIndex(), CMatrix33::IDENTITY);      
        dF_dX.SetElement(v_i.GetIndex(), v_i.GetIndex(), CMatrix33::ZERO);
        dF_dV.SetElement(v_i.GetIndex(), v_i.GetIndex(), CMatrix33::ZERO);

        // stretch forces
        for ( std::vector<int>::const_iterator iter = v_i.m_StrechSpringIndexes.begin(); iter != v_i.m_StrechSpringIndexes.end(); iter++ )
        {
            int index_j = m_StrechSpringArray[(*iter)].GetTheOtherVertexIndex(v_i.GetIndex());
            const CVertexCloth3D& v_j = m_VertexArray[index_j];

            double L0 = m_StrechSpringArray[(*iter)].GetRestLength();

            const CVector3D Xij = v_j.m_Pos - v_i.m_Pos;

            if ( Xij.Length() - L0 <= 0 )
                continue;

            CVector3D XijN = Xij;
            XijN.Normalize();
            const CVector3D Vij = v_j.m_Vel - v_i.m_Vel;

            //-------------
            // F0
            //-------------
            const CMatrix33 out_X = Xij.Out(Xij);
            const double inn_X = Xij.Dot(Xij);
            const CMatrix33 outOverInn = out_X / inn_X;

            // stretch spring force
            CVector3D fi = m_Kst * ( Xij.Length() - L0 ) * XijN;            
            CVector3D di = m_Kd * (Vij.Dot(Xij))*XijN;          
            CVector3D fsum = fi + di;

            F0[v_i.GetIndex()] += fsum * h_over_mass;
            f0[v_i.GetIndex()] += fsum;

            //-------------
            // df_dX
            //-------------
            CMatrix33 df_dx = m_Kst * (CMatrix33::IDENTITY - (L0/Xij.Length())*(CMatrix33::IDENTITY-outOverInn));

            df_dx = df_dx * h_h_over_mass;

            dF_dX.SetElement(v_i.GetIndex(), v_j.GetIndex(), df_dx);
            dF_dX.AddToElement(v_i.GetIndex(), v_i.GetIndex(), -df_dx);

            //-------------
            // df_dV
            //-------------
            //CMatrix33 df_dv = m_Kd * CMatrix33::IDENTITY;
            CMatrix33 df_dv = m_Kd * outOverInn;
            df_dv = df_dv * h_over_mass;

            dF_dV.SetElement(v_i.GetIndex(), v_j.GetIndex(), df_dv);
            dF_dV.AddToElement(v_i.GetIndex(), v_i.GetIndex(), -df_dv);
        }       
    }

Original comment by saggitas...@gmail.com on 21 Mar 2012 at 5:50

[GoogleCodeExporter]

OK. One more comment. Just to explain one line from my code.

if ( Xij.Length() - L0 <= 0 )
   continue;

Well, it is two lines :-)

In case the spring is compressed, if you ignore this case and don't put it into 
matrix, the linear solver becomes much stable. This is very simple but quite 
handy. Compressed spring can cause ill-conditioned matrix and CG solver can 
fail to converge. It is well explained in the paper below. 

http://graphics.snu.ac.kr/~kjchoi/publication/cloth.pdf

Thank you.

Original comment by saggitas...@gmail.com on 21 Mar 2012 at 5:57

[GoogleCodeExporter]

Hi saggita,
  Thanks for the feedback and comments. It will take me some while to revert to your questions since I am extremely busy these days. I will definitely reply as soon as I am free.

Thanks,
Mobeen

Original comment by mmmova...@gmail.com on 21 Mar 2012 at 1:34

[GoogleCodeExporter]

No hurry. 

Just to be clear. I was talking about the following solvers.

Implicit integration (Baraff & Witkin's model)
Implicit Euler integration

Original comment by saggitas...@gmail.com on 21 Mar 2012 at 4:30

[GoogleCodeExporter]

Hi saggita,
  Thanks for your input. I would surely read through these comments once I am free. It is always nice to have feedback from others. I hope that we could refined OpenCloth so that we have a solid base for cloth simulation and soft bodies simulation algorithms.

Take care and thanks for the input once again.
Mobeen

Original comment by mmmova...@gmail.com on 25 Mar 2012 at 6:31