yade-users team mailing list archive

[behzad <question258258@xxxxxxxxxxxxxxxxxxxxx>]

Question #258258 on Yade changed:
https://answers.launchpad.net/yade/+question/258258

    Status: Answered => Open

behzad is still having a problem:
It calculates the contact forces:

bool computeForceTorqueCohBurgers(shared_ptr<IGeom>& _geom,
shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r &
torque1, Vector3r & torque2) {

  const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
  CohBurgersPhys& phys=*static_cast<CohBurgersPhys*>(_phys.get());
  Scene* scene=Omega::instance().getScene().get();
	
  const int id1 = I->getId1();
  const int id2 = I->getId2();

  if (geom.penetrationDepth<0) {
    return false;
  } else {
    const BodyContainer& bodies = *scene->bodies;
    
    const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
    const State& de2 = *static_cast<State*>(bodies[id2]->state.get());

    const Real& dt = scene->dt;

    if (I->isFresh(scene)) {
      phys.ukn = Vector3r(0,0,0);
      phys.uks = Vector3r(0,0,0);
      phys.shearForce = Vector3r(0,0,0);
      phys.normalForce = phys.kmn * geom.penetrationDepth * geom.normal;
    } else {
      Real An = 1.0 + phys.kkn * dt / (2.0 * phys.ckn);
      Real Bn = 1.0 - phys.kkn * dt / (2.0 * phys.ckn);
      Real Cn = dt / (2.0 * phys.ckn * An) + 1.0 / phys.kmn + dt / (2.0 * phys.cmn);
      Real Dn = dt / (2.0 * phys.ckn * An) - 1.0 / phys.kmn + dt / (2.0 * phys.cmn);

      Real As = 1.0 + phys.kks * dt / (2.0 * phys.cks);
      Real Bs = 1.0 -  phys.kks * dt / (2.0 * phys.cks);
      Real Cs = dt / (2.0 * phys.cks * As) + 1.0 / phys.kms + dt / (2.0 * phys.cms);
      Real Ds = dt / (2.0 * phys.cks * As) - 1.0 / phys.kms + dt / (2.0 * phys.cms);

      const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero(); 
      const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero(); 

      const Vector3r c1x = (geom.contactPoint - de1.pos);
      const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);

      const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel;
      const Vector3r normalVelocity = geom.normal.dot(relativeVelocity) * geom.normal;
      const Vector3r shearVelocity = relativeVelocity-normalVelocity;

      const Vector3r normalForceT = (normalVelocity * dt + phys.ukn * (1.0 - Bn / An) - phys.normalForce * Dn) / Cn;
      phys.ukn = (phys.ukn * Bn + (normalForceT + phys.normalForce) * dt / (2.0 * phys.ckn)) / An;
      phys.normalForce = normalForceT;

      const Vector3r shearForceT = -1 * (shearVelocity * dt + phys.uks * (1.0 - Bs / As) - phys.shearForce * Ds) / Cs;
      phys.uks = (phys.uks * Bs + (shearForceT + phys.shearForce) * dt / (2.0 * phys.cks)) / As;
      phys.shearForce = shearForceT;

      const Real maxFs = phys.normalForce.squaredNorm() * std::pow(phys.tangensOfFrictionAngle,2);
      if( phys.shearForce.squaredNorm() > maxFs ) {
		  
        const Real ratio = sqrt(maxFs) / phys.shearForce.norm();
        phys.shearForce *= ratio;
      }

      force = phys.normalForce + phys.shearForce;
      torque1 = -c1x.cross(force);
      torque2 =  c2x.cross(force);
      return true;
    }
  }
}

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.