yade-dev team mailing list archive

[Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>]

> OK, I see what you mean. You implemented what was already working,
> sorry. (I was rotating the whole cell around about a week ago, with
> matrix that was antisymmetric on non-diagonal.) You could've tried
> before that it worked...
>
>   

I'm not very sure what you mean. What can I try?
An example of working collider in refSize is attached (diff vs. 1935).

Bruno




--

_______________
Chareyre Bruno
Maître de Conférences

Grenoble INP
Laboratoire 3SR - bureau E145
BP 53 - 38041, Grenoble cedex 9 - France
Tél : 33 4 56 52 86 21
Fax : 33 4 76 82 70 43
________________

=== modified file 'core/Cell.cpp'
--- core/Cell.cpp	2009-12-30 17:39:02 +0000
+++ core/Cell.cpp	2010-01-05 00:24:20 +0000
@@ -24,13 +24,66 @@
 	}
 	// pure shear trsf: ones on diagonal
 	_shearTrsf=Hnorm;
 	// pure unshear transformation
 	_unshearTrsf=_shearTrsf.Inverse();
+	_invTrsf=trsf.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);
+	Matrix3r U2 (trsf.Transpose()*trsf);
+	Matrix3r R,U,Q;
+	U2.EigenDecomposition(Q,U);
+			
+	//Be carefull, wm3 returns nan's sometimes when we factorize Identity... In that case, do nothing anyway and keep unit aabb.
+	if (trsf!=Matrix3r::IDENTITY) {
+		for (int i=0; i<3; i++) U(i,i)=sqrt(U(i,i));		
+		U=Q*U*(Q.Transpose());
+		R=trsf*(U.Inverse());
+		//Decompose U again into rotation and stretch? no use for now
+		//Matrix3r U3,Q3;
+		//U.EigenDecomposition(Q3,U3);
+		Matrix3r Unorm; // normalized transformed base vectors
+		for(int i=0; i<3; i++){
+			Vector3r base(U.GetColumn(i));
+			base.Normalize(); //returns length
+			Unorm.SetColumn(i,base);};
+		Matrix3r Uortho;
+		for(int i=0; i<3; i++){
+			int i1=(i+1)%3, i2=(i+2)%3;
+			Vector3r base(U.GetColumn(i1).Cross(U.GetColumn(i2)));
+			base.Normalize(); //returns length
+			Uortho.SetColumn(i,base);};
+		
+		//Define half-size of U-transformed aabb
+		_unitTrsfHalfSize=Vector3r::ZERO;
+		for (int i=0;i<3;i++) _unitTrsfHalfSize+=Unorm.GetColumn(i)/(Unorm.GetColumn(i).Dot(Uortho.GetColumn(i)));
+		_unitTrsfHalfSize=U.Inverse()*_unitTrsfHalfSize;
+	}
+	else _unitTrsfHalfSize=Vector3r(1,1,1);
 }
 
 void Cell::fillGlShearTrsfMatrix(double m[16]){

=== modified file 'core/Cell.hpp'
--- core/Cell.hpp	2009-12-30 17:38:37 +0000
+++ core/Cell.hpp	2010-01-03 22:12:03 +0000
@@ -38,6 +38,10 @@
 	const Matrix3r& getShearTrsf() const { return _shearTrsf; }
 	//! inverse of getShearTrsfMatrix().
 	const Matrix3r& getUnshearTrsf() const {return _unshearTrsf;}
+	//! return vector halfSize of a transformed Aabb (parallelepiped) with inscribed sphere of radius 1, so that the untransformed Aabb is axis aligned in reference frame
+	const Vector3r& getUnitTrsfHalfSize() const {return _unitTrsfHalfSize;}
 	//! transformation increment matrix applying arbitrary field (remove if not used in NewtonIntegrator!)
 	// const Matrix3r& getTrsfInc() const { return _trsfInc; }
 	
@@ -60,30 +64,37 @@
 	// caches; private
 	private:
 		Matrix3r _trsfInc;
-		Vector3r _size, _cos;
+		Vector3r _size, _cos, _unitTrsfHalfSize;
 		bool _hasShear;
-		Matrix3r _shearTrsf, _unshearTrsf;
+		Matrix3r _shearTrsf, _unshearTrsf, _invTrsf;
 		double _glShearTrsfMatrix[16];
 		void fillGlShearTrsfMatrix(double m[16]);
+		
 	public:
 
 	//! "integrate" velGrad, update cached values used by public getter
 	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 for arbitrary transformation applied */
+	Vector3r wrapTrsfdPt(const Vector3r& pt){ return trsfPt(wrapPt(unTrsfPt(pt))); }
 	/*! Return point inside periodic cell, even if shear is applied; store cell coordinates in period. */
 	Vector3r wrapShearedPt(const Vector3r& pt, Vector3<int>& period){ return shearPt(wrapPt(unshearPt(pt),period)); }
 	/*! Apply inverse shear on point; to put it inside (unsheared) periodic cell, apply wrapPt on the returned value. */
 	Vector3r unshearPt(const Vector3r& pt){ return _unshearTrsf*pt; }
 	//! Apply shear on point. 
 	Vector3r shearPt(const Vector3r& pt){ return _shearTrsf*pt; }
+	/*! Apply inverse transformation on point; to put it inside (initial) periodic cell, apply wrapPt on the returned value. */
+	Vector3r unTrsfPt(const Vector3r& pt){ return _invTrsf*pt; }
+	//! Apply transformation on point. 
+	Vector3r trsfPt(const Vector3r& pt){ return trsf*pt; }
 	/*! Wrap point to inside the periodic cell; don't compute number of periods wrapped */
 	Vector3r wrapPt(const Vector3r pt)const{
 		Vector3r ret; for(int i=0;i<3;i++) ret[i]=wrapNum(pt[i],_size[i]); return ret;
 	}
 	/*! Wrap point to inside the periodic cell; period will contain by how many cells it was wrapped. */
 	Vector3r wrapPt(const Vector3r pt, Vector3<int>& period)const{
-		Vector3r ret; for(int i=0; i<3; i++){ ret[i]=wrapNum(pt[i],_size[i],period[i]); } return ret;
+		Vector3r ret; for(int i=0; i<3; i++){ ret[i]=wrapNum(pt[i],refSize[i],period[i]); } return ret;
 	}
 	/*! Wrap number to interval 0…sz */
 	static Real wrapNum(const Real& x, const Real& sz){

=== modified file 'pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp'
--- pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp	2009-12-25 14:46:48 +0000
+++ pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp	2010-01-02 22:03:02 +0000
@@ -91,6 +91,11 @@
 				Vector3r shift2(I->cellDist[0]*cellSize[0],I->cellDist[1]*cellSize[1],I->cellDist[2]*cellSize[2]);
 				// in sheared cell, apply shear on the mutual position as well
 				shift2=scene->cell->shearPt(shift2);
+				// Hsize will contain colums with transformed base vectors
+				Matrix3r Hsize(scene->cell->refSize[0],scene->cell->refSize[1],scene->cell->refSize[2]); Hsize=scene->cell->trsf*Hsize;
+				Vector3r shift3((Real) I->cellDist[0]*Hsize.GetColumn(0)+(Real)I->cellDist[1]*Hsize.GetColumn(1)+(Real)I->cellDist[2]*Hsize.GetColumn(2));
+				if ((Omega::instance().getCurrentIteration() % 100 == 0)) LOG_DEBUG(shift2 << " vs " << shift3);
+				
 				geomCreated=I->functorCache.geom->go(b1->shape,b2->shape,*b1->state,*b2->state,shift2,/*force*/false,I);
 			}
 			if(!geomCreated){

=== modified file 'pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp'
--- pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp	2009-12-30 20:10:59 +0000
+++ pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp	2010-01-04 23:36:11 +0000
@@ -10,6 +10,9 @@
 #include<yade/pkg-common/Sphere.hpp>
 #include<yade/pkg-common/Aabb.hpp>
 
+CREATE_LOGGER(Bo1_Sphere_Aabb);
+
+
 void Bo1_Sphere_Aabb::go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r& se3, const Body* b){
 	Sphere* sphere = static_cast<Sphere*>(cm.get());
 	Aabb* aabb = static_cast<Aabb*>(bv.get());
@@ -19,19 +22,59 @@
 		return;
 	}
 	// adjust box size along axes so that sphere doesn't stick out of the box even if sheared (i.e. parallelepiped)
-	if(scene->cell->hasShear()) {
-		Vector3r refHalfSize(halfSize);
-		const Vector3r& cos=scene->cell->getCos();
-		for(int i=0; i<3; i++){
-			//cerr<<"cos["<<i<<"]"<<cos[i]<<" ";
-			int i1=(i+1)%3,i2=(i+2)%3;
-			halfSize[i1]+=.5*refHalfSize[i1]*(1/cos[i]-1);
-			halfSize[i2]+=.5*refHalfSize[i2]*(1/cos[i]-1);
-		}
-	}
-	//cerr<<" || "<<halfSize<<endl;
-	aabb->min = scene->cell->unshearPt(se3.position)-halfSize;
-	aabb->max = scene->cell->unshearPt(se3.position)+halfSize;	


+	Vector3r halfSize2 = (aabbEnlargeFactor>0?aabbEnlargeFactor:1.)*sphere->radius*scene->cell->getUnitTrsfHalfSize();
+	Vector3r min2 = scene->cell->unTrsfPt(se3.position)-halfSize2;
+	Vector3r max2 = scene->cell->unTrsfPt(se3.position)+halfSize2;
+	aabb->min = min2;
+	aabb->max = max2;
 }
 	
 YADE_PLUGIN((Bo1_Sphere_Aabb));


