yade-dev team mailing list archive

[gladky_anton@xxxxxxxxxxxxxxx]

Author: gladky_anton
Date: 2009-07-13 17:30:36 +0200 (Mon, 13 Jul 2009)
New Revision: 1861

Modified:
   trunk/pkg/dem/RockPM.cpp
   trunk/pkg/dem/RockPM.hpp
Log:
Changes on RockPM.* files


Modified: trunk/pkg/dem/RockPM.cpp
===================================================================
--- trunk/pkg/dem/RockPM.cpp	2009-07-13 09:41:46 UTC (rev 1860)
+++ trunk/pkg/dem/RockPM.cpp	2009-07-13 15:30:36 UTC (rev 1861)
@@ -28,30 +28,34 @@
 
 
 /********************** Law2_Dem3DofGeom_RockPMPhys_Rpm ****************************/
+CREATE_LOGGER(Law2_Dem3DofGeom_RockPMPhys_Rpm);
 
 void Law2_Dem3DofGeom_RockPMPhys_Rpm::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, MetaBody* rootBody){
 	Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get());
-	
 	RpmPhys* phys=static_cast<RpmPhys*>(ip.get());
-	
+
 	Real displN=geom->displacementN();
-	
-	if((displN<=0)){
-/*Normal Interaction*/
+	if(displN<=0){
+		/*Normal Interaction*/
 		phys->normalForce=phys->kn*displN*geom->normal;
 		Real maxFsSq=phys->normalForce.SquaredLength()*pow(phys->tanFrictionAngle,2);
 		Vector3r trialFs=phys->ks*geom->displacementT();
 		if(trialFs.SquaredLength()>maxFsSq){ geom->slipToDisplacementTMax(sqrt(maxFsSq)); trialFs*=sqrt(maxFsSq/(trialFs.SquaredLength()));} 
 		applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,rootBody);
-/*Normal Interaction_____*/
-		return;
-	} else if ((displN>0)&&(phys->isCohesive)){
-		if (displN<(phys->lengthMaxTension)) {
-			rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return;
+		/*Normal Interaction_____*/
+		if ((phys->isCohesive)&&(displN<(-phys->lengthMaxCompression))) {
+			LOG_WARN(displN<<"____"<<phys->normalForce);
 		}
-		rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return;
 	} else {
-		rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return;
+		if (phys->isCohesive) {
+			phys->normalForce=phys->kn*displN*geom->normal;
+			//LOG_WARN(displN<<"____"<<phys->normalForce);
+			applyForceAtContactPoint(phys->normalForce,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,rootBody);
+			return;
+		} else {
+			rootBody->interactions->requestErase(contact->getId1(),contact->getId2());
+			return;
+		}
 	}
 	
 }
@@ -63,6 +67,7 @@
 	if(interaction->interactionPhysics) return; 
 
 	Dem3DofGeom* contGeom=YADE_CAST<Dem3DofGeom*>(interaction->interactionGeometry.get());
+	
 	assert(contGeom);
 	
 	const shared_ptr<RpmMat>& rpm1=YADE_PTR_CAST<RpmMat>(pp1);
@@ -71,6 +76,8 @@
 	const shared_ptr<BodyMacroParameters>& elast1=static_pointer_cast<BodyMacroParameters>(pp1);
 	const shared_ptr<BodyMacroParameters>& elast2=static_pointer_cast<BodyMacroParameters>(pp2);
 
+	bool initCohesive = rpm1->initCohesive*rpm2->initCohesive;
+	
 	Real E12=2*elast1->young*elast2->young/(elast1->young+elast2->young);
 	Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
 	Real S12=Mathr::PI*pow(minRad,2);
@@ -82,13 +89,18 @@
 	contPhys->kn=contPhys->E*contPhys->crossSection;
 	contPhys->ks=contPhys->G*contPhys->crossSection;
 	
