yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2489
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Thu 2010-10-14 19:21:34 +0200
message:
  1. Enhance LawTester to work with L3Geom, fix some bugs
  2. Add uN and uT references inside L3Geom, for respective components of the deformation vector
modified:
  pkg/common/Dispatching.hpp
  pkg/dem/DomainLimiter.cpp
  pkg/dem/L3Geom.cpp
  pkg/dem/L3Geom.hpp
  scripts/test/law-test.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

=== modified file 'pkg/common/Dispatching.hpp'
--- pkg/common/Dispatching.hpp	2010-10-14 13:58:47 +0000
+++ pkg/common/Dispatching.hpp	2010-10-14 17:21:34 +0000
@@ -73,7 +73,7 @@
 	/*! 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);
+		addForce(id2,-force,scene); addTorque(id2,cx2.cross(-force),scene);
 	}
 	YADE_CLASS_BASE_DOC(LawFunctor,Functor,"Functor for applying constitutive laws on :yref:`interactions<Interaction>`.");
 };

=== modified file 'pkg/dem/DomainLimiter.cpp'
--- pkg/dem/DomainLimiter.cpp	2010-10-14 13:58:47 +0000
+++ pkg/dem/DomainLimiter.cpp	2010-10-14 17:21:34 +0000
@@ -66,12 +66,14 @@
 	//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 && !l3Geom) throw std::invalid_argument("LawTester: IGeom of "+strIds+" is neither ScGeom, nor Dem3DofGeom, nor L3Geom.");
-	assert(!((bool)scGeom && (bool)d3dGeom) && (bool)l3Geom); // nonsense
+	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();
-	if(state1->blockedDOFs!=State::DOF_RXRYRZ) { LOG_INFO("Blocking rotational DOFs for #"<<id1); state1->blockedDOFs=State::DOF_RXRYRZ;}
-	if(state2->blockedDOFs!=State::DOF_RXRYRZ) { LOG_INFO("Blocking rotational DOFs for #"<<id2); state2->blockedDOFs=State::DOF_RXRYRZ;}
+	//if(state1->blockedDOFs!=State::DOF_RXRYRZ) { LOG_INFO("Blocking rotational DOFs for #"<<id1); state1->blockedDOFs=State::DOF_RXRYRZ;}
+	//if(state2->blockedDOFs!=State::DOF_RXRYRZ) { LOG_INFO("Blocking rotational DOFs for #"<<id2); state2->blockedDOFs=State::DOF_RXRYRZ;}
+	if(state1->blockedDOFs!=State::DOF_ALL) { LOG_INFO("Blocking all DOFs for #"<<id1); state1->blockedDOFs=State::DOF_ALL;}
+	if(state2->blockedDOFs!=State::DOF_ALL) { LOG_INFO("Blocking all DOFs for #"<<id2); state2->blockedDOFs=State::DOF_ALL;}
 
 	if(step>*(_pathT.rbegin())){
 		LOG_INFO("Last step done, setting zero velocities on #"<<id1<<", #"<<id2<<".");
@@ -104,9 +106,11 @@
 		}
 		// 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;
+		trsf.col(0)=axX; trsf.col(1)=axY; trsf.col(2)=axZ;
 	} else {
 		trsf=l3Geom->trsf;
+		axX=trsf.col(0); axY=trsf.col(1); axZ=trsf.col(2);
+		ptGeom=l3Geom->u;
 	}
 	trsfQ=Quaternionr(trsf);
 	contPt=gsc->contactPoint;

=== modified file 'pkg/dem/L3Geom.cpp'
--- pkg/dem/L3Geom.cpp	2010-10-14 13:58:47 +0000
+++ pkg/dem/L3Geom.cpp	2010-10-14 17:21:34 +0000
@@ -15,8 +15,9 @@
 
 	// contact exists, go ahead
 
-	Real uN=relPos.norm()-(r1+r2);
-	Vector3r normal=relPos.normalized();
+	Real dist=relPos.norm();
+	Real uN=dist-(r1+r2);
+	Vector3r normal=relPos/dist;
 	Vector3r contPt=se31.position+(r1+0.5*uN)*normal;
 
 	// create geometry
@@ -111,6 +112,8 @@
 	Vector3r localF=geom->u.cwise()*Vector3r(phys->kn,phys->ks,phys->ks);
 	// transform back to global coords
 	Vector3r globalF=Quaternionr(geom->trsf).conjugate()*localF;
+	// only for keeping track of the force
+	phys->normalForce=geom->normal*globalF.dot(geom->normal); phys->shearForce=globalF-phys->normalForce;
 	// 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 13:58:47 +0000
+++ pkg/dem/L3Geom.hpp	2010-10-14 17:21:34 +0000
@@ -8,11 +8,12 @@
 #include<yade/pkg-dem/FrictPhys.hpp>
 
 struct L3Geom: public GenericSpheresContact{
+	Real& uN;
+	Vector2r& uT; 
 	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]",
+	YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(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
@@ -20,7 +21,13 @@
 		*/
 		((Matrix3r,trsf,Matrix3r::Identity(),,"Transformation (rotation) from global to local coordinates. (the translation part is in :yref:`GenericSpheresContact.contactPoint`)"))
 		,