=== modified file 'pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp'
--- pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp	2009-12-30 20:10:59 +0000
+++ pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp	2010-01-04 18:18:40 +0000
@@ -402,7 +402,8 @@
 	assert(periodic);
 	assert(id1!=id2); // programming error, or weird bodies (too large?)
 	for(int axis=0; axis<3; axis++){
-		Real dim=scene->cell->getSize()[axis];
+//		Real dim=scene->cell->getSize()[axis];
+ 		Real dim=scene->cell->refSize[axis];//TEST
 		// 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);

=== modified file 'pkg/common/Engine/GlobalEngine/InsertionSortCollider.hpp'
--- pkg/common/Engine/GlobalEngine/InsertionSortCollider.hpp	2009-12-25 14:46:48 +0000
+++ pkg/common/Engine/GlobalEngine/InsertionSortCollider.hpp	2010-01-04 18:18:40 +0000
@@ -160,6 +160,7 @@
 			assert(scene->isPeriodic);
 			assert(axis>=0 && axis <=2);
 			cellDim=scene->cell->getSize()[axis];
+ 			cellDim=scene->cell->refSize[axis];//TEST
 		}
 		// normalize given index to the right range (wraps around)
 		long norm(long i) const { if(i<0) i+=size; long ret=i%size; assert(ret>=0 && ret<size); return ret;}