yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2488
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Thu 2010-10-14 15:58:47 +0200
message:
  1. Make forgotten class members public (sorry, Anton)
modified:
  pkg/common/Dispatching.hpp
  pkg/dem/DomainLimiter.cpp
  pkg/dem/L3Geom.cpp
  pkg/dem/L3Geom.hpp


--
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

=== modified file 'pkg/common/Dispatching.hpp'
--- pkg/common/Dispatching.hpp	2010-10-13 16:23:08 +0000
+++ pkg/common/Dispatching.hpp	2010-10-14 13:58:47 +0000
@@ -70,6 +70,11 @@
 		addForce(id1, force,scene); addTorque(id1, (contactPoint-pos1).cross(force),scene);
 		addForce(id2,-force,scene); addTorque(id2,-(contactPoint-pos2).cross(force),scene);
 	}
+	/*! Convenience function to apply force and torque from one force, given branch vectors from body centroids (x1,x2) to contact point cp, i.e. cp-x1, cp-x2. */
+	void applyForceAtBranch(const Vector3r& force, const Body::id_t id1, const Vector3r& cx1, const Body::id_t id2, const Vector3r& cx2){
+		addForce(id1, force,scene); addTorque(id1,cx1.cross( force),scene);
+		addForce(id1,-force,scene); addTorque(id1,cx2.cross(-force),scene);
+	}
 	YADE_CLASS_BASE_DOC(LawFunctor,Functor,"Functor for applying constitutive laws on :yref:`interactions<Interaction>`.");
 };
 REGISTER_SERIALIZABLE(LawFunctor);

=== modified file 'pkg/dem/DomainLimiter.cpp'
--- pkg/dem/DomainLimiter.cpp	2010-10-14 08:18:42 +0000
+++ pkg/dem/DomainLimiter.cpp	2010-10-14 13:58:47 +0000
@@ -22,6 +22,7 @@
 
 #include<yade/pkg-dem/DemXDofGeom.hpp>
 #include<yade/pkg-dem/ScGeom.hpp>
+#include<yade/pkg-dem/L3Geom.hpp>
 #include<yade/pkg-common/NormShearPhys.hpp>
 #include<yade/pkg-common/LinearInterpolate.hpp>
 #include<yade/lib-pyutil/gil.hpp>
@@ -61,10 +62,11 @@
 	GenericSpheresContact* gsc=dynamic_cast<GenericSpheresContact*>(I->geom.get());
 	ScGeom* scGeom=dynamic_cast<ScGeom*>(I->geom.get());
 	Dem3DofGeom* d3dGeom=dynamic_cast<Dem3DofGeom*>(I->geom.get());
+	L3Geom* l3Geom=dynamic_cast<L3Geom*>(I->geom.get());
 	//NormShearPhys* phys=dynamic_cast<NormShearPhys*>(I->phys.get());			//Disabled because of warning
 	if(!gsc) throw std::invalid_argument("LawTester: IGeom of "+strIds+" not a GenericSpheresContact.");
-	if(!scGeom && !d3dGeom) throw std::invalid_argument("LawTester: IGeom of "+strIds+" is neither ScGeom, nor Dem3DofGeom.");
-	assert(!((bool)scGeom && (bool)d3dGeom)); // nonsense
+	if(!scGeom && !d3dGeom && !l3Geom) throw std::invalid_argument("LawTester: IGeom of "+strIds+" is neither ScGeom, nor Dem3DofGeom, nor L3Geom.");
+	assert(!((bool)scGeom && (bool)d3dGeom) && (bool)l3Geom); // nonsense
 	// get body objects
 	State *state1=Body::byId(id1,scene)->state.get(), *state2=Body::byId(id2,scene)->state.get();
 	scene->forces.sync();
@@ -78,35 +80,38 @@
 		else{ LOG_INFO("Running doneHook: "<<doneHook);	pyRunString(doneHook);}
 		return;
 	}
