← Back to branch summary

~ubuntu-branches/debian/jessie/yade/jessie

« back to all changes in this revision

Viewing changes to pkg/common/Cylinder.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2014-08-04 19:34:58 UTC
  • mfrom: (1.1.11)
  • Revision ID: package-import@ubuntu.com-20140804193458-cw8qhnujxe9wzi15
Tags: 1.11.0-1
* [a0600ae] Imported Upstream version 1.11.0
* [a3055e0] Do not use parallel build on kfreebsd-amd64 and s390x.
* [f86b405] Remove applied patches.

Show diffs side-by-side

added added

removed removed

Lines of Context:
pkg/common/Cylinder.cpp
396
396
                Vector3r N=a.cross(b);
397
397
                Vector3r normal;
398
398
                if(N.norm()>1e-14){
399
 
                        dist=abs(N.dot(B-A)/(N.norm()));        //distance between the two LINES.
 
399
                        dist=std::abs(N.dot(B-A)/(N.norm()));   //distance between the two LINES.
400
400
                        //But we need to check that the intersection point is inside the two SEGMENTS ...
401
401
                        //Projection of B to have a common plan between the two segments.
402
402
                        Vector3r projB1=B+dist*(N/(N.norm())) , projB2=B-dist*(N/(N.norm()));
403
403
                        Real distB1A=(projB1-A).norm() , distB2A=(projB2-A).norm() ;
404
404
                        Vector3r projB=(distB1A<=distB2A)*projB1 + (distB1A>distB2A)*projB2;
405
405
                        int b1=0, b2=1; //base vectors used to compute the segment intersection (if N is aligned with an axis, we can't use this axis)
406
 
                        if(abs(N[1])<1e-14 && abs(N[2])<1e-14){b1=1;b2=2;}
407
 
                        if(abs(N[0])<1e-14 && abs(N[2])<1e-14){b1=0;b2=2;}
 
406
                        if(std::abs(N[1])<1e-14 && std::abs(N[2])<1e-14){b1=1;b2=2;}
 
407
                        if(std::abs(N[0])<1e-14 && std::abs(N[2])<1e-14){b1=0;b2=2;}
408
408
                        Real det=a[b1]*b[b2]-a[b2]*b[b1];
409
 
                        if(abs(det)>1e-14){     //Check if the two segments are intersected (using k and m)
 
409
                        if(std::abs(det)>1e-14){        //Check if the two segments are intersected (using k and m)
410
410
                                k = (b[b2]*(projB[b1]-A[b1])+b[b1]*(A[b2]-projB[b2]))/det;
411
411
                                m = (a[b1]*(-projB[b2]+A[b2])+a[b2]*(projB[b1]-A[b1]))/det;
412
412
                                if( k<0.0 || k>=1.0 || m<0.0 || m>=1.0 ) {      //so they are not intersected
658
658
        }
659
659
}
660
660
 
661
 
void Law2_CylScGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
 
661
bool Law2_CylScGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
662
662
        int id1 = contact->getId1(), id2 = contact->getId2();
663
663
 
664
664
        CylScGeom* geom= static_cast<CylScGeom*>(ig.get());
667
667
                if (neverErase) {
668
668
                        phys->shearForce = Vector3r::Zero();
669
669
                        phys->normalForce = Vector3r::Zero();}
670
 
                else    scene->interactions->requestErase(contact);
671
 
                return;}
 
670
                else return false;}
672
671
        if (geom->isDuplicate) {
673
672
                if (id2!=geom->trueInt) {
674
673
                        //cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
675
 
                        if (geom->isDuplicate==2) {/*cerr<<"erase duplicate "<<id1<<" "<<id2<<endl;*/scene->interactions->requestErase(contact);}
676
 
                return;}
 
674
                        if (geom->isDuplicate==2) return false;}
677
675
        }
678
676
        Real& un=geom->penetrationDepth;
679
677
        phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
721
719
                scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
722
720
                scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
723
721
        }
 
722
        return true;
724
723
}
725
724
 
