yade-users team mailing list archive

[焦利青 <jiaoliqing@xxxxxxx>]

Dear All,

I tried to understand the partial codes in shop.cpp. I can't understand them. If you know the algorithms, could you tell me? Maybe you know where are the related algorithms from, could you tell me the information of the related papers. thank you so much.
The confusing codes are as follow:
void Shop::getStressForEachBody(vector<Shop::bodyState>& bodyStates){
    const shared_ptr<Scene>& scene=Omega::instance().getScene();
    bodyStates.resize(scene->bodies->size());
    FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
        if(!I->isReal()) continue;
        const NormShearPhys* phys = YADE_CAST<NormShearPhys*>(I->interactionPhysics.get());
        Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());    //For the moment only for Dem3DofGeom!!!
        if(!phys) continue;
        const body_id_t id1=I->getId1(), id2=I->getId2();
       
        Real minRad=(geom->refR1<=0?geom->refR2:(geom->refR2<=0?geom->refR1:min(geom->refR1,geom->refR2)));
        Real crossSection=Mathr::PI*pow(minRad,2);
       
        Vector3r normalStress=((1./crossSection)*geom->normal.dot(phys->normalForce))*geom->normal;
        Vector3r shearStress;
        for(int i=0; i<3; i++){
            int ix1=(i+1)%3,ix2=(i+2)%3;
            shearStress[i]=geom->normal[ix1]*phys->shearForce[ix1]+geom->normal[ix2]*phys->shearForce[ix2];
            shearStress[i]/=crossSection;
        }
       
        bodyStates[id1].normStress+=normalStress;
        bodyStates[id2].normStress+=normalStress;
       
        bodyStates[id1].shearStress+=shearStress;
        bodyStates[id2].shearStress+=shearStress;
    }
}
How to get the shearstress? I know maybe it is from thet=s-sigma(n), but I can't get the coherent expression. I think it should be "  shearStress[i]=phys->shearForce[i]-(phys->shearForce[i]*geom->normal[i]+geom->normal[ix1]*phys->shearForce[ix1]+geom->normal[ix2]*phys->shearForce[ix2])*geom->normal[i];". If you have the free time, could you help me to check it? thank you.

Regards
Liqing