yade-users team mailing list archive

[Saeed <question681991@xxxxxxxxxxxxxxxxxxxxx>]

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

Saeed gave more information on the question:
Hi there
I think there is a problem with file 'TwoPhaseFlowEngine.cpp' to calculate rmin and rmax :
A portion of this file (engine) that performs radius calculations is listed below:

----------------------
double TwoPhaseFlowEngine::computeMSPRcByPosRadius(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC)
{
    double e[3]; //edges of triangulation
    double g[3]; //gap radius between solid
        
    e[0] = (posB-posC).norm();
    e[1] = (posC-posA).norm();
    e[2] = (posB-posA).norm();
    g[0] = ((e[0]-rB-rC)>0) ? 0.5*(e[0]-rB-rC):0 ;
    g[1] = ((e[1]-rC-rA)>0) ? 0.5*(e[1]-rC-rA):0 ;
    g[2] = ((e[2]-rA-rB)>0) ? 0.5*(e[2]-rA-rB):0 ;
    
    double rmin= (std::max(g[0],std::max(g[1],g[2]))==0) ? 1.0e-11:std::max(g[0],std::max(g[1],g[2])) ;
    double rmax= computeEffRcByPosRadius(posA, rA, posB, rB, posC, rC);
    if(rmin>rmax) { cerr<<"WARNING! rmin>rmax. rmin="<<rmin<<" ,rmax="<<rmax<<endl; }
    
    double deltaForceRMin = computeDeltaForce(posA, rA, posB, rB, posC, rC, rmin);
    double deltaForceRMax = computeDeltaForce(posA, rA, posB, rB, posC, rC, rmax);
    double effPoreRadius;
    
    if(deltaForceRMin>deltaForceRMax) { effPoreRadius=rmax; }
    else if(deltaForceRMax<0) { effPoreRadius=rmax; }
    else if(deltaForceRMin>0) { effPoreRadius=rmin; }
    else { effPoreRadius=bisection(posA, rA, posB, rB, posC, rC, rmin,rmax); }
    return effPoreRadius;
}
----------------------------------------

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.