yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 1929
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Wed 2009-12-30 18:38:37 +0100
message:
  1. Separate rotation&shear from scaling in cell. Fix possible error (but also possibly not?) error in collider. Spheres' Aabb's are still probably not correctly enlarged.
modified:
  core/Cell.hpp
  core/SConscript
  pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp
  pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp
  pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp
  scripts/test/periodic-shear.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 'core/Cell.hpp'
--- core/Cell.hpp	2009-12-30 15:45:32 +0000
+++ core/Cell.hpp	2009-12-30 17:38:37 +0000
@@ -9,7 +9,15 @@
 * Matrix3r *trsf* is "deformation gradient tensor" F (written as matrix) 
 * Matrix3r *velGrad* is …
 
+The transformation has normal part and rotation/shear part. the shearPt, unshearPt, getShearTrsf etc functions refer to both shear and rotation.
+
 */
+
+#pragma once
+
+#include<yade/lib-serialization/Serializable.hpp>
+#include<yade/lib-base/yadeWm3Extra.hpp>
+
 class Cell: public Serializable{
 	public:
 		Cell(): refSize(Vector3r(1,1,1)), trsf(Matrix3r::IDENTITY), velGrad(Matrix3r::ZERO){ integrateAndUpdate(0); }
@@ -24,11 +32,9 @@
 	const Vector3r& getSize() const { return _size; }
 	//! Return copy of the current size (used only by the python wrapper)
 	Vector3r getSize_copy() const { return _size; }
-	//! stretch ratio Λ(n) (http://en.wikipedia.org/wiki/Finite_strain_theory#Stretch_ratio)
-	Vector3r getStretchRatio() const { return Vector3r(trsf[0][0],trsf[1][1],trsf[2][2]); }
 	//! return vector of consines of skew angle in yz, xz, xy planes between respective transformed base vectors
 	const Vector3r& getCos() const {return _cos;}
-	//! transformation matrix applying puse shear (normal strain removed: ones on the diagonal)
+	//! transformation matrix applying pure shear&rotation (scaling removed)
 	const Matrix3r& getShearTrsf() const { return _shearTrsf; }
 	//! inverse of getShearTrsfMatrix().
 	const Matrix3r& getUnshearTrsf() const {return _unshearTrsf;}
@@ -58,48 +64,11 @@
 		bool _hasShear;
 		Matrix3r _shearTrsf, _unshearTrsf;
 		double _glShearTrsfMatrix[16];
-		void fillGlShearTrsfMatrix(double m[16]){
-			m[0]=trsf[0][0]; m[4]=trsf[0][1]; m[8]=trsf[0][2]; m[12]=0;
-			m[1]=trsf[1][0]; m[5]=trsf[1][1]; m[9]=trsf[1][2]; m[13]=0;
-			m[2]=trsf[2][0]; m[6]=trsf[2][1]; m[10]=trsf[2][2];m[14]=0;
-			m[3]=0;          m[7]=0;          m[11]=0;         m[15]=1;
-		}
-
+		void fillGlShearTrsfMatrix(double m[16]);
 	public:
 
 	//! "integrate" velGrad, update cached values used by public getter
-	void integrateAndUpdate(Real dt){
-		//incremental displacement gradient
-		_trsfInc=dt*velGrad;
-		// total transformation; M = (Id+G).M = F.M
-		trsf+=_trsfInc*trsf;
-		// Hsize will contain colums with transformed base vectors
-		Matrix3r Hsize(refSize[0],refSize[1],refSize[2]);
-		Hsize=trsf*Hsize;
-		// lengths of transformed cell vectors, skew cosines
-		Vector3r Hnorm[3]; // normalized transformed base vectors
-		for(int i=0; i<3; i++){
-			_size[i]=Vector3r(Hsize[i][0],Hsize[i][1],Hsize[i][2]).Length();
-			Hnorm[i]=Vector3r(Hsize[i][0],Hsize[i][1],Hsize[i][2]); Hnorm[i].Normalize();
-		};
-		//cerr<<"Hsize="<<endl<<Hsize[0][0]<<" "<<Hsize[0][1]<<" "<<Hsize[0][2]<<endl<<Hsize[1][0]<<" "<<Hsize[1][1]<<" "<<Hsize[1][2]<<endl<<Hsize[2][0]<<" "<<Hsize[2][1]<<" "<<Hsize[2][2]<<endl;
-		//cerr<<"Hnorm="<<Hnorm[0]<<" | "<<Hnorm[1]<<" | "<<Hnorm[2]<<endl;
-		// skew cosines
-		for(int i=0; i<3; i++){
-			int i1=(i+1)%3, i2=(i+2)%3;
-			// sin between axes is cos of skew
-			_cos[i]=(Hnorm[i1].Cross(Hnorm[i2])).SquaredLength();
-		}
-		// pure shear trsf: ones on diagonal
-		_shearTrsf=trsf; _shearTrsf[0][0]=_shearTrsf[1][1]=_shearTrsf[2][2]=1.;
-		// pure unshear transformation
-		_unshearTrsf=_shearTrsf.Inverse();
-		// some parts can branch depending on presence/absence of shear
-		_hasShear=(trsf[0][1]!=0 || trsf[0][2]!=0 || trsf[1][0]!=0 || trsf[1][2]!=0 || trsf[2][0]!=0 || trsf[2][1]!=0);
-		// OpenGL shear matrix (used frequently)
-		fillGlShearTrsfMatrix(_glShearTrsfMatrix);
-	}
-	
+	void integrateAndUpdate(Real dt);
 	/*! Return point inside periodic cell, even if shear is applied */
 	Vector3r wrapShearedPt(const Vector3r& pt){ return shearPt(wrapPt(unshearPt(pt))); }
 	/*! Return point inside periodic cell, even if shear is applied; store cell coordinates in period. */

=== modified file 'core/SConscript'
--- core/SConscript	2009-12-09 17:11:51 +0000
+++ core/SConscript	2009-12-30 17:38:37 +0000
@@ -23,6 +23,7 @@
 			'Body.cpp',
 			'BodyContainer.cpp',
 			'Bound.cpp',
+			'Cell.cpp',
 			'Collider.cpp',
 			'PartialEngine.cpp',
 			'FileGenerator.cpp',

=== modified file 'pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp'
--- pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp	2009-12-30 15:45:32 +0000
+++ pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp	2009-12-30 17:38:37 +0000
@@ -14,8 +14,12 @@
 	Sphere* sphere = static_cast<Sphere*>(cm.get());
 	Aabb* aabb = static_cast<Aabb*>(bv.get());
 	Vector3r halfSize = (aabbEnlargeFactor>0?aabbEnlargeFactor:1.)*Vector3r(sphere->radius,sphere->radius,sphere->radius);
+	if(!scene->isPeriodic){
+		aabb->min=se3.position-halfSize; aabb->max=se3.position+halfSize;
+		return;
+	}
 	// adjust box size along axes so that sphere doesn't stick out of the box even if sheared (i.e. parallelepiped)
-	if(scene->isPeriodic && scene->cell->hasShear()) {
+	if(scene->cell->hasShear()) {
 		Vector3r refHalfSize(halfSize);
 		const Vector3r& cos=scene->cell->getCos();
 		for(int i=0; i<3; i++){

=== modified file 'pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp'
--- pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp	2009-12-25 14:46:48 +0000
+++ pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp	2009-12-30 17:38:37 +0000
@@ -403,14 +403,15 @@
 	assert(id1!=id2); // programming error, or weird bodies (too large?)
 	for(int axis=0; axis<3; axis++){
 		Real dim=scene->cell->getSize()[axis];
+		// LOG_DEBUG("dim["<<axis<<"]="<<dim);
 		// too big bodies in interaction
 		assert(maxima[3*id1+axis]-minima[3*id1+axis]<.99*dim); assert(maxima[3*id2+axis]-minima[3*id2+axis]<.99*dim);
-		// find body of which when taken as period start will make the gap smaller
+		// find body of which minimum when taken as period start will make the gap smaller
 		Real m1=minima[3*id1+axis],m2=minima[3*id2+axis];
 		Real wMn=(cellWrapRel(m1,m2,m2+dim)<cellWrapRel(m2,m1,m1+dim)) ? m2 : m1;
 		#ifdef PISC_DEBUG
 		if(watchIds(id1,id2)){
-			TRVAR3(id1,id2,axis);
+			TRVAR4(id1,id2,axis,dim);
 			TRVAR4(minima[3*id1+axis],maxima[3*id1+axis],minima[3*id2+axis],maxima[3*id2+axis]);
 			TRVAR2(cellWrapRel(m1,m2,m2+dim),cellWrapRel(m2,m1,m1+dim));
 			TRVAR3(m1,m2,wMn);
@@ -425,7 +426,10 @@
 				TRVAR4(pmn1,pmx1,pmn2,pmx2);
 			}
 		#endif
-		if((pmn1!=pmx1) || (pmn2!=pmx2)){
+		/* FIXME: condition temporary disabled
+			if wMn==m1 and minima[id2]<m1, then it might give spurious error; why should it be error, though?
+		*/
+		if(false && ((pmn1!=pmx1) || (pmn2!=pmx2))){
 			LOG_FATAL("Body #"<<(pmn1!=pmx1?id1:id2)<<" spans over half of the cell size "<<dim<<" (axis="<<axis<<", min="<<(pmn1!=pmx1?mn1:mn2)<<", max="<<(pmn1!=mn1?mx1:mx2)<<")");
 			throw runtime_error(__FILE__ ": Body larger than half of the cell size encountered.");
 		}
@@ -440,6 +444,6 @@
 
 #ifdef PISC_DEBUG
 	bool InsertionSortCollider::watchIds(body_id_t id1, body_id_t id2) const{
-		return true; //id1==1 || id2==1;
+		return id1==3 || id2==3; //true; //id1==1 || id2==1;
 	}
 #endif

=== modified file 'pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp'
--- pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp	2009-12-30 15:45:32 +0000
+++ pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp	2009-12-30 17:38:37 +0000
@@ -138,12 +138,11 @@
 	if(!scene->isPeriodic) return;
 	glColor3v(Vector3r(1,1,0));
 	glPushMatrix();
-		Vector3r refSize=scene->cell->refSize;
-		if(scaleDisplacements) refSize=diagMult(displacementScale,refSize);
-		//glTranslatev(scene->cell->shearPt(.5*size)); // shear center (moves when sheared)
-		glTranslatev(scene->cell->trsf*(.5*refSize));
+		Vector3r size=scene->cell->getSize();
+		if(scaleDisplacements) size=diagMult(displacementScale,size);
+		glTranslatev(scene->cell->shearPt(.5*size)); // shear center (moves when sheared)
 		glMultMatrixd(scene->cell->getGlShearTrsfMatrix());
-		glScalev(refSize);
+		glScalev(size);
 		glutWireCube(1);
 	glPopMatrix();
 }

=== modified file 'scripts/test/periodic-shear.py'
--- scripts/test/periodic-shear.py	2009-12-26 21:57:37 +0000
+++ scripts/test/periodic-shear.py	2009-12-30 17:38:37 +0000
@@ -40,6 +40,6 @@
 v.axes=True
 v.grid=(True,True,True)
 
-#from yade import log
+from yade import log
 #log.setLevel('Shop',log.TRACE)
 #log.setLevel('InsertionSortCollider',log.TRACE)