-	contPhys->lengthMaxCompression=S12*(rpm2->stressCompressMax)/(contPhys->kn);
-	contPhys->lengthMaxTension=S12*(rpm2->stressTensionMax)/(contPhys->kn);
+	contPhys->lengthMaxCompression=S12*(0.5*(rpm1->stressCompressMax+rpm2->stressCompressMax))/(contPhys->kn);
+	contPhys->lengthMaxTension=S12*(0.5*(rpm1->stressTensionMax+rpm2->stressTensionMax))/(contPhys->kn);
 	
-	//If 2 bodies were not damaged previously, and they belong to 1 example (stone), they can be cohesive
-	if ((rpm1->exampleNumber==rpm2->exampleNumber)&&(!(rpm1->isDamaged||rpm2->isDamaged))) {
+	initDistance = contGeom->displacementN();
+
+	if ((rpm1->exampleNumber==rpm2->exampleNumber)&&(initDistance<(contPhys->lengthMaxTension))&&(initCohesive)) {
 		contPhys->isCohesive=true;
+		//LOG_WARN("InitDistance="<<initDistance);
+		//LOG_WARN("lengthMaxCompression="<<contPhys->lengthMaxCompression);
+		//LOG_WARN("lengthMaxTension="<<contPhys->lengthMaxTension);
 	}
+
 	
 	interaction->interactionPhysics=contPhys;
 }

Modified: trunk/pkg/dem/RockPM.hpp
===================================================================
--- trunk/pkg/dem/RockPM.hpp	2009-07-13 09:41:46 UTC (rev 1860)
+++ trunk/pkg/dem/RockPM.hpp	2009-07-13 15:30:36 UTC (rev 1861)
@@ -39,6 +39,7 @@
 		FUNCTOR2D(Dem3DofGeom,RpmPhys);
 		REGISTER_CLASS_AND_BASE(Law2_Dem3DofGeom_RockPMPhys_Rpm,ConstitutiveLaw);
 		REGISTER_ATTRIBUTES(ConstitutiveLaw,/*nothing here*/);
+		DECLARE_LOGGER;	
 };
 REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_RockPMPhys_Rpm);
 
@@ -46,11 +47,11 @@
 class RpmMat: public BodyMacroParameters {
 	public:
 		int exampleNumber; //Number of "stone"
-		bool isDamaged;
+		bool initCohesive;
 		Real stressCompressMax, stressTensionMax; //Parameters for damaging
 
-		RpmMat(): exampleNumber(0.), isDamaged(false), stressCompressMax(0), stressTensionMax(0) {createIndex();};
-		REGISTER_ATTRIBUTES(BodyMacroParameters, (exampleNumber) (isDamaged) (stressCompressMax) (stressTensionMax));
+		RpmMat(): exampleNumber(0.), initCohesive(false), stressCompressMax(0), stressTensionMax(0) {createIndex();};
+		REGISTER_ATTRIBUTES(BodyMacroParameters, (exampleNumber) (initCohesive) (stressCompressMax) (stressTensionMax));
 		REGISTER_CLASS_AND_BASE(RpmMat,BodyMacroParameters);
 };
 REGISTER_SERIALIZABLE(RpmMat);
@@ -59,16 +60,16 @@
 class Ip2_RpmMat_RpmMat_RpmPhys: public InteractionPhysicsEngineUnit{
 	private:
 	public:
-		Real sigmaT;
+		Real sigmaT, initDistance;
 
 		Ip2_RpmMat_RpmMat_RpmPhys(){
-			// init to signaling_NaN to force crash if not initialized (better than unknowingly using garbage values)
-			sigmaT=3;
+			initDistance = 0;
 		}
 
 		virtual void go(const shared_ptr<PhysicalParameters>& pp1, const shared_ptr<PhysicalParameters>& pp2, const shared_ptr<Interaction>& interaction);
 		REGISTER_ATTRIBUTES(InteractionPhysicsEngineUnit,
 			(sigmaT)
+			(initDistance)
 		);
 
 		FUNCTOR2D(RpmMat,RpmMat);