yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 1828
committer: CWBoon <booncw@xxxxxxxxxxx>
branch nick: yade
timestamp: Mon 2009-11-30 17:17:35 -0500
message:
  BasicPM.cpp hpp BasicPMTest.py in pkg/dem/meta
added:
  pkg/dem/meta/BasicPM.cpp
  pkg/dem/meta/BasicPM.hpp
  pkg/dem/meta/BasicPMTest.py


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.

=== added file 'pkg/dem/meta/BasicPM.cpp'
--- pkg/dem/meta/BasicPM.cpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/meta/BasicPM.cpp	2009-11-30 22:17:35 +0000
@@ -0,0 +1,86 @@
+/*C.W.Boon @ NOV 2009 */
+
+#include"BasicPM.hpp"
+#include<yade/core/MetaBody.hpp>
+#include<yade/pkg-dem/BodyMacroParameters.hpp>
+#include<yade/pkg-dem/DemXDofGeom.hpp>
+
+
+
+YADE_PLUGIN((CundallStrackLaw)(Ip2_CSMat_CSMat_CSPhys));
+
+
+/********************** CundallStrackLaw ****************************/
+CREATE_LOGGER(CundallStrackLaw);
+
+void CundallStrackLaw::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, MetaBody* rootBody){
+	
+	Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get());
+	CSPhys* phys=static_cast<CSPhys*>(ip.get());
+	
+	/*Normal Force */
+	Real displN=geom->displacementN();
+	if(displN>0){rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return; }
+	phys->normalForce=phys->kn*displN*geom->normal; //normalForce is a vector
+	
+	/*ShearForce*/
+	Real maxFsSq=phys->normalForce.SquaredLength()*pow(phys->tangensOfFrictionAngle,2);
+	Vector3r trialFs=phys->ks*geom->displacementT();
+	if(trialFs.SquaredLength()>maxFsSq){ geom->slipToDisplacementTMax(sqrt(maxFsSq)); trialFs*=sqrt(maxFsSq/(trialFs.SquaredLength()));}
+	phys->shearForce=trialFs;
+	
+	applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,rootBody);
+}
+/********************CundallStrackLaw <end> ************************/
+
+
+
+/*****************************Ip2_CSMat_CSMat_CSPhys********************************/
+CREATE_LOGGER(Ip2_CSMat_CSMat_CSPhys);
+
+void Ip2_CSMat_CSMat_CSPhys::go(const shared_ptr<PhysicalParameters>& b1, const shared_ptr<PhysicalParameters>& b2 , const shared_ptr<Interaction>& interaction)
+{
+
+		
+		if(interaction->interactionPhysics) return;
+		//it could be other interactions by using if-else statements
+		Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(interaction->interactionGeometry.get());
+		assert(d3dg);
+
+			//const shared_ptr<CSMat>& csmat1=YADE_PTR_CAST<CSMat>(b1);
+			//const shared_ptr<CSMat>& csmat2=YADE_PTR_CAST<CSMat>(b2);
+
+			const shared_ptr<BodyMacroParameters>& sdec1 = YADE_PTR_CAST<BodyMacroParameters>(b1);
+			const shared_ptr<BodyMacroParameters>& sdec2 = YADE_PTR_CAST<BodyMacroParameters>(b2);
+			
+			shared_ptr<CSPhys> contactPhysics(new CSPhys()); //what if I declared it as constant?
+			
+			/* From interaction physics */
+			Real Ea 	= sdec1->young;
+			Real Eb 	= sdec2->young;
+			Real Va 	= sdec1->poisson;
+			Real Vb 	= sdec2->poisson;
+			Real fa 	= sdec1->frictionAngle;
+			Real fb 	= sdec2->frictionAngle;
+
+			/* From interaction geometry */
+			Real Da=d3dg->refR1>0?d3dg->refR1:d3dg->refR2; 
+			Real Db=d3dg->refR2>0?d3dg->refR2:d3dg->refR1; 
+			Vector3r normal=d3dg->normal;		  			
+			Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
+			Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
+
+			/* Pass values calculated from above to CSPhys */
+			contactPhysics->kn = Kn;
+			contactPhysics->ks = Ks;
+			contactPhysics->frictionAngle			= std::min(fa,fb); 
+			contactPhysics->tangensOfFrictionAngle		= std::tan(contactPhysics->frictionAngle); 
+			contactPhysics->prevNormal 			= normal;
+			
+
+}
+/*****************************Ip2_CSMat_CSMat_CSPhys <end>********************************/
+
+
+
+