+		/*init*/
+		((uN,u[0])) ((uT,*((Vector2r*)&u[1])))
+		,
 		/*ctor*/ createIndex();
+		, /*py*/
+		//.def_readonly("uN",&L3Geom::uN,"Normal component of *u*")
+		//.def_readonly("uT",&L3Geom::uT,"Shear component of *u*")
 	);
 	REGISTER_CLASS_INDEX(L3Geom,GenericSpheresContact);
 };

=== modified file 'scripts/test/law-test.py'
--- scripts/test/law-test.py	2010-10-12 16:46:41 +0000
+++ scripts/test/law-test.py	2010-10-14 17:21:34 +0000
@@ -20,8 +20,8 @@
 # place sphere 1 at the origin
 pt1=Vector3(0,0,0)
 # random orientation of the interaction
-normal=Vector3(random.random()-.5,random.random()-.5,random.random()-.5)
-#normal=Vector3(1,0,0)
+#normal=Vector3(random.random()-.5,random.random()-.5,random.random()-.5)
+normal=Vector3(1,0,0)
 O.bodies.append([
 	utils.sphere(pt1,r1,wire=True,color=(.7,.7,.7)),
 	utils.sphere(pt1+.999999*(r1+r2)*normal.normalized(),r2,wire=True,color=(0,0,0))
@@ -32,12 +32,11 @@
 	PyRunner(iterPeriod=1,command='import time; time.sleep(.01)'),
 	InsertionSortCollider([Bo1_Sphere_Aabb()]),
 	InteractionLoop(
-		[Ig2_Sphere_Sphere_ScGeom()],
-		[Ip2_FrictMat_FrictMat_FrictPhys()],
-		[Law2_ScGeom_FrictPhys_Basic()]
+		[Ig2_Sphere_Sphere_ScGeom()],	[Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_Basic()]
+		# [Ig2_Sphere_Sphere_L3Geom_Inc()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L3Geom_FrictPhys_Linear()] # use this line for L3Geom
 	),
 	LawTester(ids=[0,1],path=[
-		(-.1,0,0),(-.1,.1,0),(-1e-5,.1,0), # towards, shear, back to intial normal distance
+		(-1e-5,0,0),(-.1,0,0),(-.1,.1,0),(-1e-5,.1,0), # towards, shear, back to intial normal distance
 		(-.02,.1,.1),(-.02,-.1,.1),(-.02,-.1,-.1),(-.02,.1,-.1),(-.02,.1,.1), # go in square in the shear plane without changing normal deformation
 		(-1e-4,0,0) # back to the origin, but keep some overlap to not delete the interaction
 		],pathSteps=[50],doneHook='tester.dead=True; O.pause()',label='tester',rotWeight=.2,idWeight=.2),
@@ -60,4 +59,4 @@
 O.dt=1
 
 O.saveTmp()
-O.run()
+#O.run()