726
725
 
727
 
void Law2_CylScGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact) {
 
726
bool Law2_CylScGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact) {
728
727
 
729
728
    int id1 = contact->getId1(), id2 = contact->getId2();
730
729
 
738
737
    if (geom->isDuplicate) {
739
738
                if (id2!=geom->trueInt) {
740
739
                        //cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
741
 
                        if (geom->isDuplicate==2) {/*cerr<<"erase duplicate coh "<<id1<<" "<<id2<<endl;*/scene->interactions->requestErase(contact);}
742
 
                return;}
 
740
                        if (geom->isDuplicate==2) return false;}
743
741
        }
744
742
 
745
743
    if (currentContactPhysics->fragile && (-Fn)> currentContactPhysics->normalAdhesion) {
746
744
        // BREAK due to tension
747
 
        scene->interactions->requestErase(contact);
 
745
        return false;
748
746
    } else {
749
747
        if ((-Fn)> currentContactPhysics->normalAdhesion) {//normal plasticity
750
748
            Fn=-currentContactPhysics->normalAdhesion;
751
749
            currentContactPhysics->unp = un+currentContactPhysics->normalAdhesion/currentContactPhysics->kn;
752
 
            if (currentContactPhysics->unpMax && currentContactPhysics->unp<currentContactPhysics->unpMax)
753
 
                scene->interactions->requestErase(contact);
 
750
            if (currentContactPhysics->unpMax && currentContactPhysics->unp<currentContactPhysics->unpMax) return false;
754
751
        }
755
752
        currentContactPhysics->normalForce = Fn*geom->normal;
756
753
        Vector3r& shearForce = geom->rotate(currentContactPhysics->shearForce);
794
791
        }
795
792
        //applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
796
793
    }
 
794
    return true;
797
795
}
798
796
 
799
797
 
800
798
 
801
 
void Law2_ChCylGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
 
799
bool Law2_ChCylGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
802
800
  int id1 = contact->getId1(), id2 = contact->getId2();
803
801
 
804
802
    ChCylGeom6D* geom= YADE_CAST<ChCylGeom6D*>(ig.get());
816
814
    if (contact->isFresh(scene)) shearForce   = Vector3r::Zero();                       //contact nouveau => force tengentielle = 0,0,0
817
815
    Real un     = geom->penetrationDepth;                               //un : interpenetration
818
816
    Real Fn    = currentContactPhysics->kn*(un-currentContactPhysics->unp);             //Fn : force normale
819
 
    /*if (geom->isDuplicate) {
820
 
                if (id2!=geom->trueInt) {
821
 
                        //cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
822
 
                        if (geom->isDuplicate==2) {cerr<<"erase duplicate coh "<<id1<<" "<<id2<<endl;scene->interactions->requestErase(contact);}
823
 
                return;}
824
 
        }
825
 
*/
826
817
    
827
 
    if (currentContactPhysics->fragile && (-Fn)> currentContactPhysics->normalAdhesion) {
828
 
        // BREAK due to tension
829
 
        scene->interactions->requestErase(contact);
830
 
    } else {
 
818
    if (currentContactPhysics->fragile && (-Fn)> currentContactPhysics->normalAdhesion) return false; // BREAK due to tension
 
819
    else {
831
820
        if ((-Fn)> currentContactPhysics->normalAdhesion) {//normal plasticity
832
821
            Fn=-currentContactPhysics->normalAdhesion;
833
822
            currentContactPhysics->unp = un+currentContactPhysics->normalAdhesion/currentContactPhysics->kn;
834
823
            if (currentContactPhysics->unpMax && currentContactPhysics->unp<currentContactPhysics->unpMax)
835
 
                scene->interactions->requestErase(contact);
836
 
        }
 
824
                return false;
 
825
        }
837
826
    
838
827
        
839
828
        currentContactPhysics->normalForce = Fn*geom->normal;
885
874
                scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
886
875
        }
887
876
    }
 
877
    return true;
888
878
}
889
879
 
890
880