-
-	axX=gsc->normal; /* just in case */ axX.normalize();
-	if(doInit){ // initialization of the new interaction
-		// take vector in the y or z direction, depending on its length; arbitrary, but one of them is sure to be non-zero
-		axY=axX.cross(axX[1]>axX[2]?Vector3r::UnitY():Vector3r::UnitZ());
-		axY.normalize();
-		axZ=axX.cross(axY);
-		LOG_DEBUG("Initial axes x="<<axX<<", y="<<axY<<", z="<<axZ);
-		renderLength=.5*(gsc->refR1+gsc->refR2);
-		refLength=gsc->refR1+gsc->refR2;
-	} else { // udpate of an existing interaction
-		if(scGeom){
-			scGeom->rotate(axY); scGeom->rotate(axZ);
-			scGeom->rotate(shearTot);
-			shearTot+=scGeom->shearIncrement();
-			ptGeom=Vector3r(-scGeom->penetrationDepth,shearTot.dot(axY),shearTot.dot(axZ));
-		}
-		else{ // d3dGeom
-			throw runtime_error("LawTester: Dem3DofGeom not yet supported.");
-			// essentially copies code from ScGeom, which is not very nice indeed; oh well…
-			Vector3r vRel=(state2->vel+state2->angVel.cross(-gsc->refR2*gsc->normal))-(state1->vel+state1->angVel.cross(gsc->refR1*gsc->normal));
-		}
+	/* initialize or update local axes and trsf */
+	if(!l3Geom){
+		axX=gsc->normal; /* just in case */ axX.normalize();
+		if(doInit){ // initialization of the new interaction
+			// take vector in the y or z direction, depending on its length; arbitrary, but one of them is sure to be non-zero
+			axY=axX.cross(axX[1]>axX[2]?Vector3r::UnitY():Vector3r::UnitZ());
+			axY.normalize();
+			axZ=axX.cross(axY);
+			LOG_DEBUG("Initial axes x="<<axX<<", y="<<axY<<", z="<<axZ);
+		} else { // udpate of an existing interaction
+			if(scGeom){
+				scGeom->rotate(axY); scGeom->rotate(axZ);
+				scGeom->rotate(shearTot);
+				shearTot+=scGeom->shearIncrement();
+				ptGeom=Vector3r(-scGeom->penetrationDepth,shearTot.dot(axY),shearTot.dot(axZ));
+			}
+			else{ // d3dGeom
+				throw runtime_error("LawTester: Dem3DofGeom not yet supported.");
+				// essentially copies code from ScGeom, which is not very nice indeed; oh well…
+				Vector3r vRel=(state2->vel+state2->angVel.cross(-gsc->refR2*gsc->normal))-(state1->vel+state1->angVel.cross(gsc->refR1*gsc->normal));
+			}
+		}
+		// update the transformation
+		// the matrix is orthonormal, since axX, axY are normalized and and axZ is their corss-product
+		trsf.row(0)=axX; trsf.row(1)=axY; trsf.row(2)=axZ;
+	} else {
+		trsf=l3Geom->trsf;
 	}
-	// update the transformation
-	// the matrix is orthonormal, since axX, axY are normalized and and axZ is their corss-product
-	trsf.row(0)=axX; trsf.row(1)=axY; trsf.row(2)=axZ;
 	trsfQ=Quaternionr(trsf);
 	contPt=gsc->contactPoint;
-	Real refLength=gsc->refR1+gsc->refR2;
+	refLength=gsc->refR1+gsc->refR2;
+	renderLength=.5*refLength;
 	
 	// here we go ahead, finally
 	Vector3r val=linearInterpolate<Vector3r,int>(step,_pathT,_pathV,_interpPos);

=== modified file 'pkg/dem/L3Geom.cpp'
--- pkg/dem/L3Geom.cpp	2010-10-14 12:36:04 +0000
+++ pkg/dem/L3Geom.cpp	2010-10-14 13:58:47 +0000
@@ -2,7 +2,7 @@
 #include<yade/pkg-dem/L3Geom.hpp>
 #include<yade/pkg-common/Sphere.hpp>
 
-YADE_PLUGIN((L3Geom)(Ig2_Sphere_Sphere_L3Geom_Inc));
+YADE_PLUGIN((L3Geom)(Ig2_Sphere_Sphere_L3Geom_Inc)(Law2_L3Geom_FrictPhys_Linear));
 
 L3Geom::~L3Geom(){}
 
@@ -105,4 +105,13 @@
 };
 
 
+void Law2_L3Geom_FrictPhys_Linear::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
+	const shared_ptr<L3Geom> geom=static_pointer_cast<L3Geom>(ig); const shared_ptr<FrictPhys> phys=static_pointer_cast<FrictPhys>(ip);
+	// compute force:
+	Vector3r localF=geom->u.cwise()*Vector3r(phys->kn,phys->ks,phys->ks);
+	// transform back to global coords
+	Vector3r globalF=Quaternionr(geom->trsf).conjugate()*localF;
+	// apply force and torque
+	applyForceAtBranch(globalF,I->getId1(),geom->normal*(geom->refR1+.5*geom->u[0]),I->getId2(),-geom->normal*(geom->refR1+.5*geom->u[0]));
+}
 

=== modified file 'pkg/dem/L3Geom.hpp'
--- pkg/dem/L3Geom.hpp	2010-10-14 12:36:04 +0000
+++ pkg/dem/L3Geom.hpp	2010-10-14 13:58:47 +0000
@@ -1,15 +1,18 @@
 
 #include<yade/core/IGeom.hpp>
