Truss Simulation - ProkopHapala/SimpleSimulationEngine GitHub Wiki

Linearized Truss Elasticity

• See class SoftBodyLinearized in SoftBody.h and
  • also function ( SoftBody::evalForceLinearizedBonds( ) and SoftBody::linearizedBonds( ) may be relevant
• Simple C++ implementation of LBFGS

Normally we calculate spring force due to bond between two points as

   d  = p[j] - p[i];
   l  = length(d);
   h  = h/l; 
   dl = l-l0;
   f = h*(dl*k);  

There are several problem with this approach:

1. numerical stability - the calculation of the length $l$ is numerically inaccurate because we subtract two large close number at two places.
   • d=p[j]-p[i]; is numerically inaccurate when the points {\vec p}_i and {\vec p}_j are close to each other but far from the origin ${\vec p}_0=(0,0,0)$. This can happen easily e.g. for space-ship simulation which has short sticks with lengh $~1[m]$ at the distant end of girders $~10000[m]$ from the center. Here we immediaterly loose 4-orders of numerical precission.
   • dl = l-l0 - is numerically inaccurate especially assuming that stick are rather stiff and usually cannot be strained more than $~1%$. Then this inaccurate value is magnified by multiplication by stiffness $k$ which is typically very large.
2. cost of sqrt() evaluation hidden in the calculation of length(d){ l2=dot(d,d); return sqrt(l2); };.

Both these problems can be minimized by using modified linearized algorithm where we first precalculate normalized direction vectors and than use deflections of points form the center.

void linearizedDynamics(int nsub=10){
    preCalc();
    // --- this low accuracy loop can be done e.g. on GPU in single-precission
    for(int i=0;i<nsub;i++){
       evalTrussForces_Linearized();
       movePoints_Linearized     ();
    }
    postCalc();
} 


double3 points[nPoint];
int2    bond2point[nBond]; 
float  l0s [nBond];  // lenght    of each bond

float  dl0s[nBond];
float3 hbs [nBond];  // direction of each bond
float3 dpos[nPoint]; // deflections of points from equlibrium

// --- high-precision pre-calculation

void preCalc(){ 
    for(int ib=0;ib<nBond;ib++){
        int2    b   = bond2point[ib];
        double3 d   = points[b.y]-points[b.x];
        double  l   = length(d);
        hbs[ib]     = (float3)(  d/l        );
        dl0[ib]     = (float )( l - l0s[ib] );
    }
    for(int i=0;i<nPoint;i++){ dpos[i]={0.f,0.f,0.f}; dvel[i]={0.f,0.f,0.f}; }
}

void postCalc(){
    for(int i=0;i<nPoint;i++){
        vel[i] += dvel[i];
        pos[i] += dpo [i]t;
    } 
}

// --- low-precision relaxation

void evalTrussForcesLinearized(){ 
    for(int ib=0;ib<Bond;ib++){
        int2   b  = bond2point[ib];
        float3 h  = hbs[ib];
        float3 dl = dl0s[ib] + dot( h, dpos[b.y]-dpos[b.x] );
        float3 f  = h*(dl*k);
        forces[b.x] += f;
        forces[b.y] -= f;
    }
}
void movePointsLinearized(){
    for(int i=0;i<nPoint;i++){
        // ToDo: here may be done another step of linearization ( assuming movement by constant velocity? )
        //       e.g. nerby point may move with fast velocity in the same direction. We may subtract this velocity to avoid large dpos[] ?
        dvel[i] += forces[i]*(dt/m);
        dpos[i] += dvel  [i]*dt;
    } 
}