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);
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