+#include<yade/core/IPhys.hpp>
 #include<yade/core/Shape.hpp>
 #include<yade/core/State.hpp>
 #include<yade/pkg-common/Dispatching.hpp>
 #include<yade/pkg-dem/DemXDofGeom.hpp>
+#include<yade/pkg-dem/FrictPhys.hpp>
 
-class L3Geom: public GenericSpheresContact{
+struct L3Geom: public GenericSpheresContact{
 	virtual ~L3Geom();
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR(L3Geom,GenericSpheresContact,"Geometry of contact given in local coordinates with 3 degress of freedom: normal and two in shear plane. [experimental]",
 		((Vector3r,u,Vector3r::Zero(),,"Displacement components, in local coordinates. |yupdate|"))
-		((Vector3r,u0,Vector3r::Zero(),,"Zero displacement value; u0 will be always subtracted from the *geometrical* displacement by appropriate :yref:`IGeomFunctor`, resulting in *u*. This value can be changed for instance\n\n#. by :yref:`IGeomFunctor`, e.g. to take in account large shear displacement value unrepresentable by underlying geomeric algorithm based on quaternions)\n#. by :yref:`LawFunctor`, to account for normal equilibrium position different from zero geometric overlap (set once, just after the interaction is created)\n#. by :yref:`LawFunctor` to account for plastic slip.\n\n.. note:: Never set an absolute value of *u0*, only increment, since both :yref:`IGeomFunctor` and :yref:`LawFunctor` use it. If you need to keep track of plastic deformation, store it in :yref:`IPhys` isntead (this might be changed: have *u0* for :yref:`LawFunctor` exclusively, and a separate value stored (when that is needed) inside classes deriving from :yref:`L3Geom`."))
+		/*((Vector3r,u0,Vector3r::Zero(),,"Zero displacement value; u0 will be always subtracted from the *geometrical* displacement by appropriate :yref:`IGeomFunctor`, resulting in *u*. This value can be changed for instance\n\n#. by :yref:`IGeomFunctor`, e.g. to take in account large shear displacement value unrepresentable by underlying geomeric algorithm based on quaternions)\n#. by :yref:`LawFunctor`, to account for normal equilibrium position different from zero geometric overlap (set once, just after the interaction is created)\n#. by :yref:`LawFunctor` to account for plastic slip.\n\n.. note:: Never set an absolute value of *u0*, only increment, since both :yref:`IGeomFunctor` and :yref:`LawFunctor` use it. If you need to keep track of plastic deformation, store it in :yref:`IPhys` isntead (this might be changed: have *u0* for :yref:`LawFunctor` exclusively, and a separate value stored (when that is needed) inside classes deriving from :yref:`L3Geom`."))
+		*/
 		/* Is it better to store trsf as Matrix3 or Quaternion?
 			* Quaternions are much easier to re-normalize, which we should do to avoid numerical drift.
 			* Multiplication of vector with quaternion is internally done by converting to matrix first, anyway
@@ -31,7 +34,7 @@
 }
 */
 
-class Ig2_Sphere_Sphere_L3Geom_Inc: public IGeomFunctor{
+struct Ig2_Sphere_Sphere_L3Geom_Inc: public IGeomFunctor{
 		virtual bool go(const shared_ptr<Shape>& s1, const shared_ptr<Shape>& s2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& I);
 
 	YADE_CLASS_BASE_DOC_ATTRS(Ig2_Sphere_Sphere_L3Geom_Inc,IGeomFunctor,"Incrementally compute :yref:`L3Geom` for contact of 2 spheres.\n\n.. note:: The initial value of *u[0]* (normal displacement) might be non-zero, with or without *distFactor*, since it is given purely by sphere's geometry. If you want to set \"equilibrium distance\", do it in the contact law as explained in :yref:`L3Geom.u0`.",
@@ -42,3 +45,10 @@
 	DEFINE_FUNCTOR_ORDER_2D(Sphere,Sphere);
 };
 REGISTER_SERIALIZABLE(Ig2_Sphere_Sphere_L3Geom_Inc);
+
+struct Law2_L3Geom_FrictPhys_Linear: public LawFunctor{
+	virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+	FUNCTOR2D(L3Geom,FrictPhys);
+	YADE_CLASS_BASE_DOC(Law2_L3Geom_FrictPhys_Linear,LawFunctor,"Basic law for testing :yref:`L3Geom` -- linear in both normal and shear sense, without slip or breakage.");
+};
+REGISTER_SERIALIZABLE(Law2_L3Geom_FrictPhys_Linear);