=== added file 'pkg/dem/meta/BasicPM.hpp'
--- pkg/dem/meta/BasicPM.hpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/meta/BasicPM.hpp	2009-11-30 22:17:35 +0000
@@ -0,0 +1,64 @@
+/* C.W.Boon @ NOV 2009 */
+
+#pragma once
+
+#include<yade/pkg-dem/BodyMacroParameters.hpp> //Superclass for body properties
+#include<yade/pkg-common/InteractionPhysicsEngineUnit.hpp> //Superclass to link body and interaction properties
+#include<yade/pkg-common/NormalShearInteractions.hpp> //Superclass for interaction properties 
+#include<yade/pkg-common/ConstitutiveLaw.hpp> //Superclass for contact laws
+
+/* Contact Law */
+class CundallStrackLaw: public ConstitutiveLaw{
+	public:
+		
+		virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody);
+		FUNCTOR2D(Dem3DofGeom,CSPhys);
+		REGISTER_CLASS_AND_BASE(CundallStrackLaw,ConstitutiveLaw);
+		REGISTER_ATTRIBUTES(ConstitutiveLaw,/*nothing here*/);
+		DECLARE_LOGGER;	
+};
+REGISTER_SERIALIZABLE(CundallStrackLaw);
+
+/* This class stores body properties */
+/* Superclass:BodyMacroParameters has poisson, frictionAngle; ElasticBodyParameters has young's modulus */
+class CSMat: public BodyMacroParameters {
+	public:
+		CSMat(){createIndex();};
+		REGISTER_ATTRIBUTES(BodyMacroParameters, /*nothing here*/);
+		REGISTER_CLASS_AND_BASE(CSMat,BodyMacroParameters);
+		REGISTER_CLASS_INDEX(CSMat,BodyMacroParameters);
+};
+REGISTER_SERIALIZABLE(CSMat);
+
+/* This class links body and interaction properties */
+/* It does not store variables */
+class Ip2_CSMat_CSMat_CSPhys: public InteractionPhysicsEngineUnit{
+	public:
+		Ip2_CSMat_CSMat_CSPhys(){};
+		virtual void go(const shared_ptr<PhysicalParameters>& b1, const shared_ptr<PhysicalParameters>& b2, const shared_ptr<Interaction>& interaction);
+		REGISTER_ATTRIBUTES(InteractionPhysicsEngineUnit,/*nothing here*/);
+		FUNCTOR2D(CSMat,CSMat);
+		REGISTER_CLASS_AND_BASE(Ip2_CSMat_CSMat_CSPhys,InteractionPhysicsEngineUnit);
+		DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Ip2_CSMat_CSMat_CSPhys);
+
+
+/* This class stores interaction properties */
+/* Superclass:Normal NormalShearInteraction has kn, <vector> normalForce, ks, <vector> shearForce */
+class CSPhys: public NormalShearInteraction {
+	public:
+		CSPhys(){createIndex();};
+
+		// kn,ks,normal inherited from NormalShearInteraction
+		Real 		frictionAngle 			// angle of friction, according to Coulumb criterion
+				,tangensOfFrictionAngle
+				;	
+		Vector3r	prevNormal;			// unit normal of the contact plane.
+		
+	virtual ~CSPhys(){};
+	REGISTER_ATTRIBUTES(NormalShearInteraction,(frictionAngle)(tangensOfFrictionAngle));
+	REGISTER_CLASS_AND_BASE(CSPhys,NormalShearInteraction);
+	REGISTER_CLASS_INDEX(CSPhys,NormalShearInteraction);
+};
+REGISTER_SERIALIZABLE(CSPhys);

=== added file 'pkg/dem/meta/BasicPMTest.py'
--- pkg/dem/meta/BasicPMTest.py	1970-01-01 00:00:00 +0000
+++ pkg/dem/meta/BasicPMTest.py	2009-11-30 22:17:35 +0000
@@ -0,0 +1,65 @@
+#! file to test BasicPM.cpp
+# -*- coding: utf-8 -*-
+
+#!/usr/local/bin/yade-trunk -x
+# -*- coding: utf-8 -*-
+# -*- encoding=utf-8 -*-
+
+
+o=Omega()
+o.initializers=[
+	BoundingVolumeMetaEngine([InteractingSphere2AABB(),MetaInteractingGeometry2AABB()])
+]
+
+o.engines=[
+	
+	BexResetter(),
+	
+	BoundingVolumeMetaEngine([
+		InteractingSphere2AABB(),
+		MetaInteractingGeometry2AABB()
+	]),
+	
+	InsertionSortCollider(),
+	
+	InteractionGeometryMetaEngine([
+		ef2_Sphere_Sphere_Dem3DofGeom()
+	]),
+	
+	InteractionPhysicsMetaEngine([Ip2_CSMat_CSMat_CSPhys()]),
+	
+	CundallStrackLaw(),
+
+	
+	GravityEngine(gravity=[0,0,-9.81]),
+	
+	
+	
+	PhysicalActionApplier([
+		NewtonsForceLaw(),
+		NewtonsMomentumLaw(),
+	]),
+
+	PhysicalParametersMetaEngine([LeapFrogPositionIntegrator()]),
+
+	PhysicalParametersMetaEngine([LeapFrogOrientationIntegrator()])
+
+]
+
+
+from yade import utils
+
+for b in o.bodies:
+  b.mask = 1
+
+o.bodies.append(utils.sphere([0,0,2],1,dynamic = False, color=[0,1,0],young=30e9,poisson=.3,density=2400))
+o.bodies.append(utils.sphere([0,0,6],1,color=[0,0,1],young=30e9,poisson=.3,density=2400))
+
+o.dt=.2*utils.PWaveTimeStep()
+
+
+from yade import qt
+qt.Controller()
+qt.View()
+
+