search

jacktreado / dpm

Simulate deformable particles with OOP framework. Supports cell linked-list, header files partitioned into different particle types.
2 stars 4 forks source link

Negative eigenvalue in bending energy dynamical matrix #3

Closed jacktreado closed 3 years ago

jacktreado commented 3 years ago

negEval_config_issue negEval_spectrum_issue

While the above configuration is jammed, the newly minted function dpmBendingHessian2D creates one single negative eigenvalue (red square). When kb = 0 no negative eigenvalues were found, likely an issue in the implementation:

            // -- Stiffness Matrix

        // main diagonal block
        Hb(kx,kx) = Kb * (pow(dtiim1x,2.0) + pow(dtiix,2.0) + pow(dtiip1x,2.0));
        Hb(ky,ky) = Kb * (pow(dtiim1y,2.0) + pow(dtiiy,2.0) + pow(dtiip1y,2.0));
        Hb(kx,ky) = Kb * ((dtiim1x*dtiim1y) + (dtiix*dtiiy) + (dtiip1x*dtiip1y));
        Hb(ky,kx) = Hb(kx,ky);

        // 1off diagonal (enforce symmetry)
        Hb(kx,kxp1) = Kb * (dtiix*dtiip1x + dtip1ix*dtip1ip1x);
        Hb(ky,kyp1) = Kb * (dtiiy*dtiip1y + dtip1iy*dtip1ip1y);
        Hb(kx,kyp1) = Kb * (dtiix*dtiip1y + dtip1ix*dtip1ip1y);
        Hb(ky,kxp1) = Kb * (dtiiy*dtiip1x + dtip1iy*dtip1ip1x);

        Hb(kxp1,kx) = Hb(kx,kxp1);
        Hb(kyp1,ky) = Hb(ky,kyp1);
        Hb(kyp1,kx) = Hb(kx,kyp1);
        Hb(kxp1,ky) = Hb(ky,kxp1);

        // 2off diagonal (enforce symmetry)
        Hb(kx,kxp2) = Kb * dtip1ix * dtip1ip2x;
        Hb(ky,kyp2) = Kb * dtip1iy * dtip1ip2y;
        Hb(kx,kyp2) = Kb * dtip1ix * dtip1ip2y;
        Hb(ky,kxp2) = Kb * dtip1iy * dtip1ip2x;

        Hb(kxp2,kx) = Hb(kx,kxp2);
        Hb(kyp2,ky) = Hb(ky,kyp2);
        Hb(kyp2,kx) = Hb(kx,kyp2);
        Hb(kxp2,ky) = Hb(ky,kxp2);

        // -- Stress Matrix

        // main diagonal block
        Sb(kx,kx)       = (2.0*Kb*(dtip1 - dti)*ulxy[gi]/pow(l[gi],2.0)) + (2.0*Kb*(dti - dtim1)*ulxy[im1[gi]]/(pow(l[im1[gi]],2.0)));
        Sb(ky,ky)       = -Sb(kx,kx);
        Sb(kx,ky)       = Kb*(dtip1 - dti)*uld[gi]/pow(l[gi],2.0) + Kb*(dti - dtim1)*uld[im1[gi]]/pow(l[im1[gi]],2.0);
        Sb(ky,kx)       = Sb(kx,ky);

        // 1off diagonal
        Sb(kx,kxp1)     = 2.0*Kb*(dti - dtip1)*ulxy[gi]/pow(l[gi],2.0);
        Sb(ky,kyp1)     = -Sb(kx,kxp1);
        Sb(kx,kyp1)     = Kb*(dti - dtip1)*uld[gi]/pow(l[gi],2.0);
        Sb(ky,kxp1)     = Sb(kx,kyp1);

        // enforce symmetry in lower triangle
        Sb(kxp1,kx)     = Sb(kx,kxp1);
        Sb(kyp1,ky)     = Sb(ky,kyp1);
        Sb(kyp1,kx)     = Sb(kx,kyp1);
        Sb(kxp1,ky)     = Sb(ky,kxp1);

The correct version is implemented in MATLAB, check there first.

jacktreado commented 3 years ago

Fixed! line 2725 had an error, was dth[ci] = th - t0[gi];, changed to dth[gi] = th - t0[gi]; need to index with gi