yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 1916
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Tue 2009-12-22 20:01:20 +0100
message:
  Code using the velocity gradient to define periodic deformations (experimental) + script to test it.
  Needs compilation with #defined VELGRAD (currently needs uncommenting "#define VELGRAD" in 4 files in total). 
added:
  scripts/test/periodic-triax-velgrad.py
modified:
  core/Cell.hpp
  core/Scene.cpp
  core/Scene.hpp
  pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp
  pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp
  py/yadeWrapper/yadeWrapper.cpp


--
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 'core/Cell.hpp'
--- core/Cell.hpp	2009-12-20 22:03:17 +0000
+++ core/Cell.hpp	2009-12-22 19:01:20 +0000
@@ -16,6 +16,9 @@
 #endif
 // end yade compatibility
 
+	
+//#define VELGRAD
+	
 /*! Periodic cell parameters and routines. Usually instantiated as Scene::cell.
 
 The cell has size and shear, which are independent. 
@@ -23,7 +26,16 @@
 */
 class Cell: public Serializable{
 	public:
-		Cell(): refSize(Vector3r(1,1,1)), strain(Matrix3r::ZERO){ updateCache(); }
+		Cell(): refSize(Vector3r(1,1,1)), strain(Matrix3r::ZERO)
+#ifdef VELGRAD
+				, velGrad(Matrix3r::ZERO) 
+				, Hsize(Matrix3r::IDENTITY)
+				,_shearTrsfMatrix(Matrix3r::IDENTITY)
+				{ updateCache(0); }
+#else
+				{ updateCache(); }
+#endif
+				
 		//! size of the cell, projected on axes (for non-zero shear)
 		Vector3r refSize;
 		//! 3 non-diagonal components of the shear matrix, ordered by axis normal to the shear plane,
@@ -33,6 +45,11 @@
 		//! for details
 		// Vector3r shear;
 		Matrix3r strain;
+#ifdef VELGRAD
+		Matrix3r velGrad;
+		Matrix3r Hsize;
+		Matrix3r _shearIncrt;
+#endif
 
 		//! reference values of size and shear (for rendering, mainly)
 		// Vector3r refShear, refCenter;
@@ -49,6 +66,11 @@
 	const Matrix3r& getShearTrsfMatrix() const { return _shearTrsfMatrix; }
 	//! inverse of getShearTrsfMatrix().
 	const Matrix3r& getUnshearTrsfMatrix() const {return _unshearTrsfMatrix;}
+#ifdef VELGRAD	
+	//! transformation increment matrix applying arbitrary field
+	const Matrix3r& getShearIncrMatrix() const { return _shearIncrt; }
+#endif	
+	
 	/*! return pointer to column-major OpenGL 4x4 matrix applying pure shear. This matrix is suitable as argument for glMultMatrixd.
 
 	Note: the order of OpenGL transoformations matters; for instance, if you draw sheared wire box of size *size*,
@@ -83,6 +105,21 @@
 	public:
 
 	// should be called before every step
+#ifdef VELGRAD
+	void updateCache(double dt){	
+		//incremental disp gradient
+		_shearIncrt=dt*velGrad;
+		
+		//update Hsize (redundant with total transformation perhaps)
+		Hsize=Hsize+_shearIncrt*Hsize;
+		//total transformation
+		_shearTrsfMatrix = _shearTrsfMatrix + _shearIncrt*_shearTrsfMatrix;// M = (Id+G).M = F.M
+		//compatibility with Vaclav's code :
+		_size[0]=Hsize[0][0]; _size[1]=Hsize[1][1]; _size[2]=Hsize[2][2];
+		_hasShear = false;
+		_unshearTrsfMatrix=_shearTrsfMatrix.Inverse();
+		strain=_shearTrsfMatrix;
+#else
 	void updateCache(){
 		/*
 		for(int i=0; i<3; i++) {
@@ -98,10 +135,14 @@
 		_unshearTrsfMatrix=_shearTrsfMatrix.Inverse();
 		// _shearTrsf=Matrix3r(1,shear[2],shear[1],shear[2],1,shear[0],shear[1],shear[0],1);
 		// _unshearTrsf=strain.Inverse();
+#endif
 		fillGlShearTrsfMatrix(_glShearTrsfMatrix);
 		for(int i=0; i<3; i++) for(int j=0; j<3; j++) _cosMatrix[i][j]=(i==j?0:cos(atan(strain[i][j])));
 	}
-
+	
+// #ifdef VELGRAD		
+// 	void updateCache(){ updateCache(0);}
+// #endif
 
 	// doesn't seem to be really useful
 	#if 0
@@ -144,9 +185,16 @@
 	static Real wrapNum(const Real& x, const Real& sz, int& period){
 		Real norm=x/sz; period=(int)floor(norm); return (norm-period)*sz;
 	}
-
+#ifdef VELGRAD
+	void postProcessAttributes(bool deserializing){ if(deserializing) updateCache(0); }
+#else
 	void postProcessAttributes(bool deserializing){ if(deserializing) updateCache(); }
-	REGISTER_ATTRIBUTES(Serializable,(refSize)(strain));
+#endif
+	REGISTER_ATTRIBUTES(Serializable,(refSize)(strain)
+#ifdef VELGRAD
+			(velGrad)(Hsize)
+#endif
+			   );
 	REGISTER_CLASS_AND_BASE(Cell,Serializable);
 };
 REGISTER_SERIALIZABLE(Cell);

=== modified file 'core/Scene.cpp'
--- core/Scene.cpp	2009-12-21 22:19:11 +0000
+++ core/Scene.cpp	2009-12-22 19:01:20 +0000
@@ -81,7 +81,11 @@
 		needsInitializers=false;
 	}
 	// update cell data
+#ifndef VELGRAD
 	if(isPeriodic) cell->updateCache();
+#else
+	if(isPeriodic) cell->updateCache(dt);
+#endif
 	//bex.reset(); // uncomment if PhysicalActionContainerReseter is removed
 	bool TimingInfo_enabled=TimingInfo::enabled; // cache the value, so that when it is changed inside the step, the engine that was just running doesn't get bogus values
 	TimingInfo::delta last=TimingInfo::getNow(); // actually does something only if TimingInfo::enabled, no need to put the condition here

=== modified file 'core/Scene.hpp'
--- core/Scene.hpp	2009-12-20 22:03:17 +0000
+++ core/Scene.hpp	2009-12-22 19:01:20 +0000
@@ -23,6 +23,8 @@
 #define HOST_NAME_MAX 255 
 #endif
 
+//#define VELGRAD
+
 class Bound;
 
 class Scene: public Serializable{

=== modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp'
--- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp	2009-12-20 22:03:17 +0000
+++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp	2009-12-22 19:01:20 +0000
@@ -12,6 +12,9 @@
 #include<yade/pkg-common/VelocityBins.hpp>
 #include<yade/lib-base/yadeWm3Extra.hpp>
 
+//#define VELGRAD
+
+
 YADE_PLUGIN((NewtonIntegrator));
 CREATE_LOGGER(NewtonIntegrator);
 void NewtonIntegrator::cundallDamp(const Real& dt, const Vector3r& N, const Vector3r& V, Vector3r& A){
@@ -62,7 +65,9 @@
 void NewtonIntegrator::action(Scene*)
 {
 	// only temporary
+#ifndef VELGRAD
 	if(homotheticCellResize) throw runtime_error("Homothetic resizing not yet implemented with new, core/Cell.hpp based cell (extension+shear).");
+#endif
 
 	scene->bex.sync();
 	Real dt=scene->dt;
@@ -166,10 +171,20 @@
 	//   or position by (x-x') (homotheticCellResize==2)
 	// FIXME: wrap state->pos first, then apply the shift; otherwise result might be garbage
 	Vector3r dPos(Vector3r::ZERO); // initialize to avoid warning; find way to avoid it in a better way
+#ifndef VELGRAD
 	if(cellChanged && homotheticCellResize){ for(int i=0; i<3; i++) dPos[i]=(state->pos[i]/prevCellSize[i])*scene->cell->getSize()[i]-state->pos[i]; }
+#else
+	if(homotheticCellResize){ dPos = scene->cell->getShearIncrMatrix()*scene->cell->wrapShearedPt(state->pos);
+	}
+#endif	
 	blockTranslateDOFs(state->blockedDOFs, state->accel);
-	state->vel+=dt*state->accel;                             if(homotheticCellResize==1) state->vel+=dPos/dt;
-	state->pos += state->vel*dt + scene->bex.getMove(id);    if(homotheticCellResize==2) state->pos+=dPos;
+	state->vel+=dt*state->accel;
+	state->pos += state->vel*dt + scene->bex.getMove(id);
+	//apply cell deformation
+	if(homotheticCellResize>=1) state->pos+=dPos;
+	//update velocity for usage in rate dependant equations (e.g. viscous law) FIXME : it is not recommended to do that because it impacts the dynamics (this modified velocity will be used as reference in the next time-step)
+	if(homotheticCellResize==2) state->vel+=dPos/dt;
+
 }
 inline void NewtonIntegrator::leapfrogSphericalRotate(Scene* scene, State* state, const body_id_t& id, const Real& dt )
 {

=== modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp'
--- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp	2009-12-15 09:41:07 +0000
+++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp	2009-12-22 19:01:20 +0000
@@ -73,7 +73,7 @@
 		/// velocity bins (not used if not created)
 		shared_ptr<VelocityBins> velocityBins;
 		virtual void action(Scene *);		
-		NewtonIntegrator(): prevCellSize(Vector3r::ZERO),damping(0.2), maxVelocitySq(-1), exactAsphericalRot(false), homotheticCellResize(0){
+		NewtonIntegrator(): prevCellSize(Vector3r::ZERO),damping(0.2), maxVelocitySq(-1), exactAsphericalRot(false), homotheticCellResize(1){
 			#ifdef YADE_OPENMP
 				threadMaxVelocitySq.resize(omp_get_max_threads());
 			#endif

=== modified file 'py/yadeWrapper/yadeWrapper.cpp'
--- py/yadeWrapper/yadeWrapper.cpp	2009-12-20 22:03:17 +0000
+++ py/yadeWrapper/yadeWrapper.cpp	2009-12-22 19:01:20 +0000
@@ -56,6 +56,8 @@
 using namespace boost;
 using namespace std;
 
+//#define VELGRAD
+
 #include<yade/extra/boost_python_len.hpp>
 
 
@@ -803,6 +805,10 @@
 	python::class_<Cell,shared_ptr<Cell>, python::bases<Serializable>, noncopyable>("Cell",python::no_init)
 		.def_readwrite("refSize",&Cell::refSize)
 		.def_readwrite("strain",&Cell::strain)
+#ifdef VELGRAD
+		.def_readwrite("velGrad",&Cell::velGrad)
+		.def_readwrite("Hsize",&Cell::Hsize)
+#endif
 		.add_property("extension",&Cell::getExtensionalStrain)
 		//.add_property("size",&Cell::getSize,python::return_value_policy<python::return_internal_referece>()
 	;

=== added file 'scripts/test/periodic-triax-velgrad.py'
--- scripts/test/periodic-triax-velgrad.py	1970-01-01 00:00:00 +0000
+++ scripts/test/periodic-triax-velgrad.py	2009-12-22 19:01:20 +0000
@@ -0,0 +1,50 @@
+# coding: utf-8
+# 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+"Test and demonstrate use of PeriTriaxController."
+from yade import *
+from yade import pack,log,qt
+#log.setLevel('PeriTriaxController',log.DEBUG)
+O.periodic=True
+O.cell.refSize=Vector3(.1,.1,.1)
+O.cell.Hsize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
+sp=pack.SpherePack()
+radius=5e-3
+num=sp.makeCloud(Vector3().ZERO,O.cell.refSize,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
+O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])
+
+
+O.engines=[
+	BexResetter(),
+	BoundDispatcher([Bo1_Sphere_Aabb()]),
+	InsertionSortCollider(nBins=5,sweepLength=.05*radius),
+	InteractionDispatchers(
+		[Ig2_Sphere_Sphere_Dem3DofGeom()],
+		[SimpleElasticRelationships()],
+		[Law2_Dem3Dof_Elastic_Elastic()]
+	),
+	PeriTriaxController(goal=[-1e5,-1e5,0],stressMask=3,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
+	NewtonIntegrator(damping=.6),
+]
+O.dt=0.1*utils.PWaveTimeStep()
+O.run(1)
+qt.View()
+O.cell.velGrad=Matrix3(0,5,0,0,0,0, 0,0,-5)
+O.run();
+
+
+phase=0
+def triaxDone():
+	global phase
+	if phase==0:
+		print 'Here we are: stress',triax['stress'],'strain',triax['strain'],'stiffness',triax['stiff']
+		print 'Now εz will go from 0 to .2 while σx and σy will be kept the same.'
+		triax['goal']=[-1e5,-1e5,.2]
+		O.cell.velGrad=Matrix3(0,0,0,5,0,0, 0,0,0)
+		phase+=1
+	elif phase==1:
+		print 'Here we are: stress',triax['stress'],'strain',triax['strain'],'stiffness',triax['stiff']
+		print 'Done, pausing now.'
+		O.pause()
+		
+	
+