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;}
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
