yade-users team mailing list archive

[behzad <question255780@xxxxxxxxxxxxxxxxxxxxx>]

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

behzad posted a new comment:

IMPLEMENTATION


#include<yade/pkg/dem/Burgers.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>


//********************** Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys ****************************/


void Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys::go(const shared_ptr<Material>& b1    
                                        , const shared_ptr<Material>& b2 
                                        , const shared_ptr<Interaction>& interaction)
{
	CohBurgersMat* sdec1 = static_cast<CohBurgersMat*>(b1.get());
	CohBurgersMat* sdec2 = static_cast<CohBurgersMat*>(b2.get());
	ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());

	if (setCohesionNow && cohesionDefinitionIteration==-1) cohesionDefinitionIteration=scene->iter;
	if (setCohesionNow && cohesionDefinitionIteration!=-1 && cohesionDefinitionIteration!=scene->iter) {
		cohesionDefinitionIteration = -1;
		setCohesionNow = 0;}

	if (geom) {
		if (!interaction->phys) {
			interaction->phys = shared_ptr<CohBurgersPhys>(new CohBurgersPhys());
			CohBurgersPhys* contactPhysics = YADE_CAST<CohBurgersPhys*>(interaction->phys.get());
			
			Real Da 		= geom->radius1;
			Real Db 		= geom->radius2;

			Real BurEmn		=sdec1->Em;
			Real BurEkn		=sdec1->Ek;
			Real BurCmn		=sdec1->Cm;
			Real BurCkn		=sdec1->Ck;

			Real BurEms		=sdec1->poissonRatio*BurEmn;
			Real BurEks		=sdec1->poissonRatio*BurEkn;
			Real BurCms		=sdec1->poissonRatio*BurCmn;
			Real BurCks		=sdec1->poissonRatio*BurCkn;

			Real poissonRatio	=sdec1->poissonRatio
			Real frictionAngle	=sdec1->frictionAngle
			Real normalAdhesion	=sdec1->normalAdhesion
			Real shearAdhesion	=sdec1->shearAdhesion



//********************** Ip2_FrictMat_CohBurgersMat_CohFrictPhys ****************************/


CREATE_LOGGER(Ip2_FrictMat_CohBurgersMat_CohFrictPhys);
void Ip2_FrictMat_CohBurgersMat_CohFrictPhys::go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction){
	TIMING_DELTAS_START();
	const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(pp1);
	const shared_ptr<CohBurgersMat>& mat2 = YADE_PTR_CAST<CohBurgersMat>(pp2);
	Ip2_FrictMat_CohBurgersMat_CohFrictPhys().go(mat1,mat2,interaction);
	TIMING_DELTAS_CHECKPOINT("end of Ip2_CohFritPhys");




//********************** CohesiveBurgersContactLaw ****************************/


YADE_PLUGIN((CohesiveBurgersContactLaw)(Law2_ScGeom_CohBurgersPhys)(CohBurgerstMat)(CohBurgersPhys)(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys));
CREATE_LOGGER(Law2_ScGeom_CohBurgersPhys);


void CohesiveBurgersContactLaw::action()

bool Law2_ScGeom_CohBurgersPhys::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
	int id1 = contact->getId1(), id2 = contact->getId2();

	ScGeom*    geom= static_cast<ScGeom*>(ig.get());
	FrictPhys* phys = static_cast<CohBurgersPhys*>(ip.get());
	if(geom->penetrationDepth <0){
		if (neverErase) {
			Real Da  = geom->radius1;
                  	Real Db  = geom->radius2;
                  
			Vector3r& shearForce = Phys->shearForce;
                  	Vector3r& normalForce = Phys->normalForce;
                 	Vector3r& ukkn = Phys->ukkn;
                 	Vector3r& shearForceT1 = Phys->shearForceT1;
                 	Vector3r& normalForceT1 = Phys->normalForceT1;
                  	Vector3r& ukks = Phys->ukks;
		else return false;


                  int& Ite        = Phys->Ite;
                  
                  Real Dinit = Da+Db;
                  Real knk = BurEkn*Dinit;
                  Real cnk = BurCkn*Dinit;
                  Real knm = BurEmn*Dinit;
                  Real cnm = BurCmn*Dinit;
                  Real ksk = BurEks*Dinit/(2.0+2.0*BurPoi);
                  Real csk = BurCks*Dinit/(2.0+2.0*BurPoi);
                  Real ksm = BurEms*Dinit/(2.0+2.0*BurPoi);
                  Real csm = BurCms*Dinit/(2.0+2.0*BurPoi);
                  
                  Real un=geom->penetrationDepth;
                  Real burnA=1.0+knk*dt/(2.0*cnk);
                  Real burnB=1.0-knk*dt/(2.0*cnk);
                  Real burnC=dt/(2.0*cnk*burnA)+1.0/knm+dt/(2.0*cnm);
                  Real burnD=dt/(2.0*cnk*burnA)-1.0/knm+dt/(2.0*cnm);
                  Real bursA=1.0+ksk*dt/(2.0*csk);
                  Real bursB=1.0-ksk*dt/(2.0*csk);
                  Real bursC=dt/(2.0*csk*bursA)+1.0/ksm+dt/(2.0*csm);
                  Real bursD=dt/(2.0*csk*bursA)-1.0/ksm+dt/(2.0*csm);

                  if (contact->isNew)
                        {   ukkn=Vector3r(0,0,0);
                            ukks=Vector3r(0,0,0);
                            shearForce=Vector3r(0,0,0);
                            normalForce=knm*std::max(un,(Real)
0)*geom->normal; 
                         }
                 Vector3r axis;
                 Real angle;
                 axis = Phys->prevNormal.Cross(geom->normal);
                 shearForce -= shearForce.Cross(axis);
                 Vector3r summaryAngularVelocity(0,0,0);
                 if (isDynamic1) summaryAngularVelocity += de1->angularVelocity;
                 if (isDynamic2) summaryAngularVelocity += de2->angularVelocity;
                 angle = dt*0.5*geom->normal.Dot(summaryAngularVelocity);
                 axis = angle*geom->normal;
                 shearForce -= shearForce.Cross(axis);
                 Vector3r x = geom->contactPoint;
                 Vector3r c1x = (x - de1-se3.position);
                 Vector3r c2x = (x - de2-se3.position);
                 Vector3r _c1x_= (isDynamic1) ?
geom->radius1*geom->normal : x -de1->zeroPoint;
                 Vector3r _c2x_= (isDynamic2) ? -
geom->radius2*geom->normal : x -
de2->zeroPoint;
                 Vector3r relativeVelocity= (de2-
>velocity+de2->angularVelocity.Cross(_c2x_))-(de1->velocity+de1-
>angularVelocity.Cross(_c1x_));
                 Vector3r normalVelocity = 
std::abs(geom->normal.Dot(relativeVelocity))*geom->normal;
                 Vector3r shearVelocity =relativeVelocity-normalVelocity;
                 Vector3r shearDisplacement = shearVelocity*dt;
                 Vector3r normalDisplacement= normalVelocity*dt;

                       if (contact->isNew)
Ite=Omega::instance().getCurrentIteration();
                  normalForceT1 = (normalDisplacement+ukkn*(1.0-
burnB/burnA)-normalForce*burnD)/burnC;
                  ukkn =
(ukkn*burnB+(normalForceT1+normalForce)*dt/(2.0*cnk))/burnA;
                  normalForce = normalForceT1;
                  
                  shearForceT1 = -1.0*(shearDisplacement+ukks*(1.0-
bursB/bursA)+bursD*shearForce)/bursC;
                  ukks = (ukks*bursB-
(shearForceT1+shearForce)*dt/(2.0*csk))/bursA;
                  shearForce = shearForceT1;
                  
                  Real maxFs = normalForce.SquaredLength() *
std::pow(Phys->tangensOfFrictionAngle,2);
                  if( shearForce.SquaredLength() > maxFs )
                  {
					   maxFs = Mathr::Sqrt(maxFs) / shearForce.Length();
                       shearForce *= maxFs;
                  }
                  
                  Vector3r f = normalForce +
shearForce;
                  static_cast<Force*> ( ncb->physicalActions-
>find( id1 , actionForceIndex).get() )->force -= f;
                  static_cast<Force*> ( ncb->physicalActions-
>find( id2 , actionForceIndex ).get() )->force += f;
                  static_cast<Momentum*>( ncb->physicalActions-
>find( id1 , actionMomentumIndex ).get() )->momentum -= c1x.Cross(f);
                  static_cast<Momentum*>( ncb->physicalActions-
>find( id2 , actionMomentumIndex ).get() )->momentum += c2x.Cross(f);
                  Phys->prevNormal =geom->normal;
            }
      }
}
YADE_PLUGIN();

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